diff --git a/CHANGES b/CHANGES index 7f220b5..7129b7e 100644 --- a/CHANGES +++ b/CHANGES @@ -1,6 +1,17 @@ PySurfer Changes ================ +Version 0.4 +------------- + +Enhancements +~~~~~~~~~~~~ + +- Display data from both hemispheres simultaneously +- Display multiple views simultaneously +- Toggling Mayavi toolbars +- Use nibabel for IO functions + Version 0.3.1 ------------- diff --git a/bin/pysurfer b/bin/pysurfer index b9237ff..c602311 100755 --- a/bin/pysurfer +++ b/bin/pysurfer @@ -24,10 +24,16 @@ if __name__ == '__main__': # boot up mlab/IPython if len(sys.argv) > 3: subjects_dir = os.environ['SUBJECTS_DIR'] - surf_file = os.path.join(subjects_dir, - "%s/surf/%s.%s" % tuple(sys.argv[1:4])) - if not os.path.exists(surf_file): - sys.exit("ERROR: Could not find %s" % surf_file) + if sys.argv[2] in ['both', 'split']: + hemi_checks = ['lh', 'rh'] + else: + hemi_checks = [sys.argv[2]] + for h in hemi_checks: + surf_file = os.path.join(subjects_dir, + "%s/surf/%s.%s" % (sys.argv[1], h, + sys.argv[3])) + if not os.path.exists(surf_file): + sys.exit("ERROR: Could not find %s" % surf_file) if not is_ipython: # Parse the args so that --help exits back to the shell @@ -62,10 +68,11 @@ if __name__ == '__main__': argdict = args.__dict__ config_opts = dict([(k, v) for k, v in argdict.items() if k in confkeys and v]) + views = args.views # Load up the figure and underlying brain object b = Brain(args.subject_id, args.hemi, args.surf, args.curv, - args.title, config_opts=config_opts) + args.title, views=views, config_opts=config_opts) # Maybe load some morphometry if args.morphometry is not None: diff --git a/doc/_static/split_view.png b/doc/_static/split_view.png new file mode 100644 index 0000000..3aefb00 Binary files /dev/null and b/doc/_static/split_view.png differ diff --git a/doc/conf.py b/doc/conf.py index 274af03..43b7af1 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -45,7 +45,8 @@ # General information about the project. project = u'PySurfer' -copyright = u'2011-2012, Michael Waskom, Alexandre Gramfort, Scott Burns, Martin Luessi' +copyright = (u'2011-2013, Michael Waskom, Alexandre Gramfort, Scott Burns, ' + 'Martin Luessi, Eric Larson') # Generate the plots for the gallery plot_gallery = True @@ -222,5 +223,5 @@ # (source start file, name, description, authors, manual section). man_pages = [ ('index', 'pysurfer', u'PySurfer Documentation', - [u'Michael Waskom, Alexandre Gramfort, Scott Burns'], 1) + [u'Michael Waskom, Alexandre Gramfort, Scott Burns, Eric Larson'], 1) ] diff --git a/doc/documentation/config_file.rst b/doc/documentation/config_file.rst index 8a41420..871af8d 100644 --- a/doc/documentation/config_file.rst +++ b/doc/documentation/config_file.rst @@ -7,7 +7,7 @@ There are several aspects of how PySurfer looks and behaves that you may wish to customize but do not want to have to alter every time you instantiate a Brain object or use the command-line interface. To facilitate this, PySurfer allows you to set certain options with a -standard Python config file. +standard Python config file. When a new Brain object is created (either in a script or via the command-line), it can read configuration options from one of two places: @@ -27,9 +27,14 @@ Visual *size* How large the sides of the display window should be (measured in - pixels.) The window is always square, so just give one value, and - it will be used for the height and width. (Possible values: any - positive number.) + pixels.) When using this option, the window will be square, and + "size" it will be used for the height and width. (Possible values: + any positive number.) + +*width* and *height* + The size of the window. This is an alternative to size, allowing + the width and height to be set independently. (Possible values: + any positive number.) *default_view* Which view should be shown at the beginning of a visualization @@ -49,3 +54,8 @@ Overlay provided when calling the :meth:`Brain.add_overlay` method. (Possible values: any float, ``robust_max``, ``actual_max``.) +Options +------- +*logging_level* + Amount of verbosity to use in outputting status messages. + (Possible values: ``DEBUG``, ``INFO``, ``WARNING``, ``ERROR``.) diff --git a/doc/documentation/index.rst b/doc/documentation/index.rst index 01f90f0..5de1ba2 100644 --- a/doc/documentation/index.rst +++ b/doc/documentation/index.rst @@ -7,5 +7,6 @@ Documentation ./command_line ./custom_viz + ./split_brain ./config_file diff --git a/doc/documentation/split_brain.rst b/doc/documentation/split_brain.rst new file mode 100644 index 0000000..71cf0c5 --- /dev/null +++ b/doc/documentation/split_brain.rst @@ -0,0 +1,48 @@ +.. _split_brain: + +Working with a split-screen brain +================================= + +The split-screen view can be activated by using the argument ``hemi='split'``. +Using this option will put views of the left hemisphere in consecutive +vertical frames on the left, and views of the right hemisphere in +consecutive vertical frames on the right. For example, running the following:: + + brain = Brain('fsaverage', 'split', 'inflated', views=['lat', 'med']) + +Will produce a window with two columns (hemispheres) and two rows (the +lateral and medial views, respectively), shown below. + +.. image:: ../../_static/split_view.png + +Adding and displaying data +-------------------------- + +Data can be added to either hemisphere using the same functions that are +normally used, e.g. ``add_data``, ``add_overlay``, ``add_morphometry``. +The objects are automatically shown on all views of the brain. When +calling these functions, the ``hemi`` keyword argument can be set to +``hemi='lh'`` or ``hemi='rh'`` to specify the hemisphere to plot to. +In some instances (e.g., ``add_morphometry``), if no keyword argument +is provided, PySurfer will attempt to load data or both hemispheres +automtically. + +Note that the ``show_view`` method accepts arguments for the ``row`` and +``col`` values, which allow the user to control which ``Brain`` panel +gets the updated view. + +Caveats +------- +The multi-view support is available thanks to the capabilities of the +TraitsUI framework. However, due to some limitations in the implementation +of TraitsUI, there is no guarantee that a set of scripted commands will +result in a painted window when the user may think it will. For +example, making a series of calls to ``brain.add_label()`` followed by +``brain.save_image('out.png')`` may result in some or all of the labels +being absent from the saved ``out.png``. While we have implemented some +workarounds to help prevent this occurrance, we cannot guarantee it will +work. Thus we recommend that for critical non-interactive plotting (e.g., +if scripting figure generation for a paper) only a single view is used +with ``hemi`` set to ``'lh'``, ``'rh'``, or ``'both'``. This will use a single, +pure Mayavi window, thereby bypassing TraisUI entirely -- this helps +guarantee that drawing commands result in updated visual display. diff --git a/doc/index.rst b/doc/index.rst index 24b479b..d81c7ab 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -10,14 +10,14 @@ MRI and MEG data. PySurfer offers both a command-line interface designed to broadly replicate Freesurfer's Tksurfer program as well as a Python library -for writing scripts to efficiently explore complex datasets. +for writing scripts to efficiently explore complex datasets. Contents -------- .. toctree:: :maxdepth: 2 - + install examples/index.rst documentation/index.rst @@ -33,6 +33,8 @@ Scott Burns, Harvard Med. School MGH Martinos Center Martin Luessi, Harvard Med. School MGH Martinos Center +Eric Larson, University of Washington ILABS + License ------- diff --git a/doc/sphinxext/gen_rst.py b/doc/sphinxext/gen_rst.py index 60a887f..743f102 100755 --- a/doc/sphinxext/gen_rst.py +++ b/doc/sphinxext/gen_rst.py @@ -148,9 +148,11 @@ def generate_file_rst(fname, target_dir, src_dir, plot_gallery): this_template = rst_template last_dir = os.path.split(src_dir)[-1] # to avoid leading . in file names - if last_dir == '.': last_dir = '' - else: last_dir += '_' - short_fname = last_dir + fname + if last_dir == '.': + last_dir = '' + else: + last_dir += '_' + short_fname = last_dir + fname src_file = os.path.join(src_dir, fname) example_file = os.path.join(target_dir, fname) shutil.copyfile(src_file, example_file) @@ -158,56 +160,39 @@ def generate_file_rst(fname, target_dir, src_dir, plot_gallery): # generate the plot as png image if file name # starts with plot and if it is more recent than an # existing image. - if not os.path.exists( - os.path.join(target_dir, 'images')): + if not os.path.exists(os.path.join(target_dir, 'images')): os.makedirs(os.path.join(target_dir, 'images')) image_file = os.path.join(target_dir, 'images', image_name) if (not os.path.exists(image_file) or - os.stat(image_file).st_mtime <= - os.stat(src_file).st_mtime): + os.stat(image_file).st_mtime <= os.stat(src_file).st_mtime): print 'plotting %s' % fname import matplotlib.pyplot as plt plt.close('all') try: - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - mlab.close(all=True) - except: - pass - - try: - execfile(example_file, {'pl' : plt}) - facecolor = plt.gcf().get_facecolor() # hack to keep black bg + brain = None + global plt + global brain + execfile(example_file, globals()) + facecolor = plt.gcf().get_facecolor() # hack to keep black bg if facecolor == (0.0, 0.0, 0.0, 1.0): plt.savefig(image_file, facecolor='black') else: plt.savefig(image_file) - try: - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - e = mlab.get_engine() - if len(e.scenes) > 0: - mlab.savefig(image_file, size=(400, 400)) - except: - pass + brain.save_image(image_file) + brain.close() except: - print 80*'_' + print 80 * '_' print '%s is not compiling:' % fname traceback.print_exc() - print 80*'_' + print 80 * '_' this_template = plot_rst_template docstring, short_desc, end_row = extract_docstring(example_file) - f = open(os.path.join(target_dir, fname[:-2] + 'rst'),'w') - f.write( this_template % locals()) + f = open(os.path.join(target_dir, fname[:-2] + 'rst'), 'w') + f.write(this_template % locals()) f.flush() diff --git a/doc/surfer.cfg b/doc/surfer.cfg index 0139887..f381397 100644 --- a/doc/surfer.cfg +++ b/doc/surfer.cfg @@ -1,4 +1,6 @@ [visual] +width=600 +height=600 background=black cortex=classic default_view=lateral diff --git a/examples/example_data/index.rst b/examples/example_data/index.rst deleted file mode 100644 index e69de29..0000000 diff --git a/examples/example_data/meg_source_estimate-rh.stc b/examples/example_data/meg_source_estimate-rh.stc new file mode 100644 index 0000000..6becbf7 Binary files /dev/null and b/examples/example_data/meg_source_estimate-rh.stc differ diff --git a/examples/plot_from_volume.py b/examples/plot_from_volume.py index b55aacf..f5c4b3c 100644 --- a/examples/plot_from_volume.py +++ b/examples/plot_from_volume.py @@ -12,31 +12,37 @@ from surfer import Brain, io """Bring up the visualization""" -brain = Brain("fsaverage", "lh", "inflated", +brain = Brain("fsaverage", "split", "inflated", views=['lat', 'med'], config_opts=dict(background="white")) """Project the volume file and return as an array""" mri_file = "example_data/resting_corr.nii.gz" reg_file = "example_data/register.dat" -surf_data = io.project_volume_data(mri_file, "lh", reg_file) +surf_data_lh = io.project_volume_data(mri_file, "lh", reg_file) +surf_data_rh = io.project_volume_data(mri_file, "rh", reg_file) """ You can pass this array to the add_overlay method for a typical activation overlay (with thresholding, etc.) """ -brain.add_overlay(surf_data, min=.3, max=.7, name="ang_corr") +brain.add_overlay(surf_data_lh, min=.3, max=.7, name="ang_corr_lh", hemi='lh') +brain.add_overlay(surf_data_rh, min=.3, max=.7, name="ang_corr_rh", hemi='rh') """ You can also pass it to add_data for more control over the visualzation. Here we'll plot the whole range of correlations """ -brain.overlays["ang_corr"].remove() -brain.add_data(surf_data, -.7, .7, colormap="jet", alpha=.7) +for overlay in brain.overlays_dict["ang_corr_lh"]: + overlay.remove() +for overlay in brain.overlays_dict["ang_corr_rh"]: + overlay.remove() +brain.add_data(surf_data_lh, -.7, .7, colormap="jet", alpha=.7, hemi='lh') +brain.add_data(surf_data_rh, -.7, .7, colormap="jet", alpha=.7, hemi='rh') """ This overlay represents resting-state correlations with a seed in left angular gyrus. Let's plot that seed. """ seed_coords = (-45, -67, 36) -brain.add_foci(seed_coords, map_surface="white") +brain.add_foci(seed_coords, map_surface="white", hemi='lh') diff --git a/examples/plot_meg_inverse_solution.py b/examples/plot_meg_inverse_solution.py index 7973b71..4d8e531 100644 --- a/examples/plot_meg_inverse_solution.py +++ b/examples/plot_meg_inverse_solution.py @@ -14,46 +14,50 @@ from surfer.io import read_stc """ -define subject, surface and hemisphere +define subject, surface and hemisphere(s) to plot """ -subject_id, surface, hemi = 'fsaverage', 'inflated', 'lh' +subject_id, surface = 'fsaverage', 'inflated' +hemi = 'split' """ create Brain object for visualization """ -brain = Brain(subject_id, hemi, surface) +brain = Brain(subject_id, hemi, surface, + config_opts=dict(width=800, height=400)) """ read MNE dSPM inverse solution """ -stc_fname = os.path.join('example_data', - 'meg_source_estimate-' + hemi + '.stc') -stc = read_stc(stc_fname) +for hemi in ['lh', 'rh']: + stc_fname = os.path.join('example_data', + 'meg_source_estimate-' + hemi + '.stc') + stc = read_stc(stc_fname) -""" -data and vertices for which the data is defined -""" -data = stc['data'] -vertices = stc['vertices'] + """ + data and vertices for which the data is defined + """ + data = stc['data'] + vertices = stc['vertices'] -""" -time points in milliseconds -""" -time = 1e3 * np.linspace(stc['tmin'], - stc['tmin'] + data.shape[1] * stc['tstep'], - data.shape[1]) -""" -colormap to use -""" -colormap = 'hot' + """ + time points in milliseconds + """ + time = 1e3 * np.linspace(stc['tmin'], + stc['tmin'] + data.shape[1] * stc['tstep'], + data.shape[1]) + """ + colormap to use + """ + colormap = 'hot' -""" -label for time annotation -""" -time_label = 'time=%0.2f ms' + """ + label for time annotation + """ + time_label = 'time=%0.2f ms' -brain.add_data(data, colormap=colormap, vertices=vertices, smoothing_steps=10, - time=time, time_label=time_label) + brain.add_data(data, colormap=colormap, vertices=vertices, + smoothing_steps=10, time=time, time_label=time_label, + hemi=hemi) """ scale colormap and set time (index) to display diff --git a/examples/plot_morphometry.py b/examples/plot_morphometry.py index 324ac81..010c758 100644 --- a/examples/plot_morphometry.py +++ b/examples/plot_morphometry.py @@ -11,7 +11,7 @@ from surfer import Brain -brain = Brain("fsaverage", "lh", "pial", +brain = Brain("fsaverage", "both", "pial", views="frontal", config_opts=dict(background="dimgray")) """ @@ -19,16 +19,18 @@ recon-all live in a predicatble location, all you need to call the add_morphometry method with is the name of the measure you want. -Here, we'll look at cortical curvatuve values. +Here, we'll look at cortical curvatuve values, +and plot them for both hemispheres. """ brain.add_morphometry("curv") """ Each of the possible values is displayed in an appropriate full-color map, but you can also -display in grayscale. +display in grayscale. Here we only plot the +left hemisphere. """ -brain.add_morphometry("sulc", grayscale=True) +brain.add_morphometry("sulc", hemi='lh', grayscale=True) """ The Brain object can only hold one morphometry diff --git a/examples/plot_parcellation.py b/examples/plot_parcellation.py index 63b57e3..f370876 100644 --- a/examples/plot_parcellation.py +++ b/examples/plot_parcellation.py @@ -13,13 +13,15 @@ subject_id = 'fsaverage' -hemi = 'lh' +hemi = 'both' surface = 'inflated' +view = 'frontal' + """ Bring up the visualization """ -brain = Brain(subject_id, hemi, surface, +brain = Brain(subject_id, hemi, surface, views=view, config_opts={"cortex": "bone", "background": "ivory"}) @@ -37,8 +39,11 @@ """ You may also provide a full path to an annotation file -at an arbitray location on the disc. +at an arbitray location on the disc. You can also +plot things separately for the left and right hemispheres. """ subjects_dir = os.environ["SUBJECTS_DIR"] annot_path = pjoin(subjects_dir, subject_id, "label", "lh.aparc.annot") -brain.add_annotation(annot_path) +brain.add_annotation(annot_path, hemi='lh', borders=False) +annot_path = pjoin(subjects_dir, subject_id, "label", "rh.aparc.a2009s.annot") +brain.add_annotation(annot_path, hemi='rh', remove_existing=False) diff --git a/examples/plot_show_view.py b/examples/plot_show_view.py index 2dd7422..8b417cf 100644 --- a/examples/plot_show_view.py +++ b/examples/plot_show_view.py @@ -11,7 +11,7 @@ from surfer import Brain sub = 'fsaverage' -hemi = 'lh' +hemi = 'both' surf = 'inflated' brain = Brain(sub, hemi, surf) diff --git a/examples/plot_to_montage.py b/examples/plot_to_montage.py index 7054b43..54a4b2a 100644 --- a/examples/plot_to_montage.py +++ b/examples/plot_to_montage.py @@ -7,27 +7,27 @@ """ print __doc__ - from surfer import Brain sub = 'fsaverage' hemi = 'lh' surf = 'inflated' +bgcolor = 'w' -brain = Brain(sub, hemi, surf) +brain = Brain(sub, hemi, surf, config_opts={'background': bgcolor}) ############################################################################### -# Save a set of images as a montage -brain.save_montage('/tmp/fsaverage_h_montage.png', - ['l', 'v', 'm'], orientation='v') +# Get a set of images as a montage, note the data could be saved if desired +image = brain.save_montage(None, ['l', 'v', 'm'], orientation='v') brain.close() ############################################################################### # View created image -import Image import pylab as pl -image = Image.open('/tmp/fsaverage_h_montage.png') -fig = pl.figure(figsize=(5, 3)) -pl.imshow(image, origin='lower') +fig = pl.figure(figsize=(5, 3), facecolor=bgcolor) +ax = pl.axes(frameon=False) +ax.imshow(image, origin='upper') pl.xticks(()) pl.yticks(()) +pl.draw() +pl.show() diff --git a/surfer/__init__.py b/surfer/__init__.py index ca489da..99a6997 100644 --- a/surfer/__init__.py +++ b/surfer/__init__.py @@ -1,4 +1,7 @@ from io import Surface from viz import Brain, TimeViewer +from utils import verbose, set_log_level, set_log_file __version__ = "0.4.dev" + +set_log_level() diff --git a/surfer/_commandline.py b/surfer/_commandline.py index 162b741..2252d0a 100644 --- a/surfer/_commandline.py +++ b/surfer/_commandline.py @@ -32,7 +32,8 @@ description=help_text) parser.add_argument("subject_id", help="subject id as in subjects dir") -parser.add_argument("hemi", metavar="hemi", choices=["lh", "rh"], +parser.add_argument("hemi", metavar="hemi", choices=["lh", "rh", + "both", "split"], help="hemisphere to load") parser.add_argument("surf", help="surface mesh (e.g. 'pial', 'inflated')") @@ -66,3 +67,5 @@ help="colormap for binary cortex curvature") parser.add_argument("-title", help="title to use for the figure") +parser.add_argument("-views", nargs="*", default=['lat'], + help="view list (space-separated) to use") diff --git a/surfer/config.py b/surfer/config.py index 3822ad3..30a3557 100644 --- a/surfer/config.py +++ b/surfer/config.py @@ -7,13 +7,17 @@ [visual] background = black foreground = white -size = 800 +height = 800 +width = 800 cortex = classic default_view = lateral [overlay] min_thresh = 2.0 max_thresh = robust_max + +[options] +logging_level = INFO """) config = ConfigParser.ConfigParser() diff --git a/surfer/config_options.txt b/surfer/config_options.txt deleted file mode 100644 index af29e30..0000000 --- a/surfer/config_options.txt +++ /dev/null @@ -1,34 +0,0 @@ -Placeholder for config options until online docs are up and running. -May 17, 2011 - -Visual: -background: - - black - - white - - slate - - midnight - - sand - -size: - - any float - - default: 800 - -cortex: - - classic - - high_contrast - - low_contrast - - bone - ? - -min_thresh: - - any float - - robust_min - - actual_min - -max_thresh: - - any float - - robust_max - - actual_max - -default_view: - - any view string diff --git a/surfer/io.py b/surfer/io.py index c7eed08..9dad36a 100644 --- a/surfer/io.py +++ b/surfer/io.py @@ -8,6 +8,9 @@ import nibabel as nib from nibabel.spatialimages import ImageFileError +import logging +logger = logging.getLogger('surfer') + def _get_subjects_dir(subjects_dir=None): """Get the subjects directory from parameter or environment variable @@ -29,8 +32,8 @@ def _get_subjects_dir(subjects_dir=None): if 'SUBJECTS_DIR' in os.environ: subjects_dir = os.environ['SUBJECTS_DIR'] else: - raise ValueError('The subjects directory has to be specified using ' - 'either the subjects_dir parameter or the ' + raise ValueError('The subjects directory has to be specified ' + 'using either the subjects_dir parameter or the ' 'SUBJECTS_DIR environment variable.') if not os.path.exists(subjects_dir): @@ -217,11 +220,9 @@ def project_volume_data(filepath, hemi, reg_file=None, subject_id=None, Path to label file to constrain projection; otherwise uses cortex target_subject : string Subject to warp data to in surface space after projection - verbose : boolean - If true, print command line - + verbose : bool + If True, print the command used """ - # Set the basic commands cmd_list = ["mri_vol2surf", "--mov", filepath, @@ -301,7 +302,8 @@ class Surface(object): instead of the value set using the SUBJECTS_DIR environment variable. """ - def __init__(self, subject_id, hemi, surf, subjects_dir=None): + def __init__(self, subject_id, hemi, surf, subjects_dir=None, + offset=None): """Surface Parameters @@ -312,10 +314,17 @@ def __init__(self, subject_id, hemi, surf, subjects_dir=None): Which hemisphere to load surf : string Name of the surface to load (eg. inflated, orig ...) + offset : float | None + If 0.0, the surface will be offset such that the medial + wall is aligned with the origin. If None, no offset will + be applied. If != 0.0, an additional offset will be used. """ + if hemi not in ['lh', 'rh']: + raise ValueError('hemi must be "lh" or "rh') self.subject_id = subject_id self.hemi = hemi self.surf = surf + self.offset = offset subjects_dir = _get_subjects_dir(subjects_dir) self.data_path = pjoin(subjects_dir, subject_id) @@ -324,6 +333,11 @@ def load_geometry(self): surf_path = pjoin(self.data_path, "surf", "%s.%s" % (self.hemi, self.surf)) self.coords, self.faces = nib.freesurfer.read_geometry(surf_path) + if self.offset is not None: + if self.hemi == 'lh': + self.coords[:, 0] -= (np.max(self.coords[:, 0]) + self.offset) + else: + self.coords[:, 0] -= (np.min(self.coords[:, 0]) + self.offset) def save_geometry(self): surf_path = pjoin(self.data_path, "surf", diff --git a/surfer/tests/test_viz.py b/surfer/tests/test_viz.py index 4218c7f..3171861 100644 --- a/surfer/tests/test_viz.py +++ b/surfer/tests/test_viz.py @@ -34,10 +34,25 @@ def has_freesurfer(): 'Requires FreeSurfer command line tools') +def test_image(): + """Test image saving + """ + mlab.options.backend = 'auto' + brain = Brain(*std_args, config_opts=small_brain) + tmp_name = mktemp() + '.png' + brain.save_image(tmp_name) + brain.save_imageset(tmp_name, ['med', 'lat'], 'jpg') + brain.save_montage(tmp_name, ['l', 'v', 'm'], orientation='v') + brain.screenshot() + brain.close() + + def test_brains(): """Test plotting of Brain with different arguments """ - mlab.options.backend = 'test' + # testing backend breaks when passing in a figure, so we use 'auto' here + # (shouldn't affect usability, but it makes testing more annoying) + mlab.options.backend = 'auto' surfs = ['inflated', 'sphere'] hemis = ['lh', 'rh'] curvs = [True, False] @@ -47,7 +62,9 @@ def test_brains(): subj_dirs = [None, subj_dir] for surf, hemi, curv, title, co, fig, sd \ in zip(surfs, hemis, curvs, titles, config_opts, figs, subj_dirs): - Brain(subject_id, hemi, surf, curv, title, co, fig, sd) + print 'hello' + brain = Brain(subject_id, hemi, surf, curv, title, co, fig, sd) + brain.close() assert_raises(ValueError, Brain, subject_id, 'lh', 'inflated', subjects_dir='') @@ -62,6 +79,7 @@ def test_annot(): brain = Brain(*std_args) for a, b, p in zip(annots, borders, alphas): brain.add_annotation(a, b, p) + brain.close() def test_contour(): @@ -75,6 +93,7 @@ def test_contour(): line_width=2) brain.contour['surface'].actor.property.line_width = 1 brain.contour['surface'].contour.number_of_contours = 10 + brain.close() @requires_fs @@ -87,6 +106,7 @@ def test_data(): reg_file = pjoin(data_dir, 'register.dat') surf_data = io.project_volume_data(mri_file, "lh", reg_file) brain.add_data(surf_data, -.7, .7, colormap="jet", alpha=.7) + brain.close() def test_foci(): @@ -106,6 +126,7 @@ def test_foci(): scale_factor = 0.7 brain.add_foci(coords, coords_as_verts=True, scale_factor=scale_factor, color="#A52A2A") + brain.close() def test_label(): @@ -127,6 +148,7 @@ def test_label(): brain.add_label("V1", color="steelblue", alpha=.6) brain.add_label("V2", color="#FF6347", alpha=.6) brain.add_label("entorhinal", color=(.2, 1, .5), alpha=.6) + brain.close() def test_meg_inverse(): @@ -148,6 +170,7 @@ def test_meg_inverse(): brain.set_data_time_index(2) brain.scale_data_colormap(fmin=13, fmid=18, fmax=22, transparent=True) # viewer = TimeViewer(brain) + brain.close() def test_morphometry(): @@ -158,6 +181,7 @@ def test_morphometry(): brain.add_morphometry("curv") brain.add_morphometry("sulc", grayscale=True) brain.add_morphometry("thickness") + brain.close() def test_overlay(): @@ -189,6 +213,7 @@ def test_overlay(): brain.add_overlay(conjunct, 4, 30, name="conjunct") brain.overlays["conjunct"].pos_bar.lut_mode = "Purples" brain.overlays["conjunct"].pos_bar.visible = False + brain.close() def test_probabilistic_labels(): @@ -213,6 +238,7 @@ def test_probabilistic_labels(): brain.data["colorbar"].number_of_colors = 10 brain.data["colorbar"].number_of_labels = 11 + brain.close() def test_text(): @@ -221,6 +247,7 @@ def test_text(): mlab.options.backend = 'test' brain = Brain(*std_args) brain.add_text(0.1, 0.1, 'Hello', 'blah') + brain.close() def test_animate(): @@ -237,18 +264,6 @@ def test_animate(): brain.close() -def test_image(): - """Test image saving - """ - mlab.options.backend = 'auto' - brain = Brain(*std_args, config_opts=small_brain) - tmp_name = mktemp() + '.png' - brain.save_image(tmp_name) - brain.save_imageset(tmp_name, ['med', 'lat'], 'jpg') - brain.save_montage(tmp_name, ['l', 'v', 'm'], orientation='v') - brain.screenshot() - - def test_views(): """Test showing different views """ @@ -264,3 +279,4 @@ def test_views(): brain.show_view('dor') brain.show_view({'distance': 432}) brain.show_view({'azimuth': 135, 'elevation': 79}, roll=107) + brain.close() diff --git a/surfer/utils.py b/surfer/utils.py index 87bc77b..958962e 100644 --- a/surfer/utils.py +++ b/surfer/utils.py @@ -1,10 +1,157 @@ -import os +import logging +import warnings +import sys +from os import path as op +import inspect +from functools import wraps import numpy as np from scipy import sparse from scipy.spatial.distance import cdist from .io import Surface +from .config import config +logger = logging.getLogger('surfer') + + +############################################################################### +# LOGGING (courtesy of mne-python) + +def set_log_level(verbose=None, return_old_level=False): + """Convenience function for setting the logging level + + Parameters + ---------- + verbose : bool, str, int, or None + The verbosity of messages to print. If a str, it can be either DEBUG, + INFO, WARNING, ERROR, or CRITICAL. Note that these are for + convenience and are equivalent to passing in logging.DEBUG, etc. + For bool, True is the same as 'INFO', False is the same as 'WARNING'. + If None, the environment variable MNE_LOG_LEVEL is read, and if + it doesn't exist, defaults to INFO. + return_old_level : bool + If True, return the old verbosity level. + """ + if verbose is None: + verbose = config.get('options', 'logging_level') + elif isinstance(verbose, bool): + if verbose is True: + verbose = 'INFO' + else: + verbose = 'WARNING' + if isinstance(verbose, basestring): + verbose = verbose.upper() + logging_types = dict(DEBUG=logging.DEBUG, INFO=logging.INFO, + WARNING=logging.WARNING, ERROR=logging.ERROR, + CRITICAL=logging.CRITICAL) + if not verbose in logging_types: + raise ValueError('verbose must be of a valid type') + verbose = logging_types[verbose] + logger = logging.getLogger('surfer') + old_verbose = logger.level + logger.setLevel(verbose) + return (old_verbose if return_old_level else None) + + +class WrapStdOut(object): + """Ridiculous class to work around how doctest captures stdout""" + def __getattr__(self, name): + # Even more ridiculous than this class, this must be sys.stdout (not + # just stdout) in order for this to work (tested on OSX and Linux) + return getattr(sys.stdout, name) + + +def set_log_file(fname=None, output_format='%(message)s', overwrite=None): + """Convenience function for setting the log to print to a file + + Parameters + ---------- + fname : str, or None + Filename of the log to print to. If None, stdout is used. + To suppress log outputs, use set_log_level('WARN'). + output_format : str + Format of the output messages. See the following for examples: + http://docs.python.org/dev/howto/logging.html + e.g., "%(asctime)s - %(levelname)s - %(message)s". + overwrite : bool, or None + Overwrite the log file (if it exists). Otherwise, statements + will be appended to the log (default). None is the same as False, + but additionally raises a warning to notify the user that log + entries will be appended. + """ + logger = logging.getLogger('surfer') + handlers = logger.handlers + for h in handlers: + if isinstance(h, logging.FileHandler): + h.close() + logger.removeHandler(h) + if fname is not None: + if op.isfile(fname) and overwrite is None: + warnings.warn('Log entries will be appended to the file. Use ' + 'overwrite=False to avoid this message in the ' + 'future.') + mode = 'w' if overwrite is True else 'a' + lh = logging.FileHandler(fname, mode=mode) + else: + """ we should just be able to do: + lh = logging.StreamHandler(sys.stdout) + but because doctests uses some magic on stdout, we have to do this: + """ + lh = logging.StreamHandler(WrapStdOut()) + + lh.setFormatter(logging.Formatter(output_format)) + # actually add the stream handler + logger.addHandler(lh) + + +def verbose(function): + """Decorator to allow functions to override default log level + + Do not call this function directly to set the global verbosity level, + instead use set_log_level(). + + Parameters (to decorated function) + ---------------------------------- + verbose : bool, str, int, or None + The level of messages to print. If a str, it can be either DEBUG, + INFO, WARNING, ERROR, or CRITICAL. Note that these are for + convenience and are equivalent to passing in logging.DEBUG, etc. + For bool, True is the same as 'INFO', False is the same as 'WARNING'. + None defaults to using the current log level [e.g., set using + mne.set_log_level()]. + """ + arg_names = inspect.getargspec(function).args + # this wrap allows decorated functions to be pickled (e.g., for parallel) + + @wraps(function) + def dec(*args, **kwargs): + # Check if the first arg is "self", if it has verbose, make it default + if len(arg_names) > 0 and arg_names[0] == 'self': + default_level = getattr(args[0], 'verbose', None) + else: + default_level = None + verbose_level = kwargs.get('verbose', default_level) + if verbose_level is not None: + old_level = set_log_level(verbose_level, True) + # set it back if we get an exception + try: + ret = function(*args, **kwargs) + except: + set_log_level(old_level) + raise + set_log_level(old_level) + return ret + else: + return function(*args, **kwargs) + + # set __wrapped__ attribute so ?? in IPython gets the right source + dec.__wrapped__ = function + + return dec + + +############################################################################### +# USEFUL FUNCTIONS def find_closest_vertices(surface_coords, point_coords): """Return the vertices on a surface mesh closest to some given coordinates. @@ -78,7 +225,8 @@ def mesh_edges(faces): return edges -def smoothing_matrix(vertices, adj_mat, smoothing_steps=20): +@verbose +def smoothing_matrix(vertices, adj_mat, smoothing_steps=20, verbose=None): """Create a smoothing matrix which can be used to interpolate data defined for a subset of vertices onto mesh with an adjancency matrix given by adj_mat. @@ -95,6 +243,8 @@ def smoothing_matrix(vertices, adj_mat, smoothing_steps=20): N x N adjacency matrix of the full mesh smoothing_steps : int or None number of smoothing steps (Default: 20) + verbose : bool, str, int, or None + If not None, override default verbose level (see surfer.verbose). Returns ------- @@ -103,7 +253,7 @@ def smoothing_matrix(vertices, adj_mat, smoothing_steps=20): """ from scipy import sparse - print "Updating smoothing matrix, be patient.." + logger.info("Updating smoothing matrix, be patient..") e = adj_mat.copy() e.data[e.data == 2] = 1 @@ -122,11 +272,11 @@ def smoothing_matrix(vertices, adj_mat, smoothing_steps=20): smooth_mat = scale_mat * e_use[idx_use, :] * smooth_mat - print "Smoothing matrix creation, step %d" % (k + 1) + logger.info("Smoothing matrix creation, step %d" % (k + 1)) if smoothing_steps is None and len(idx_use) >= n_vertices: break - # Make sure the smooting matrix has the right number of rows + # Make sure the smoothing matrix has the right number of rows # and is in COO format smooth_mat = smooth_mat.tocoo() smooth_mat = sparse.coo_matrix((smooth_mat.data, @@ -138,8 +288,9 @@ def smoothing_matrix(vertices, adj_mat, smoothing_steps=20): return smooth_mat +@verbose def coord_to_label(subject_id, coord, label, hemi='lh', n_steps=30, - map_surface='white', coord_as_vert=False): + map_surface='white', coord_as_vert=False, verbose=None): """Create label from MNI coordinate Parameters @@ -158,6 +309,8 @@ def coord_to_label(subject_id, coord, label, hemi='lh', n_steps=30, The surface name used to find the closest point coord_as_vert : bool whether the coords parameter should be interpreted as vertex ids + verbose : bool, str, int, or None + If not None, override default verbose level (see surfer.verbose). """ geo = Surface(subject_id, hemi, map_surface) geo.load_geometry() @@ -176,7 +329,7 @@ def coord_to_label(subject_id, coord, label, hemi='lh', n_steps=30, idx = np.where(data.ravel() > 0)[0] # Write label label_fname = label + '-' + hemi + '.label' - print "Saving label : %s" % label_fname + logger.info("Saving label : %s" % label_fname) f = open(label_fname, 'w') f.write('#label at %s from subject %s\n' % (coord, subject_id)) f.write('%d\n' % len(idx)) diff --git a/surfer/viz.py b/surfer/viz.py index a296a17..1940691 100644 --- a/surfer/viz.py +++ b/surfer/viz.py @@ -3,8 +3,7 @@ from warnings import warn import numpy as np -from scipy import stats -from scipy import ndimage +from scipy import stats, ndimage, misc from matplotlib.colors import colorConverter import nibabel as nib @@ -13,31 +12,35 @@ from . import utils from .io import Surface, _get_subjects_dir from .config import config +from .utils import verbose + +import logging +logger = logging.getLogger('surfer') try: from traits.api import (HasTraits, Range, Int, Float, - Bool, Enum, on_trait_change) + Bool, Enum, on_trait_change, Instance) except ImportError: - from enthought.traits.api import (HasTraits, Range, Int, Float, \ - Bool, Enum, on_trait_change) - + from enthought.traits.api import (HasTraits, Range, Int, Float, + Bool, Enum, on_trait_change, Instance) lh_viewdict = {'lateral': {'v': (180., 90.), 'r': 90.}, - 'medial': {'v': (0., 90.), 'r': -90.}, - 'rostral': {'v': (90., 90.), 'r': -180.}, - 'caudal': {'v': (270., 90.), 'r': 0.}, - 'dorsal': {'v': (180., 0.), 'r': 90.}, - 'ventral': {'v': (180., 180.), 'r': 90.}, - 'frontal': {'v': (120., 80.), 'r': 106.739}, - 'parietal': {'v': (-120., 60.), 'r': 49.106}} + 'medial': {'v': (0., 90.), 'r': -90.}, + 'rostral': {'v': (90., 90.), 'r': -180.}, + 'caudal': {'v': (270., 90.), 'r': 0.}, + 'dorsal': {'v': (180., 0.), 'r': 90.}, + 'ventral': {'v': (180., 180.), 'r': 90.}, + 'frontal': {'v': (120., 80.), 'r': 106.739}, + 'parietal': {'v': (-120., 60.), 'r': 49.106}} rh_viewdict = {'lateral': {'v': (180., -90.), 'r': -90.}, - 'medial': {'v': (0., -90.), 'r': 90.}, - 'rostral': {'v': (-90., -90.), 'r': 180.}, - 'caudal': {'v': (90., -90.), 'r': 0.}, - 'dorsal': {'v': (180., 0.), 'r': 90.}, - 'ventral': {'v': (180., 180.), 'r': 90.}, - 'frontal': {'v': (60., 80.), 'r': -106.739}, - 'parietal': {'v': (-60., 60.), 'r': -49.106}} + 'medial': {'v': (0., -90.), 'r': 90.}, + 'rostral': {'v': (-90., -90.), 'r': 180.}, + 'caudal': {'v': (90., -90.), 'r': 0.}, + 'dorsal': {'v': (180., 0.), 'r': 90.}, + 'ventral': {'v': (180., 180.), 'r': 90.}, + 'frontal': {'v': (60., 80.), 'r': -106.739}, + 'parietal': {'v': (-60., 60.), 'r': -49.106}} +viewdicts = dict(lh=lh_viewdict, rh=rh_viewdict) def make_montage(filename, fnames, orientation='h', colorbar=None, @@ -47,9 +50,11 @@ def make_montage(filename, fnames, orientation='h', colorbar=None, Parameters ---------- filename : str - The name of the file, e.g, 'montage.png' - fnames : list of str - The images to make the montage off. + The name of the file, e.g, 'montage.png'. If None, the image + will not be saved. + fnames : list of str | list of array + The images to make the montage of. Can be a list of filenames + or a list of image data arrays. orientation : 'h' | 'v' The orientation of the montage: horizontal or vertical colorbar : None | list of int @@ -57,9 +62,21 @@ def make_montage(filename, fnames, orientation='h', colorbar=None, is present. border_size : int The size of the border to keep. + + Returns + ------- + out : array + The montage image data array. """ import Image - images = map(Image.open, fnames) + # This line is only necessary to overcome a PIL bug, see: + # http://stackoverflow.com/questions/10854903/what-is-causing- + # dimension-dependent-attributeerror-in-pil-fromarray-function + fnames = [f if isinstance(f, basestring) else f.copy() for f in fnames] + if isinstance(fnames[0], basestring): + images = map(Image.open, fnames) + else: + images = map(Image.fromarray, fnames) # get bounding box for cropping boxes = [] for ix, im in enumerate(images): @@ -118,184 +135,539 @@ def make_montage(filename, fnames, orientation='h', colorbar=None, pos = (0, x) x += i.size[1] new.paste(i, pos) - try: - new.save(filename) - except Exception: - print("Error saving %s" % filename) + if filename is not None: + try: + new.save(filename) + except Exception: + print("Error saving %s" % filename) + return np.array(new) -class Brain(object): - """Brain object for visualizing with mlab.""" +def _prepare_data(data): + """Ensure data is float64 and has proper endianness. - def __init__(self, subject_id, hemi, surf, - curv=True, title=None, config_opts={}, - figure=None, subjects_dir=None): - """Initialize a Brain object with Freesurfer-specific data. + Note: this is largely aimed at working around a Mayavi bug. - Parameters - ---------- - subject_id : str - subject name in Freesurfer subjects dir - hemi : str - hemisphere id (ie 'lh' or 'rh') - surf : geometry name - freesurfer surface mesh name (ie 'white', 'inflated', etc.) - curv : boolean - if true, loads curv file and displays binary curvature - (default: True) - title : str - title for the mayavi figure - config_opts : dict - options to override visual options in config file - figure : instance of mayavi.core.scene.Scene | None - If None, the last figure will be cleaned and a new figure will - be created. - subjects_dir : str | None - If not None, this directory will be used as the subjects directory - instead of the value set using the SUBJECTS_DIR environment - variable. - """ + """ + data = data.copy() + data = data.astype(np.float64) + if data.dtype.byteorder == '>': + data.byteswap(True) + return data + + +def _force_render(figures, backend): + """Ensure plots are updated before properties are used""" + try: + from mayavi import mlab + assert mlab + except: + from enthought.mayavi import mlab + if not isinstance(figures, list): + figures = [[figures]] + for ff in figures: + for f in ff: + f.render() + mlab.draw(figure=f) + if backend == 'TraitsUI': + from pyface.api import GUI + _gui = GUI() + orig_val = _gui.busy + _gui.set_busy(busy=True) + _gui.process_events() + _gui.set_busy(busy=orig_val) + _gui.process_events() + + +def _make_viewer(figure, n_row, n_col, title, scene_size, offscreen): + """Triage viewer creation + + If n_row == n_col == 1, then we can use a Mayavi figure, which + generally guarantees that things will be drawn before control + is returned to the command line. With the multi-view, TraitsUI + unfortunately has no such support, so we only use it if needed. + """ + try: + from mayavi import mlab + assert mlab + except: + from enthought.mayavi import mlab + if figure is None: + # spawn scenes + h, w = scene_size + if offscreen is True: + orig_val = mlab.options.offscreen + mlab.options.offscreen = True + figures = [[mlab.figure(size=(h / n_row, w / n_col)) + for _ in xrange(n_col)] for __ in xrange(n_row)] + mlab.options.offscreen = orig_val + _v = None + else: + # Triage: don't make TraitsUI if we don't have to + if n_row == 1 and n_col == 1: + figure = mlab.figure(title, size=(w, h)) + mlab.clf(figure) + figures = [[figure]] + _v = None + else: + window = _MlabGenerator(n_row, n_col, w, h, title) + figures, _v = window._get_figs_view() + else: + if not isinstance(figure, (list, tuple)): + figure = [figure] + if not len(figure) == n_row * n_col: + raise ValueError('For the requested view, figure must be a ' + 'list or tuple with exactly %i elements, ' + 'not %i' % (n_row * n_col, len(figure))) + _v = None + figures = [figure[slice(ri * n_col, (ri + 1) * n_col)] + for ri in range(n_row)] + return figures, _v + + +class _MlabGenerator(HasTraits): + """TraitsUI mlab figure generator""" + try: + from traitsui.api import View + except ImportError: try: - from mayavi import mlab + from traits.ui.api import View except ImportError: - from enthought.mayavi import mlab + from enthought.traits.ui.api import View - # Set the identifying info - self.subject_id = subject_id - self.hemi = hemi - if self.hemi == 'lh': - self.viewdict = lh_viewdict - elif self.hemi == 'rh': - self.viewdict = rh_viewdict - self.surf = surf - - # Initialize the mlab figure - if title is None: - title = subject_id - self._set_scene_properties(config_opts) - if figure is None: - figure = mlab.figure(title, **self.scene_properties) - mlab.clf() - self._f = figure + view = Instance(View) - # Get the subjects directory from parameter or env. var - self.subjects_dir = _get_subjects_dir(subjects_dir=subjects_dir) + def __init__(self, n_row, n_col, width, height, title, **traits): + HasTraits.__init__(self, **traits) + try: + from mayavi.tools.mlab_scene_model import MlabSceneModel + assert MlabSceneModel + except: + from enthought.mayavi.tools.mlab_scene_model import MlabSceneModel + + self.mlab_names = [] + self.n_row = n_row + self.n_col = n_col + self.width = width + self.height = height + for fi in xrange(n_row * n_col): + name = 'mlab_view%03g' % fi + self.mlab_names.append(name) + self.add_trait(name, Instance(MlabSceneModel, ())) + self.view = self._get_gen_view() + self._v = self.edit_traits(view=self.view) + self._v.title = title + + def _get_figs_view(self): + figures = [] + ind = 0 + for ri in range(self.n_row): + rfigs = [] + for ci in range(self.n_col): + x = getattr(self, self.mlab_names[ind]) + rfigs.append(x.mayavi_scene) + ind += 1 + figures.append(rfigs) + return figures, self._v + + def _get_gen_view(self): + try: + from mayavi.core.ui.api import SceneEditor + from mayavi.core.ui.mayavi_scene import MayaviScene + assert SceneEditor + assert MayaviScene + except: + from enthought.mayavi.core.ui.api import SceneEditor + from enthought.mayavi.core.ui.mayavi_scene import MayaviScene + try: + from traitsui.api import (View, Item, VGroup, HGroup) + except ImportError: + try: + from traits.ui.api import (View, Item, VGroup, HGroup) + except ImportError: + from enthought.traits.ui.api import (View, Item, + VGroup, HGroup) + ind = 0 + va = [] + for ri in xrange(self.n_row): + ha = [] + for ci in xrange(self.n_col): + ha += [Item(name=self.mlab_names[ind], style='custom', + resizable=True, show_label=False, + editor=SceneEditor(scene_class=MayaviScene))] + ind += 1 + va += [HGroup(*ha)] + view = View(VGroup(*va), resizable=True, + height=self.height, width=self.width) + return view - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - # Set the lights so they are oriented by hemisphere - self._orient_lights() - # Initialize a Surface object as the geometry - self._geo = Surface(subject_id, hemi, surf, - subjects_dir=self.subjects_dir) +class Brain(object): + """Class for visualizing a brain using multiple views in mlab - # Load in the geometry and (maybe) curvature - self._geo.load_geometry() - if curv: - self._geo.load_curvature() - curv_data = self._geo.bin_curv - meshargs = dict(scalars=curv_data) - else: - curv_data = None - meshargs = dict() + Parameters + ---------- + subject_id : str + subject name in Freesurfer subjects dir + hemi : str + hemisphere id (ie 'lh', 'rh', 'both', or 'split'). In the case + of 'both', both hemispheres are shown in the same window. + In the case of 'split' hemispheres are displayed side-by-side + in different viewing panes. + surf : geometry name + freesurfer surface mesh name (ie 'white', 'inflated', etc.) + curv : boolean + if true, loads curv file and displays binary curvature + (default: True) + title : str + title for the window + config_opts : dict + options to override visual options in config file + figure : list of instances of mayavi.core.scene.Scene | None + If None, a new window will be created with the appropriate + views. + subjects_dir : str | None + If not None, this directory will be used as the subjects directory + instead of the value set using the SUBJECTS_DIR environment + variable. + views : list | str + views to use + show_toolbar : bool + If True, toolbars will be shown for each view. + offscreen : bool + If True, rendering will be done offscreen (not shown). Useful + mostly for generating images or screenshots, but can be buggy. + Use at your own risk. + + Attributes + ---------- + brains : list + List of the underlying brain instances. + """ + def __init__(self, subject_id, hemi, surf, curv=True, title=None, + config_opts={}, figure=None, subjects_dir=None, + views=['lat'], show_toolbar=False, offscreen=False): + col_dict = dict(lh=1, rh=1, both=1, split=2) + n_col = col_dict[hemi] + if not hemi in col_dict.keys(): + raise ValueError('hemi must be one of [%s], not %s' + % (', '.join(col_dict.keys()), hemi)) + # Get the subjects directory from parameter or env. var + subjects_dir = _get_subjects_dir(subjects_dir=subjects_dir) - # mlab pipeline mesh for geomtery - self._geo_mesh = mlab.pipeline.triangular_mesh_source( - self._geo.x, self._geo.y, self._geo.z, - self._geo.faces, **meshargs) + self._hemi = hemi + if title is None: + title = subject_id + self.subject_id = subject_id - # mlab surface for the geometry - if curv: - colormap, vmin, vmax, reverse = self._get_geo_colors(config_opts) - self._geo_surf = mlab.pipeline.surface(self._geo_mesh, - colormap=colormap, vmin=vmin, vmax=vmax) - if reverse: - curv_bar = mlab.scalarbar(self._geo_surf) - curv_bar.reverse_lut = True - curv_bar.visible = False + if not isinstance(views, list): + views = [views] + n_row = len(views) + + # load geometry for one or both hemispheres as necessary + offset = None if hemi != 'both' else 0.0 + self.geo = dict() + if hemi in ['split', 'both']: + geo_hemis = ['lh', 'rh'] + elif hemi == 'lh': + geo_hemis = ['lh'] + elif hemi == 'rh': + geo_hemis = ['rh'] else: - self._geo_surf = mlab.pipeline.surface(self._geo_mesh, - color=(.5, .5, .5)) - + raise ValueError('bad hemi value') + for h in geo_hemis: + # Initialize a Surface object as the geometry + geo = Surface(subject_id, h, surf, subjects_dir, offset) + # Load in the geometry and (maybe) curvature + geo.load_geometry() + if curv: + geo.load_curvature() + self.geo[h] = geo + + # deal with making figures + self._set_window_properties(config_opts) + figures, _v = _make_viewer(figure, n_row, n_col, title, + self._scene_size, offscreen) + self._figures = figures + self._v = _v + self._window_backend = 'Mayavi' if self._v is None else 'TraitsUI' + for ff in self._figures: + for f in ff: + if f.scene is not None: + f.scene.background = self._bg_color + f.scene.foreground = self._fg_color + # force rendering so scene.lights exists + _force_render(self._figures, self._window_backend) + self.toggle_toolbars(show_toolbar) + _force_render(self._figures, self._window_backend) + self._toggle_render(False) + + # fill figures with brains + kwargs = dict(surf=surf, curv=curv, title=None, + config_opts=config_opts, subjects_dir=subjects_dir, + bg_color=self._bg_color, offset=offset) + brains = [] + brain_matrix = [] + for ri, view in enumerate(views): + brain_row = [] + for hi, h in enumerate(['lh', 'rh']): + if not (hemi in ['lh', 'rh'] and h != hemi): + ci = hi if hemi == 'split' else 0 + kwargs['hemi'] = h + kwargs['geo'] = self.geo[h] + kwargs['figure'] = figures[ri][ci] + kwargs['backend'] = self._window_backend + brain = _Hemisphere(subject_id, **kwargs) + brain.show_view(view) + brains += [dict(row=ri, col=ci, brain=brain, hemi=h)] + brain_row += [brain] + brain_matrix += [brain_row] + self._toggle_render(True) + self._original_views = views + self._brain_list = brains + for brain in self._brain_list: + brain['brain']._orient_lights() + self.brains = [b['brain'] for b in brains] + self.brain_matrix = np.array(brain_matrix) + self.subjects_dir = subjects_dir # Initialize the overlay and label dictionaries - self.overlays = dict() - self.labels = dict() - self.foci = dict() - self.texts = dict() + self.foci_dict = dict() + self.labels_dict = dict() + self.overlays_dict = dict() + self.contour_list = [] + self.morphometry_list = [] + self.annot_list = [] + self.data_dict = dict(lh=None, rh=None) + # note that texts gets treated differently + self.texts_dict = dict() + self.n_times = None + + ########################################################################### + # HELPERS + def _toggle_render(self, state, views=None): + """Turn rendering on (True) or off (False)""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab - # Bring up the lateral view - self.show_view(config.get("visual", "default_view")) + figs = [] + [figs.extend(f) for f in self._figures] + if views is None: + views = [None] * len(figs) + for vi, (_f, view) in enumerate(zip(figs, views)): + if state is False and view is None: + views[vi] = mlab.view(figure=_f) + + # Testing backend doesn't have this option + if mlab.options.backend != 'test': + _f.scene.disable_render = not state + + if state is True and view is not None: + mlab.draw(figure=_f) + mlab.view(*view, figure=_f) + # let's do the ugly force draw + if state is True: + _force_render(self._figures, self._window_backend) + return views + + def _set_window_properties(self, config_opts): + """Set window properties using config_opts""" + # old option "size" sets both width and height + if "size" in config_opts: + width = config_opts["size"] + height = width + + if "width" in config_opts: + width = config_opts["width"] + else: + width = config.getfloat("visual", "width") + if 'height' in config_opts: + height = config_opts['height'] + else: + height = config.getfloat("visual", "height") + self._scene_size = (height, width) - # Turn disable render off so that it displays (option doesn't exist - # in testing mode) - if mlab.options.backend != 'test': - self._f.scene.disable_render = False + try: + bg_color_name = config_opts['background'] + except KeyError: + bg_color_name = config.get("visual", "background") + if bg_color_name is not None: + bg_color_code = colorConverter.to_rgb(bg_color_name) + else: + bg_color_code = None + self._bg_color = bg_color_code - def show_view(self, view=None, roll=None): - """Orient camera to display view + try: + fg_color_name = config_opts['foreground'] + except KeyError: + fg_color_name = config.get("visual", "foreground") + fg_color_code = colorConverter.to_rgb(fg_color_name) + self._fg_color = fg_color_code - Parameters - ---------- - view : {'lateral' | 'medial' | 'rostral' | 'caudal' | - 'dorsal' | 'ventral' | 'frontal' | 'parietal' | - dict} - brain surface to view or kwargs to pass to mlab.view() + def get_data_properties(self): + """ Get properties of the data shown Returns ------- - cv: tuple - tuple returned from mlab.view - cr: float - current camera roll + props : dict + Dictionary with data properties + props["fmin"] : minimum colormap + props["fmid"] : midpoint colormap + props["fmax"] : maximum colormap + props["transparent"] : lower part of colormap transparent? + props["time"] : time points + props["time_idx"] : current time index + props["smoothing_steps"] : number of smoothing steps """ - if isinstance(view, basestring): - try: - vd = self.xfm_view(view, 'd') - view = dict(azimuth=vd['v'][0], elevation=vd['v'][1]) - roll = vd['r'] - except ValueError as v: - print(v) - raise - cv, cr = self.__view(view, roll) - return (cv, cr) + props = dict() + keys = ['fmin', 'fmid', 'fmax', 'transparent', 'time', 'time_idx', + 'smoothing_steps'] + try: + if self.data_dict['lh'] is not None: + hemi = 'lh' + else: + hemi = 'rh' + for key in keys: + props[key] = self.data_dict[hemi][key] + except KeyError: + # The user has not added any data + for key in keys: + props[key] = 0 + return props - def __view(self, viewargs=None, roll=None): - """Wrapper for mlab.view() + def toggle_toolbars(self, show=None): + """Toggle toolbar display Parameters ---------- - viewargs: dict - mapping with keys corresponding to mlab.view args - roll: num - int or float to set camera roll - - Returns - ------- - camera settings: tuple - view settings, roll setting - + show : bool | None + If None, the state is toggled. If True, the toolbar will + be shown, if False, hidden. """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - if viewargs: - viewargs['reset_roll'] = True - viewargs['figure'] = self._f - mlab.view(**viewargs) - if not roll is None: - mlab.roll(roll, figure=self._f) + # don't do anything if testing is on + if self._figures[0][0].scene is not None: + # this may not work if QT is not the backend (?), or in testing + if hasattr(self._figures[0][0].scene, 'scene_editor'): + # Within TraitsUI + bars = [f.scene.scene_editor._tool_bar + for ff in self._figures for f in ff] + else: + # Mayavi figure + bars = [f.scene._tool_bar for ff in self._figures for f in ff] - view = mlab.view(figure=self._f) - roll = mlab.roll(figure=self._f) + if show is None: + if hasattr(bars[0], 'isVisible'): + # QT4 + show = not bars[0].isVisible() + else: + # WX + show = not bars[0].Shown() + for bar in bars: + if hasattr(bar, 'setVisible'): + bar.setVisible(show) + else: + bar.Show(show) + + def _get_one_brain(self, d, name): + """Helper for various properties""" + if len(self.brains) > 1: + raise ValueError('Cannot access brain.%s when more than ' + 'one view is plotted. Use brain.brain_matrix ' + 'or brain.brains.' % name) + if isinstance(d, dict): + out = dict() + for key, value in d.iteritems(): + out[key] = value[0] + else: + out = d[0] + return out + + @property + def overlays(self): + """Wrap to overlays""" + return self._get_one_brain(self.overlays_dict, 'overlays') + + @property + def foci(self): + """Wrap to foci""" + return self._get_one_brain(self.foci_dict, 'foci') + + @property + def labels(self): + """Wrap to labels""" + return self._get_one_brain(self.labels_dict, 'labels') + + @property + def contour(self): + """Wrap to contour""" + return self._get_one_brain(self.contour_list, 'contour') + + @property + def annot(self): + """Wrap to annot""" + return self._get_one_brain(self.annot_list, 'contour') + + @property + def texts(self): + """Wrap to texts""" + self._get_one_brain([[]], 'texts') + out = dict() + for key, val in self.texts_dict.iteritems(): + out[key] = val['text'] + return out + + @property + def _geo(self): + """Wrap to _geo""" + self._get_one_brain([[]], '_geo') + if self.geo['lh'] is not None: + return self.geo['lh'] + else: + return self.geo['rh'] + + @property + def data(self): + """Wrap to data""" + self._get_one_brain([[]], 'data') + if self.data_dict['lh'] is not None: + data = self.data_dict['lh'].copy() + else: + data = self.data_dict['rh'].copy() + if 'colorbars' in data: + data['colorbar'] = data['colorbars'][0] + return data - return view, roll + def _check_hemi(self, hemi): + """Check for safe single-hemi input, returns str""" + if hemi is None: + if self._hemi not in ['lh', 'rh']: + raise ValueError('hemi must not be None when both ' + 'hemispheres are displayed') + else: + hemi = self._hemi + elif hemi not in ['lh', 'rh']: + extra = ' or None' if self._hemi in ['lh', 'rh'] else '' + raise ValueError('hemi must be either "lh" or "rh"' + extra) + return hemi + + def _check_hemis(self, hemi): + """Check for safe dual or single-hemi input, returns list""" + if hemi is None: + if self._hemi not in ['lh', 'rh']: + hemi = ['lh', 'rh'] + else: + hemi = [self._hemi] + elif hemi not in ['lh', 'rh']: + extra = ' or None' if self._hemi in ['lh', 'rh'] else '' + raise ValueError('hemi must be either "lh" or "rh"' + extra) + else: + hemi = [hemi] + return hemi - def _read_scalar_data(self, source, name=None, cast=True): + def _read_scalar_data(self, source, hemi, name=None, cast=True): """Load in scalar data from an image stored in a file or an array Parameters @@ -317,14 +689,13 @@ def _read_scalar_data(self, source, name=None, cast=True): if no name was provided, deduces the name if filename was given as a source """ - # If source is a string, try to load a file if isinstance(source, basestring): if name is None: basename = os.path.basename(source) if basename.endswith(".gz"): basename = basename[:-3] - if basename.startswith("%s." % self.hemi): + if basename.startswith("%s." % hemi): basename = basename[3:] name = os.path.splitext(basename)[0] scalar_data = io.read_scalar_data(source) @@ -333,13 +704,65 @@ def _read_scalar_data(self, source, name=None, cast=True): scalar_data = source if cast: - if (scalar_data.dtype.char == 'f' and - scalar_data.dtype.itemsize < 8): + if (scalar_data.dtype.char == 'f' + and scalar_data.dtype.itemsize < 8): scalar_data = scalar_data.astype(np.float) return scalar_data, name - def add_overlay(self, source, min=None, max=None, sign="abs", name=None): + def _get_display_range(self, scalar_data, min, max, sign): + if scalar_data.min() >= 0: + sign = "pos" + elif scalar_data.max() <= 0: + sign = "neg" + + # Get data with a range that will make sense for automatic thresholding + if sign == "neg": + range_data = np.abs(scalar_data[np.where(scalar_data < 0)]) + elif sign == "pos": + range_data = scalar_data[np.where(scalar_data > 0)] + else: + range_data = np.abs(scalar_data) + + # Get the min and max from among various places + if min is None: + try: + min = config.getfloat("overlay", "min_thresh") + except ValueError: + min_str = config.get("overlay", "min_thresh") + if min_str == "robust_min": + min = stats.scoreatpercentile(range_data, 2) + elif min_str == "actual_min": + min = range_data.min() + else: + min = 2.0 + warn("The 'min_thresh' value in your config value must be " + "a float, 'robust_min', or 'actual_min', but it is " + "%s. I'm setting the overlay min to the config " + "default of 2" % min) + + if max is None: + try: + max = config.getfloat("overlay", "max_thresh") + except ValueError: + max_str = config.get("overlay", "max_thresh") + if max_str == "robust_max": + max = stats.scoreatpercentile(scalar_data, 98) + elif max_str == "actual_max": + max = range_data.max() + else: + max = stats.scoreatpercentile(range_data, 98) + warn("The 'max_thresh' value in your config value must be " + "a float, 'robust_min', or 'actual_min', but it is " + "%s. I'm setting the overlay min to the config " + "default of robust_max" % max) + + return min, max + + ########################################################################### + # ADDING DATA PLOTS + def add_overlay(self, source, min=None, max=None, sign="abs", name=None, + hemi=None): """Add an overlay to the overlay dict from a file or array. Parameters @@ -354,43 +777,33 @@ def add_overlay(self, source, min=None, max=None, sign="abs", name=None): whether positive, negative, or both values should be displayed name : str name for the overlay in the internal dictionary - + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, an error will + be thrown. """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - scalar_data, name = self._read_scalar_data(source, name) - + hemi = self._check_hemi(hemi) + # load data here + scalar_data, name = self._read_scalar_data(source, hemi, name=name) min, max = self._get_display_range(scalar_data, min, max, sign) - - if name in self.overlays: - "%s%d" % (name, len(self.overlays) + 1) - if not sign in ["abs", "pos", "neg"]: raise ValueError("Overlay sign must be 'abs', 'pos', or 'neg'") - - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - view = mlab.view() - self.overlays[name] = Overlay(scalar_data, self._geo, min, max, sign) - for bar in ["pos_bar", "neg_bar"]: - try: - self._format_cbar_text(getattr(self.overlays[name], bar)) - except AttributeError: - pass - - # Testing backend doesn't have disable_render, option view == None - if mlab.options.backend != 'test': - mlab.view(*view) - self._f.scene.disable_render = False + old = OverlayData(scalar_data, self.geo[hemi], min, max, sign) + ol = [] + views = self._toggle_render(False) + for brain in self._brain_list: + if brain['hemi'] == hemi: + ol.append(brain['brain'].add_overlay(old)) + if name in self.overlays_dict: + name = "%s%d" % (name, len(self.overlays_dict) + 1) + self.overlays_dict[name] = ol + self._toggle_render(True, views) def add_data(self, array, min=None, max=None, thresh=None, colormap="blue-red", alpha=1, vertices=None, smoothing_steps=20, time=None, - time_label="time index=%d", colorbar=True): + time_label="time index=%d", colorbar=True, + hemi=None): """Display data from a numpy array on the surface. This provides a similar interface to add_overlay, but it displays @@ -435,22 +848,12 @@ def add_data(self, array, min=None, max=None, thresh=None, format of the time label (or None for no label) colorbar : bool whether to add a colorbar to the figure + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, an error will + be thrown. """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - view = mlab.view() - - # Possibly remove old data - if hasattr(self, "data"): - self.data["surface"].remove() - if 'colorbar' in self.data: - self.data["colorbar"].remove() + hemi = self._check_hemi(hemi) if min is None: min = array.min() @@ -458,10 +861,10 @@ def add_data(self, array, min=None, max=None, thresh=None, max = array.max() # Create smoothing matrix if necessary - if len(array) < self._geo.x.shape[0]: - if vertices == None: + if len(array) < self.geo[hemi].x.shape[0]: + if vertices is None: raise ValueError("len(data) < nvtx: need vertices") - adj_mat = utils.mesh_edges(self._geo.faces) + adj_mat = utils.mesh_edges(self.geo[hemi].faces) smooth_mat = utils.smoothing_matrix(vertices, adj_mat, smoothing_steps) else: @@ -475,23 +878,11 @@ def add_data(self, array, min=None, max=None, thresh=None, else: raise ValueError("data has to be 1D or 2D") - if smooth_mat != None: + if smooth_mat is not None: array_plot = smooth_mat * array_plot # Copy and byteswap to deal with Mayavi bug - mlab_plot = self._prepare_data(array_plot) - - # Set up the visualization pipeline - mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, - self._geo.y, - self._geo.z, - self._geo.faces, - scalars=mlab_plot) - if thresh is not None: - if array_plot.min() >= thresh: - warn("Data min is greater than threshold.") - else: - mesh = mlab.pipeline.threshold(mesh, low=thresh) + mlab_plot = _prepare_data(array_plot) # process colormap argument if isinstance(colormap, basestring): @@ -504,59 +895,52 @@ def add_data(self, array, min=None, max=None, thresh=None, raise ValueError(err) colormap = "blue-red" - surf = mlab.pipeline.surface(mesh, colormap=colormap, - vmin=min, vmax=max, - opacity=float(alpha)) - - # apply look up table if given - if lut is not None: - surf.module_manager.scalar_lut_manager.lut.table = lut - - # Get the original colormap table - orig_ctable = \ - surf.module_manager.scalar_lut_manager.lut.table.to_array().copy() - - # Fill in the data dict - self.data = dict(surface=surf, orig_ctable=orig_ctable, - array=array, smoothing_steps=smoothing_steps, - fmin=min, fmid=(min + max) / 2, fmax=max, - transparent=False, time=0, time_idx=0) - if vertices != None: - self.data["vertices"] = vertices - self.data["smooth_mat"] = smooth_mat - - # Get the colorbar - if colorbar: - bar = mlab.scalarbar(surf) - self._format_cbar_text(bar) - bar.scalar_bar_representation.position2 = .8, 0.09 - self.data['colorbar'] = bar - - - - # view = None on testing backend - if mlab.options.backend != 'test': - mlab.view(*view) + data = dict(array=array, smoothing_steps=smoothing_steps, + fmin=min, fmid=(min + max) / 2, fmax=max, + transparent=False, time=0, time_idx=0, + vertices=vertices, smooth_mat=smooth_mat) # Create time array and add label if 2D if array.ndim == 2: - if time == None: + if time is None: time = np.arange(array.shape[1]) self._times = time - self.data["time_label"] = time_label - self.data["time"] = time - self.data["time_idx"] = 0 + self.n_times = array.shape[1] + if not self.n_times == len(time): + raise ValueError('time is not the same length as ' + 'array.shape[1]') + data["time_label"] = time_label + data["time"] = time + data["time_idx"] = 0 y_txt = 0.05 + 0.05 * bool(colorbar) - if time_label is not None: - self.add_text(0.05, y_txt, time_label % time[0], name="time_label") else: self._times = None - - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = False - - def add_annotation(self, annot, borders=True, alpha=1): + self.n_times = None + + surfs = [] + bars = [] + views = self._toggle_render(False) + for bi, brain in enumerate(self._brain_list): + if brain['hemi'] == hemi: + out = brain['brain'].add_data(array, mlab_plot, vertices, + smooth_mat, min, max, thresh, + lut, colormap, alpha, time, + time_label, colorbar) + s, ct, bar = out + surfs.append(s) + bars.append(bar) + row, col = np.unravel_index(bi, self.brain_matrix.shape) + if array.ndim == 2 and time_label is not None: + self.add_text(0.05, y_txt, time_label % time[0], + name="time_label", row=row, col=col) + self._toggle_render(True, views) + data['surfaces'] = surfs + data['colorbars'] = bars + data['orig_ctable'] = ct + self.data_dict[hemi] = data + + def add_annotation(self, annot, borders=True, alpha=1, hemi=None, + remove_existing=True): """Add an annotation file. Parameters @@ -567,85 +951,90 @@ def add_annotation(self, annot, borders=True, alpha=1): Show only borders of regions alpha : float in [0, 1] Alpha level to control opacity - + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, data must exist + for both hemispheres. + remove_existing : bool + If True (default), remove old annotations. """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - view = mlab.view() + hemis = self._check_hemis(hemi) # Figure out where the data is coming from if os.path.isfile(annot): filepath = annot - annot = os.path.basename(filepath).split('.')[1] + path = os.path.split(filepath)[0] + file_hemi, annot = os.path.basename(filepath).split('.')[:2] + if len(hemis) > 1: + if annot[:2] == 'lh.': + filepaths = [filepath, pjoin(path, 'rh' + annot[2:])] + elif annot[:2] == 'rh.': + filepaths = [pjoin(path, 'lh' + annot[2:], filepath)] + else: + raise RuntimeError('To add both hemispheres ' + 'simultaneously, filename must ' + 'begin with "lh." or "rh."') + else: + filepaths = [filepath] else: - filepath = pjoin(self.subjects_dir, - self.subject_id, - 'label', - ".".join([self.hemi, annot, 'annot'])) - if not os.path.exists(filepath): - raise ValueError('Annotation file %s does not exist' - % filepath) - - # Read in the data - labels, cmap, _ = nib.freesurfer.read_annot(filepath, orig_ids=True) - - # Maybe zero-out the non-border vertices - if borders: - n_vertices = labels.size - edges = utils.mesh_edges(self._geo.faces) - border_edges = labels[edges.row] != labels[edges.col] - show = np.zeros(n_vertices, dtype=np.int) - show[np.unique(edges.row[border_edges])] = 1 - labels *= show - - # Handle null labels properly - # (tksurfer doesn't use the alpha channel, so sometimes this - # is set weirdly. For our purposes, it should always be 0. - # Unless this sometimes causes problems? - cmap[np.where(cmap[:, 4] == 0), 3] = 0 - if np.any(labels == 0) and not np.any(cmap[:, -1] == 0): - cmap = np.vstack((cmap, np.zeros(5, int))) - - # Set label ids sensibly - ord = np.argsort(cmap[:, -1]) - ids = ord[np.searchsorted(cmap[ord, -1], labels)] - cmap = cmap[:, :4] - - # Set the alpha level - alpha_vec = cmap[:, 3] - alpha_vec[alpha_vec > 0] = alpha * 255 - - # Maybe get rid of old annot - if hasattr(self, "annot"): - self.annot['surface'].remove() - - # Create an mlab surface to visualize the annot - mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, - self._geo.y, - self._geo.z, - self._geo.faces, - scalars=ids) - surf = mlab.pipeline.surface(mesh, name=annot) - - # Set the color table - surf.module_manager.scalar_lut_manager.lut.table = cmap - - # Set the brain attributes - self.annot = dict(surface=surf, name=annot, colormap=cmap) - - # Testing backend doesn't have disable_render option, view = None - if mlab.options.backend != 'test': - mlab.view(*view) - self._f.scene.disable_render = False + filepaths = [] + for hemi in hemis: + filepath = pjoin(self.subjects_dir, + self.subject_id, + 'label', + ".".join([hemi, annot, 'annot'])) + if not os.path.exists(filepath): + raise ValueError('Annotation file %s does not exist' + % filepath) + filepaths += [filepath] + + views = self._toggle_render(False) + if remove_existing is True: + # Get rid of any old annots + for a in self.annot_list: + a['surface'].remove() + self.annot_list = [] + + al = self.annot_list + for hemi, filepath in zip(hemis, filepaths): + # Read in the data + labels, cmap, _ = nib.freesurfer.read_annot(filepath, + orig_ids=True) + + # Maybe zero-out the non-border vertices + if borders: + n_vertices = labels.size + edges = utils.mesh_edges(self.geo[hemi].faces) + border_edges = labels[edges.row] != labels[edges.col] + show = np.zeros(n_vertices, dtype=np.int) + show[np.unique(edges.row[border_edges])] = 1 + labels *= show + + # Handle null labels properly + # (tksurfer doesn't use the alpha channel, so sometimes this + # is set weirdly. For our purposes, it should always be 0. + # Unless this sometimes causes problems? + cmap[np.where(cmap[:, 4] == 0), 3] = 0 + if np.any(labels == 0) and not np.any(cmap[:, -1] == 0): + cmap = np.vstack((cmap, np.zeros(5, int))) + + # Set label ids sensibly + ord = np.argsort(cmap[:, -1]) + ids = ord[np.searchsorted(cmap[ord, -1], labels)] + cmap = cmap[:, :4] + + # Set the alpha level + alpha_vec = cmap[:, 3] + alpha_vec[alpha_vec > 0] = alpha * 255 + + for brain in self._brain_list: + if brain['hemi'] == hemi: + al.append(brain['brain'].add_annotation(annot, ids, cmap)) + self.annot_list = al + self._toggle_render(True, views) def add_label(self, label, color="crimson", alpha=1, - scalar_thresh=None, borders=False): + scalar_thresh=None, borders=False, hemi=None): """Add an ROI label to the image. Parameters @@ -664,24 +1053,16 @@ def add_label(self, label, color="crimson", alpha=1, scalar >= thresh) borders : bool show only label borders + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, an error will + be thrown. Notes ----- To remove previously added labels, run Brain.remove_labels(). - """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - view = mlab.view() - - # Figure out where the data is coming from - + hemi = self._check_hemi(hemi) if isinstance(label, basestring): if os.path.isfile(label): filepath = label @@ -691,7 +1072,7 @@ def add_label(self, label, color="crimson", alpha=1, filepath = pjoin(self.subjects_dir, self.subject_id, 'label', - ".".join([self.hemi, label_name, 'label'])) + ".".join([hemi, label_name, 'label'])) if not os.path.exists(filepath): raise ValueError('Label file %s does not exist' % filepath) @@ -704,7 +1085,7 @@ def add_label(self, label, color="crimson", alpha=1, else: # try to extract parameters from label instance try: - hemi = label.hemi + lhemi = label.hemi ids = label.vertices if label.name is None: label_name = 'unnamed' @@ -718,50 +1099,42 @@ def add_label(self, label, color="crimson", alpha=1, 'must have attributes "hemi", "vertices", ' '"name", and (if scalar_thresh is not None)' '"values"') - if not hemi == self.hemi: + if not lhemi == hemi: raise ValueError('label hemisphere (%s) and brain hemisphere ' - '(%s) must match' % (label.hemi, self.hemi)) + '(%s) must match' % (lhemi, hemi)) if scalar_thresh is not None: ids = ids[scalars >= scalar_thresh] - label = np.zeros(self._geo.coords.shape[0]) + label = np.zeros(self.geo[hemi].coords.shape[0]) label[ids] = 1 # make sure we have a unique name - if label_name in self.labels: + if label_name in self.labels_dict: i = 2 name = label_name + '_%i' - while name % i in self.labels: + while name % i in self.labels_dict: i += 1 label_name = name % i if borders: n_vertices = label.size - edges = utils.mesh_edges(self._geo.faces) + edges = utils.mesh_edges(self.geo[hemi].faces) border_edges = label[edges.row] != label[edges.col] show = np.zeros(n_vertices, dtype=np.int) show[np.unique(edges.row[border_edges])] = 1 label *= show - mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, - self._geo.y, - self._geo.z, - self._geo.faces, - scalars=label) - surf = mlab.pipeline.surface(mesh, name=label_name) - - color = colorConverter.to_rgba(color, alpha) - cmap = np.array([(0, 0, 0, 0,), color]) * 255 - surf.module_manager.scalar_lut_manager.lut.table = cmap - - self.labels[label_name] = surf - - # Testing backend doesn't have disable_render option, view = None - if mlab.options.backend != 'test': - mlab.view(*view) - self._f.scene.disable_render = False - - def remove_labels(self, labels=None): + # make a list of all the plotted labels + ll = [] + views = self._toggle_render(False) + for brain in self._brain_list: + if brain['hemi'] == hemi: + ll.append(brain['brain'].add_label(label, label_name, + color, alpha)) + self.labels_dict[label_name] = ll + self._toggle_render(True, views) + + def remove_labels(self, labels=None, hemi=None): """Remove one or more previously added labels from the image. Parameters @@ -770,17 +1143,24 @@ def remove_labels(self, labels=None): Labels to remove. Can be a string naming a single label, or None to remove all labels. Possible names can be found in the Brain.labels attribute. + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, an error will + be thrown. """ + hemi = self._check_hemi(hemi) if labels is None: - labels = self.labels.keys() + labels = self.labels_dict.keys() elif isinstance(labels, str): labels = [labels] for key in labels: - label = self.labels.pop(key) - label.remove() + label = self.labels_dict.pop(key) + for ll in label: + ll.remove() - def add_morphometry(self, measure, grayscale=False): + def add_morphometry(self, measure, grayscale=False, hemi=None, + remove_existing=True): """Add a morphometry overlay to the image. Parameters @@ -789,84 +1169,72 @@ def add_morphometry(self, measure, grayscale=False): which measure to load grayscale : bool whether to load the overlay with a grayscale colormap - + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, data must exist + for both hemispheres. + remove_existing : bool + If True (default), remove old annotations. """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - # Find the source data - surf_dir = pjoin(self.subjects_dir, self.subject_id, 'surf') - morph_file = pjoin(surf_dir, '.'.join([self.hemi, measure])) - if not os.path.exists(morph_file): - raise ValueError( - 'Could not find %s in subject directory' % morph_file) - - # Preset colormaps - cmap_dict = dict(area="pink", - curv="RdBu", - jacobian_white="pink", - sulc="RdBu", - thickness="pink") - - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - - # Maybe get rid of an old overlay - if hasattr(self, "morphometry"): - self.morphometry['surface'].remove() - self.morphometry['colorbar'].visible = False - - # Save the inital view - view = mlab.view() - - # Read in the morphometric data - morph_data = nib.freesurfer.read_morph_data(morph_file) - - # Get a cortex mask for robust range - self._geo.load_label("cortex") - ctx_idx = self._geo.labels["cortex"] - - # Get the display range - if measure == "thickness": - min, max = 1, 4 - else: - min, max = stats.describe(morph_data[ctx_idx])[1] - - # Set up the Mayavi pipeline - morph_data = self._prepare_data(morph_data) + hemis = self._check_hemis(hemi) + morph_files = [] + for hemi in hemis: + # Find the source data + surf_dir = pjoin(self.subjects_dir, self.subject_id, 'surf') + morph_file = pjoin(surf_dir, '.'.join([hemi, measure])) + if not os.path.exists(morph_file): + raise ValueError( + 'Could not find %s in subject directory' % morph_file) + morph_files += [morph_file] + + views = self._toggle_render(False) + if remove_existing is True: + # Get rid of any old overlays + for m in self.morphometry_list: + m['surface'].remove() + m['colorbar'].visible = False + self.morphometry_list = [] + ml = self.morphometry_list + for hemi, morph_file in zip(hemis, morph_files): + # Preset colormaps + cmap_dict = dict(area="pink", + curv="RdBu", + jacobian_white="pink", + sulc="RdBu", + thickness="pink") + + # Read in the morphometric data + morph_data = nib.freesurfer.read_morph_data(morph_file) + + # Get a cortex mask for robust range + self.geo[hemi].load_label("cortex") + ctx_idx = self.geo[hemi].labels["cortex"] + + # Get the display range + if measure == "thickness": + min, max = 1, 4 + else: + min, max = stats.describe(morph_data[ctx_idx])[1] - mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, - self._geo.y, - self._geo.z, - self._geo.faces, - scalars=morph_data) - if grayscale: - colormap = "gray" - else: - colormap = cmap_dict[measure] - surf = mlab.pipeline.surface(mesh, colormap=colormap, - vmin=min, vmax=max, - name=measure) + # Set up the Mayavi pipeline + morph_data = _prepare_data(morph_data) - # Get the colorbar - bar = mlab.scalarbar(surf) - self._format_cbar_text(bar) - bar.scalar_bar_representation.position2 = .8, 0.09 + if grayscale: + colormap = "gray" + else: + colormap = cmap_dict[measure] - # Fil in the morphometry dict - self.morphometry = dict(surface=surf, - colorbar=bar, - measure=measure) - # Testing backend doesn't have disable_render option, view = None - if mlab.options.backend != 'test': - mlab.view(*view) - self._f.scene.disable_render = False + for brain in self._brain_list: + if brain['hemi'] == hemi: + ml.append(brain['brain'].add_morphometry(morph_data, + colormap, measure, + min, max)) + self.morphometry_list = ml + self._toggle_render(True, views) def add_foci(self, coords, coords_as_verts=False, map_surface=None, - scale_factor=1, color="white", alpha=1, name=None): + scale_factor=1, color="white", alpha=1, name=None, + hemi=None): """Add spherical foci, possibly mapping to displayed surf. The foci spheres can be displayed at the coordinates given, or @@ -891,55 +1259,47 @@ def add_foci(self, coords, coords_as_verts=False, map_surface=None, opacity of focus gylphs name : str internal name to use - + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, an error will + be thrown. """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab + hemi = self._check_hemi(hemi) # Figure out how to interpret the first parameter if coords_as_verts: - coords = self._geo.coords[coords] + coords = self.geo[hemi].coords[coords] map_surface = None # Possibly map the foci coords through a surface if map_surface is None: foci_coords = np.atleast_2d(coords) else: - foci_surf = io.Surface(self.subject_id, self.hemi, map_surface, + foci_surf = io.Surface(self.subject_id, hemi, map_surface, subjects_dir=self.subjects_dir) foci_surf.load_geometry() foci_vtxs = utils.find_closest_vertices(foci_surf.coords, coords) - foci_coords = self._geo.coords[foci_vtxs] + foci_coords = self.geo[hemi].coords[foci_vtxs] # Get a unique name (maybe should take this approach elsewhere) if name is None: - name = "foci_%d" % (len(self.foci) + 1) + name = "foci_%d" % (len(self.foci_dict) + 1) # Convert the color code if not isinstance(color, tuple): color = colorConverter.to_rgb(color) - # Create the visualization - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - view = mlab.view() - points = mlab.points3d(foci_coords[:, 0], - foci_coords[:, 1], - foci_coords[:, 2], - np.ones(foci_coords.shape[0]), - scale_factor=(10. * scale_factor), - color=color, opacity=alpha, name=name) - self.foci[name] = points - # Testing backend doesn't have disable_render option, view = None - if mlab.options.backend != 'test': - mlab.view(*view) - self._f.scene.disable_render = False + views = self._toggle_render(False) + fl = [] + for brain in self._brain_list: + if brain['hemi'] == hemi: + fl.append(brain['brain'].add_foci(foci_coords, scale_factor, + color, alpha, name)) + self.foci_dict[name] = fl + self._toggle_render(True, views) def add_contour_overlay(self, source, min=None, max=None, - n_contours=7, line_width=1.5): + n_contours=7, line_width=1.5, hemi=None): """Add a topographic contour overlay of the positive data. Note: This visualization will look best when using the "low_contrast" @@ -957,58 +1317,39 @@ def add_contour_overlay(self, source, min=None, max=None, number of contours to use in the display line_width : float width of contour lines - + hemi : str | None + If None, it is assumed to belong to the hemipshere being + shown. If two hemispheres are being shown, an error will + be thrown. """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab + hemi = self._check_hemi(hemi) # Read the scalar data - scalar_data, _ = self._read_scalar_data(source) - + scalar_data, _ = self._read_scalar_data(source, hemi) min, max = self._get_display_range(scalar_data, min, max, "pos") - # Prep the viz - # Testing backend doesn't have this option - if mlab.options.backend != 'test': - self._f.scene.disable_render = True - view = mlab.view() + # Deal with Mayavi bug + scalar_data = _prepare_data(scalar_data) # Maybe get rid of an old overlay if hasattr(self, "contour"): - self.contour['surface'].remove() - self.contour['colorbar'].visible = False - - # Deal with Mayavi bug - scalar_data = self._prepare_data(scalar_data) - - # Set up the pipeline - mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, self._geo.y, - self._geo.z, self._geo.faces, - scalars=scalar_data) - thresh = mlab.pipeline.threshold(mesh, low=min) - surf = mlab.pipeline.contour_surface(thresh, contours=n_contours, - line_width=line_width) - - # Set the colorbar and range correctly - bar = mlab.scalarbar(surf, - nb_colors=n_contours, - nb_labels=n_contours + 1) - bar.data_range = min, max - self._format_cbar_text(bar) - bar.scalar_bar_representation.position2 = .8, 0.09 - - # Set up a dict attribute with pointers at important things - self.contour = dict(surface=surf, colorbar=bar) - - # Show the new overlay - # Testing backend doesn't have disable_render option, view = None - if mlab.options.backend != 'test': - mlab.view(*view) - self._f.scene.disable_render = False - - def add_text(self, x, y, text, name, color=None, opacity=1.0): + for c in self.contour_list: + c['surface'].remove() + c['colorbar'].visible = False + + views = self._toggle_render(False) + cl = [] + for brain in self._brain_list: + if brain['hemi'] == hemi: + cl.append(brain['brain'].add_contour_overlay(scalar_data, + min, max, + n_contours, + line_width)) + self.contour_list = cl + self._toggle_render(True, views) + + def add_text(self, x, y, text, name, color=None, opacity=1.0, + row=-1, col=-1): """ Add a text to the visualization Parameters @@ -1025,200 +1366,491 @@ def add_text(self, x, y, text, name, color=None, opacity=1.0): Color of the text. Default: (1, 1, 1) opacity : Float Opacity of the text. Default: 1.0 + row : int + Row index of which brain to use + col : int + Column index of which brain to use """ - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - - text = mlab.text(x, y, text, figure=None, name=name, - color=color, opacity=opacity) + if name in self.texts_dict: + self.texts_dict[name]['text'].remove() + text = self.brain_matrix[row, col].add_text(x, y, text, + name, color, opacity) + self.texts_dict[name] = dict(row=row, col=col, text=text) - self.texts[name] = text - - def _set_scene_properties(self, config_opts): - """Set the scene_prop dict from the config parser. + def update_text(self, text, name, row=-1, col=-1): + """Update text label Parameters ---------- - config_opts : dict - dictionary of config file "visual" options - + text : str + New text for label + name : str + Name of text label """ - colors = dict(black=[0, 0, 0], - white=[256, 256, 256], - midnight=[12, 7, 32], - slate=[112, 128, 144], - charcoal=[59, 69, 79], - sand=[245, 222, 179]) - - try: - bg_color_name = config_opts['background'] - except KeyError: - bg_color_name = config.get("visual", "background") - if bg_color_name is not None: - bg_color_code = colorConverter.to_rgb(bg_color_name) - else: - bg_color_code = None - - try: - fg_color_name = config_opts['foreground'] - except KeyError: - fg_color_name = config.get("visual", "foreground") - fg_color_code = colorConverter.to_rgb(fg_color_name) + if name not in self.texts_dict: + raise KeyError('text name "%s" unknown' % name) + self.texts_dict[name]['text'].text = text + + ########################################################################### + # DATA SCALING / DISPLAY + def reset_view(self): + """Orient camera to display original view + """ + for view, brain in zip(self._original_views, self._brain_list): + brain['brain'].show_view(view) - try: - size = config_opts['size'] - except KeyError: - size = config.getfloat("visual", "size") - size = (size, size) + def show_view(self, view=None, roll=None, distance=None, row=-1, col=-1): + """Orient camera to display view - self.scene_properties = dict(fgcolor=fg_color_code, - bgcolor=bg_color_code, - size=size) + Parameters + ---------- + view : {'lateral' | 'medial' | 'rostral' | 'caudal' | + 'dorsal' | 'ventral' | 'frontal' | 'parietal' | + dict} + brain surface to view or kwargs to pass to mlab.view() - def _orient_lights(self): - """Set lights to come from same direction relative to brain.""" - if self.hemi == "rh": - if self._f.scene is not None: - for light in self._f.scene.light_manager.lights: - light.azimuth *= -1 + Returns + ------- + view : tuple + tuple returned from mlab.view + roll : float + camera roll + distance : float | 'auto' | None + distance from the origin + row : int + Row index of which brain to use + col : int + Column index of which brain to use + """ + return self.brain_matrix[row][col].show_view(view, roll, distance) - def _get_geo_colors(self, config_opts): - """Return an mlab colormap name, vmin, and vmax for binary curvature. + def set_distance(self, distance=None): + """Set view distances for all brain plots to the same value Parameters ---------- - config_opts : dict - dictionary of config file "visual" options + distance : float | None + Distance to use. If None, brains are set to the farthest + "best fit" distance across all current views; note that + the underlying "best fit" function can be buggy. Returns ------- - colormap : string - mlab colormap name - vmin : float - curv colormap minimum - vmax : float - curv colormap maximum - reverse : boolean - boolean indicating whether the colormap should be reversed + distance : float + The distance used. + """ + try: + from mayavi import mlab + assert mlab + except: + from enthought.mayavi import mlab + if distance is None: + distance = [] + for ff in self._figures: + for f in ff: + mlab.view(figure=f, distance='auto') + v = mlab.view(figure=f) + # This should only happen for the test backend + if v is None: + v = [0, 0, 100] + distance += [v[2]] + distance = max(distance) + + for ff in self._figures: + for f in ff: + mlab.view(distance=distance, figure=f) + return distance + + @verbose + def scale_data_colormap(self, fmin, fmid, fmax, transparent, verbose=None): + """Scale the data colormap. + Parameters + ---------- + fmin : float + minimum value of colormap + fmid : float + value corresponding to color midpoint + fmax : float + maximum value for colormap + transparent : boolean + if True: use a linear transparency between fmin and fmid + verbose : bool, str, int, or None + If not None, override default verbose level (see surfer.verbose). """ - colormap_map = dict(classic=("Greys", -1, 2, False), - high_contrast=("Greys", -.1, 1.3, False), - low_contrast=("Greys", -5, 5, False), - bone=("bone", -.2, 2, True)) + if not (fmin < fmid) and (fmid < fmax): + raise ValueError("Invalid colormap, we need fmin= self.n_times: + raise ValueError("time index out of range") - ftype = fname[fname.rfind('.') + 1:] - good_ftypes = ['png', 'jpg', 'bmp', 'tiff', 'ps', - 'eps', 'pdf', 'rib', 'oogl', 'iv', 'vrml', 'obj'] - if not ftype in good_ftypes: - raise ValueError("Supported image types are %s" - % " ".join(good_ftypes)) - mlab.draw(self._f) - mlab.savefig(fname, figure=self._f) + views = self._toggle_render(False) + for hemi in ['lh', 'rh']: + data = self.data_dict[hemi] + if data is not None: + plot_data = data["array"][:, time_idx] + if data["smooth_mat"] is not None: + plot_data = data["smooth_mat"] * plot_data + for surf in data["surfaces"]: + surf.mlab_source.scalars = plot_data + data["time_idx"] = time_idx + + # Update time label + if data["time_label"]: + time = data["time"][time_idx] + self.update_text(data["time_label"] % time, "time_label") + self._toggle_render(True, views) + + @verbose + def set_data_smoothing_steps(self, smoothing_steps, verbose=None): + """Set the number of smoothing steps - def screenshot(self, mode='rgb', antialiased=False): - """Generate a screenshot of current view + Parameters + ---------- + smoothing_steps : int + Number of smoothing steps + verbose : bool, str, int, or None + If not None, override default verbose level (see surfer.verbose). + """ + views = self._toggle_render(False) + for hemi in ['lh', 'rh']: + data = self.data_dict[hemi] + if data is not None: + adj_mat = utils.mesh_edges(self.geo[hemi].faces) + smooth_mat = utils.smoothing_matrix(data["vertices"], + adj_mat, smoothing_steps) + data["smooth_mat"] = smooth_mat + + # Redraw + if data["array"].ndim == 1: + plot_data = data["array"] + else: + plot_data = data["array"][:, data["time_idx"]] - Wraps to mlab.screenshot for ease of use. + plot_data = data["smooth_mat"] * plot_data + for surf in data["surfaces"]: + surf.mlab_source.scalars = plot_data + + # Update data properties + data["smoothing_steps"] = smoothing_steps + self._toggle_render(True, views) + + def set_time(self, time): + """Set the data time index to the time point closest to time Parameters ---------- - mode: string + time : scalar + Time. + """ + if self.n_times is None: + raise RuntimeError("Brain has no time axis") + times = self._times + + # Check that time is in range + tmin = np.min(times) + tmax = np.max(times) + max_diff = (tmax - tmin) / (len(times) - 1) / 2 + if time < tmin - max_diff or time > tmax + max_diff: + err = ("time = %s lies outside of the time axis " + "[%s, %s]" % (time, tmin, tmax)) + raise ValueError(err) + + idx = np.argmin(np.abs(times - time)) + self.set_data_time_index(idx) + + def _get_colorbars(self, row, col): + shape = self.brain_matrix.shape + row = row % shape[0] + col = col % shape[1] + ind = np.ravel_multi_index((row, col), self.brain_matrix.shape) + colorbars = [] + h = self._brain_list[ind]['hemi'] + if self.data_dict[h] is not None and 'colorbars' in self.data_dict[h]: + colorbars.append(self.data_dict[h]['colorbars'][row]) + if len(self.morphometry_list) > 0: + colorbars.append(self.morphometry_list[ind]['colorbar']) + if len(self.contour_list) > 0: + colorbars.append(self.contour_list[ind]['colorbar']) + if len(self.overlays_dict) > 0: + for name, obj in self.overlays_dict.items(): + for bar in ["pos_bar", "neg_bar"]: + try: + colorbars.append(getattr(obj[ind], bar)) + except AttributeError: + pass + return colorbars + + def _colorbar_visibility(self, visible, row, col): + for cb in self._get_colorbars(row, col): + cb.visible = visible + + def show_colorbar(self, row=-1, col=-1): + """Show colorbar(s) for given plot + + Parameters + ---------- + row : int + Row index of which brain to use + col : int + Column index of which brain to use + """ + self._colorbar_visibility(True, row, col) + + def hide_colorbar(self, row=-1, col=-1): + """Hide colorbar(s) for given plot + + Parameters + ---------- + row : int + Row index of which brain to use + col : int + Column index of which brain to use + """ + self._colorbar_visibility(False, row, col) + + def close(self): + """Close all figures and cleanup data structure.""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + for ri, ff in enumerate(self._figures): + for ci, f in enumerate(ff): + if f is not None: + mlab.close(f) + self._figures[ri][ci] = None + + #should we tear down other variables? + if self._v is not None: + self._v.dispose() + self._v = None + + def __del__(self): + if hasattr(self, '_v') and self._v is not None: + self._v.dispose() + self._v = None + + + ########################################################################### + # SAVING OUTPUT + def save_single_image(self, filename, row=-1, col=-1): + """Save view from one panel to disk + + Only mayavi image types are supported: + (png jpg bmp tiff ps eps pdf rib oogl iv vrml obj + + Parameters + ---------- + filename: string + path to new image file + row : int + row index of the brain to use + col : int + column index of the brain to use + + Due to limitations in TraitsUI, if multiple views or hemi='split' + is used, there is no guarantee painting of the windows will + complete before control is returned to the command line. Thus + we strongly recommend using only one figure window (which uses + a Mayavi figure to plot instead of TraitsUI) if you intend to + script plotting commands. + """ + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + brain = self.brain_matrix[row, col] + ftype = filename[filename.rfind('.') + 1:] + good_ftypes = ['png', 'jpg', 'bmp', 'tiff', 'ps', + 'eps', 'pdf', 'rib', 'oogl', 'iv', 'vrml', 'obj'] + if not ftype in good_ftypes: + raise ValueError("Supported image types are %s" + % " ".join(good_ftypes)) + mlab.draw(brain._f) + mlab.savefig(filename, figure=brain._f) + + def save_image(self, filename): + """Save view from all panels to disk + + Only mayavi image types are supported: + (png jpg bmp tiff ps eps pdf rib oogl iv vrml obj + + Parameters + ---------- + filename: string + path to new image file + + Due to limitations in TraitsUI, if multiple views or hemi='split' + is used, there is no guarantee painting of the windows will + complete before control is returned to the command line. Thus + we strongly recommend using only one figure window (which uses + a Mayavi figure to plot instead of TraitsUI) if you intend to + script plotting commands. + """ + misc.imsave(filename, self.screenshot()) + + def screenshot(self, mode='rgb', antialiased=False): + """Generate a screenshot of current view + + Wraps to mlab.screenshot for ease of use. + + Parameters + ---------- + mode: string + Either 'rgb' or 'rgba' for values to return + antialiased: bool + Antialias the image (see mlab.screenshot() for details) + row : int + row index of the brain to use + col : int + column index of the brain to use + + Returns + ------- + screenshot: array + Image pixel values + + Notes + ----- + Due to limitations in TraitsUI, if multiple views or hemi='split' + is used, there is no guarantee painting of the windows will + complete before control is returned to the command line. Thus + we strongly recommend using only one figure window (which uses + a Mayavi figure to plot instead of TraitsUI) if you intend to + script plotting commands. + """ + row = [] + for ri in range(self.brain_matrix.shape[0]): + col = [] + n_col = 2 if self._hemi == 'split' else 1 + for ci in range(n_col): + col += [self.screenshot_single(mode, antialiased, + ri, ci)] + row += [np.concatenate(col, axis=1)] + data = np.concatenate(row, axis=0) + return data + + def screenshot_single(self, mode='rgb', antialiased=False, row=-1, col=-1): + """Generate a screenshot of current view from a single panel + + Wraps to mlab.screenshot for ease of use. + + Parameters + ---------- + mode: string Either 'rgb' or 'rgba' for values to return antialiased: bool Antialias the image (see mlab.screenshot() for details) + row : int + row index of the brain to use + col : int + column index of the brain to use Returns ------- screenshot: array Image pixel values + + Notes + ----- + Due to limitations in TraitsUI, if multiple views or hemi='split' + is used, there is no guarantee painting of the windows will + complete before control is returned to the command line. Thus + we strongly recommend using only one figure window (which uses + a Mayavi figure to plot instead of TraitsUI) if you intend to + script plotting commands. """ try: from mayavi import mlab + assert mlab except ImportError: from enthought.mayavi import mlab - return mlab.screenshot(self._f, mode, antialiased) + brain = self.brain_matrix[row, col] + return mlab.screenshot(brain._f, mode, antialiased) - def save_imageset(self, prefix, views, filetype='png', colorbar='auto'): + def save_imageset(self, prefix, views, filetype='png', colorbar='auto', + row=-1, col=-1): """Convenience wrapper for save_image Files created are prefix+'_$view'+filetype Parameters ---------- - prefix: string - filename prefix for image to be created + prefix: string | None + filename prefix for image to be created. If None, a list of + arrays representing images is returned (not saved to disk). views: list desired views for images filetype: string @@ -1227,6 +1859,10 @@ def save_imageset(self, prefix, views, filetype='png', colorbar='auto'): if None no colorbar is visible. If 'auto' is given the colorbar is only shown in the middle view. Otherwise on the listed views when a list of int is passed. + row : int + row index of the brain to use + col : int + column index of the brain to use Returns ------- @@ -1241,23 +1877,24 @@ def save_imageset(self, prefix, views, filetype='png', colorbar='auto'): images_written = [] for iview, view in enumerate(views): try: - fname = "%s_%s.%s" % (prefix, view, filetype) - images_written.append(fname) if colorbar is not None and iview in colorbar: - self.show_colorbar() + self.show_colorbar(row, col) else: - self.hide_colorbar() - self.show_view(view) - - try: - self.save_image(fname) - except ValueError: - print("Bad image type") + self.hide_colorbar(row, col) + self.show_view(view, row=row, col=col) + if prefix is not None: + fname = "%s_%s.%s" % (prefix, view, filetype) + images_written.append(fname) + self.save_single_image(fname, row, col) + else: + images_written.append(self.screenshot_single(row=row, + col=col)) except ValueError: print("Skipping %s: not in view dict" % view) return images_written - def save_image_sequence(self, time_idx, fname_pattern, use_abs_idx=True): + def save_image_sequence(self, time_idx, fname_pattern, use_abs_idx=True, + row=-1, col=-1): """Save a temporal image sequence The files saved are named "fname_pattern % (pos)" where "pos" is a @@ -1273,21 +1910,23 @@ def save_image_sequence(self, time_idx, fname_pattern, use_abs_idx=True): if True the indices given by "time_idx" are used in the filename if False the index in the filename starts at zero and is incremented by one for each image (Default: True) + row : int + row index of the brain to use + col : int + column index of the brain to use Returns ------- images_written: list all filenames written """ - current_time_idx = self.data["time_idx"] - images_written = list() rel_pos = 0 for idx in time_idx: self.set_data_time_index(idx) fname = fname_pattern % (idx if use_abs_idx else rel_pos) - self.save_image(fname) + self.save_single_image(fname, row, col) images_written.append(fname) rel_pos += 1 @@ -1296,82 +1935,15 @@ def save_image_sequence(self, time_idx, fname_pattern, use_abs_idx=True): return images_written - def scale_data_colormap(self, fmin, fmid, fmax, transparent): - """Scale the data colormap. - - Parameters - ---------- - fmin : float - minimum value of colormap - fmid : float - value corresponding to color midpoint - fmax : float - maximum value for colormap - transparent : boolean - if True: use a linear transparency between fmin and fmid - """ - - if not (fmin < fmid) and (fmid < fmax): - raise ValueError("Invalid colormap, we need fmin= self.data["array"].shape[1]: - raise ValueError("time index out of range") - - plot_data = self.data["array"][:, time_idx] - - if "smooth_mat" in self.data: - plot_data = self.data["smooth_mat"] * plot_data - self.data["surface"].mlab_source.scalars = plot_data - self.data["time_idx"] = time_idx - - # Update time label - if self.data["time_label"]: - time = self.data["time"][time_idx] - self.update_text(self.data["time_label"] % time, "time_label") - - def set_data_smoothing_steps(self, smoothing_steps): - """ Set the number of smoothing steps - - Parameters - ---------- - smoothing_steps : int - Number of smoothing steps - """ - - adj_mat = utils.mesh_edges(self._geo.faces) - smooth_mat = utils.smoothing_matrix(self.data["vertices"], adj_mat, - smoothing_steps) - - self.data["smooth_mat"] = smooth_mat - - # Redraw - if self.data["array"].ndim == 1: - plot_data = self.data["array"] - else: - plot_data = self.data["array"][:, self.data["time_idx"]] - - plot_data = self.data["smooth_mat"] * plot_data - - self.data["surface"].mlab_source.scalars = plot_data - - # Update data properties - self.data["smoothing_steps"] = smoothing_steps - - def set_time(self, time): - """Set the data time index to the time point closest to time - - Parameters - ---------- - time : scalar - Time. - """ - times = getattr(self, '_times', None) - if times is None: - raise RuntimeError("Brain has no time axis") - - # Check that time is in range - tmin = np.min(times) - tmax = np.max(times) - max_diff = (tmax - tmin) / (len(times) - 1) / 2 - if time < tmin - max_diff or time > tmax + max_diff: - err = ("time = %s lies outside of the time axis " - "[%s, %s]" % (time, tmin, tmax)) - raise ValueError(err) - - idx = np.argmin(np.abs(times - time)) - self.set_data_time_index(idx) - - def update_text(self, text, name): - """ Update text label - - Parameters - ---------- - text : str - New text for label - name : str - Name of text label - """ - self.texts[name].text = text - - def min_diff(self, beg, end): - """Determine minimum "camera distance" between two views. - - Parameters - ---------- - beg: string - origin anatomical view - end: string - destination anatomical view - - Returns - ------- - diffs: tuple - (min view "distance", min roll "distance") - - """ - beg = self.xfm_view(beg) - end = self.xfm_view(end) - if beg == end: - dv = [360., 0.] - dr = 0 - else: - end_d = self.xfm_view(end, 'd') - beg_d = self.xfm_view(beg, 'd') - dv = [] - for b, e in zip(beg_d['v'], end_d['v']): - diff = e - b - # to minimize the rotation we need -180 <= diff <= 180 - if diff > 180: - dv.append(diff - 360) - elif diff < -180: - dv.append(diff + 360) - else: - dv.append(diff) - dr = np.array(end_d['r']) - np.array(beg_d['r']) - return (np.array(dv), dr) - - def animate(self, views, n_steps=180., fname=None, use_cache=False): + def animate(self, views, n_steps=180., fname=None, use_cache=False, + row=-1, col=-1): """Animate a rotation. Currently only rotations through the axial plane are allowed. @@ -1549,8 +2010,13 @@ def animate(self, views, n_steps=180., fname=None, use_cache=False): fname should end in '.avi' as only the AVI format is supported use_cache: bool Use previously generated images in ./.tmp/ + row : int + Row index of the brain to use + col : int + Column index of the brain to use """ - gviews = map(self.xfm_view, views) + brain = self.brain_matrix[row, col] + gviews = map(brain._xfm_view, views) allowed = ('lateral', 'caudal', 'medial', 'rostral') if not len([v for v in gviews if v in allowed]) == len(gviews): raise ValueError('Animate through %s views.' % ' '.join(allowed)) @@ -1564,19 +2030,19 @@ def animate(self, views, n_steps=180., fname=None, use_cache=False): for i, beg in enumerate(gviews): try: end = gviews[i + 1] - dv, dr = self.min_diff(beg, end) + dv, dr = brain._min_diff(beg, end) dv /= np.array((n_steps)) dr /= np.array((n_steps)) - self.show_view(beg) + brain.show_view(beg) for i in range(int(n_steps)): - self._f.scene.camera.orthogonalize_view_up() - self._f.scene.camera.azimuth(dv[0]) - self._f.scene.camera.elevation(dv[1]) - self._f.scene.renderer.reset_camera_clipping_range() - self._f.scene.render() + brain._f.scene.camera.orthogonalize_view_up() + brain._f.scene.camera.azimuth(dv[0]) + brain._f.scene.camera.elevation(dv[1]) + brain._f.scene.renderer.reset_camera_clipping_range() + _force_render([[brain._f]], self._window_backend) if fname is not None: if not (os.path.isfile(tmp_fname % i) and use_cache): - self.save_image(tmp_fname % i) + self.save_single_image(tmp_fname % i, row, col) except IndexError: pass if fname is not None: @@ -1595,7 +2061,86 @@ def animate(self, views, n_steps=180., fname=None, use_cache=False): if ret: print("\n\nError occured when exporting movie\n\n") - def xfm_view(self, view, out='s'): + +class _Hemisphere(object): + """Object for visualizing one hemisphere with mlab""" + def __init__(self, subject_id, hemi, surf, figure, geo, curv, title, + config_opts, subjects_dir, bg_color, offset, backend): + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + if not hemi in ['lh', 'rh']: + raise ValueError('hemi must be either "lh" or "rh"') + # Set the identifying info + self.subject_id = subject_id + self.hemi = hemi + self.subjects_dir = subjects_dir + self.viewdict = viewdicts[hemi] + self.surf = surf + self._f = figure + self._bg_color = bg_color + self._backend = backend + + # mlab pipeline mesh and surface for geomtery + self._geo = geo + if curv: + curv_data = self._geo.bin_curv + meshargs = dict(scalars=curv_data) + colormap, vmin, vmax, reverse = self._get_geo_colors(config_opts) + kwargs = dict(colormap=colormap, vmin=vmin, vmax=vmax) + else: + curv_data = None + meshargs = dict() + kwargs = dict(color=(.5, .5, .5)) + meshargs['figure'] = self._f + x, y, z, f = self._geo.x, self._geo.y, self._geo.z, self._geo.faces + self._geo_mesh = mlab.pipeline.triangular_mesh_source(x, y, z, f, + **meshargs) + self._geo_surf = mlab.pipeline.surface(self._geo_mesh, + figure=self._f, reset_zoom=True, + **kwargs) + if curv and reverse: + curv_bar = mlab.scalarbar(self._geo_surf) + curv_bar.reverse_lut = True + curv_bar.visible = False + + def show_view(self, view=None, roll=None, distance=None): + """Orient camera to display view""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + if isinstance(view, basestring): + try: + vd = self._xfm_view(view, 'd') + view = dict(azimuth=vd['v'][0], elevation=vd['v'][1]) + roll = vd['r'] + except ValueError as v: + print(v) + raise + + _force_render(self._f, self._backend) + if view is not None: + view['reset_roll'] = True + view['figure'] = self._f + view['distance'] = distance + # DO NOT set focal point, can screw up non-centered brains + #view['focalpoint'] = (0.0, 0.0, 0.0) + mlab.view(**view) + if roll is not None: + mlab.roll(roll=roll, figure=self._f) + _force_render(self._f, self._backend) + + view = mlab.view(figure=self._f) + roll = mlab.roll(figure=self._f) + + return view, roll + + def _xfm_view(self, view, out='s'): """Normalize a given string to available view Parameters @@ -1624,110 +2169,286 @@ def xfm_view(self, view, out='s'): else: return view - def _get_colorbars(self): - colorbars = [] - if hasattr(self, 'data') and 'colorbar' in self.data: - colorbars.append(self.data['colorbar']) - if hasattr(self, 'morphometry') and 'colorbar' in self.morphometry: - colorbars.append(self.morphometry['colorbar']) - if hasattr(self, 'contour') and 'colorbar' in self.contour: - colorbars.append(self.contour['colorbar']) - if hasattr(self, 'overlays'): - for name, obj in self.overlays.items(): - for bar in ["pos_bar", "neg_bar"]: - try: - colorbars.append(getattr(obj, bar)) - except AttributeError: - pass - return colorbars + def _min_diff(self, beg, end): + """Determine minimum "camera distance" between two views. - def _colorbar_visibility(self, visible): - for cb in self._get_colorbars(): - cb.visible = visible + Parameters + ---------- + beg: string + origin anatomical view + end: string + destination anatomical view - def show_colorbar(self): - "Show colorbar(s)" - self._colorbar_visibility(True) + Returns + ------- + diffs: tuple + (min view "distance", min roll "distance") + + """ + beg = self._xfm_view(beg) + end = self._xfm_view(end) + if beg == end: + dv = [360., 0.] + dr = 0 + else: + end_d = self._xfm_view(end, 'd') + beg_d = self._xfm_view(beg, 'd') + dv = [] + for b, e in zip(beg_d['v'], end_d['v']): + diff = e - b + # to minimize the rotation we need -180 <= diff <= 180 + if diff > 180: + dv.append(diff - 360) + elif diff < -180: + dv.append(diff + 360) + else: + dv.append(diff) + dr = np.array(end_d['r']) - np.array(beg_d['r']) + return (np.array(dv), dr) - def hide_colorbar(self): - "Hide colorbar(s)" - self._colorbar_visibility(False) + def add_overlay(self, old): + """Add an overlay to the overlay dict from a file or array""" + surf = OverlayDisplay(old, figure=self._f) + for bar in ["pos_bar", "neg_bar"]: + try: + self._format_cbar_text(getattr(surf, bar)) + except AttributeError: + pass + return surf + + @verbose + def add_data(self, array, mlab_plot, vertices, smooth_mat, min, max, + thresh, lut, colormap, alpha, time, time_label, colorbar): + """Add data to the brain""" + # Calculate initial data to plot + if array.ndim == 1: + array_plot = array + elif array.ndim == 2: + array_plot = array[:, 0] + else: + raise ValueError("data has to be 1D or 2D") - def close(self): - """Close the figure and cleanup data structure.""" try: from mayavi import mlab + assert mlab except ImportError: from enthought.mayavi import mlab + # Set up the visualization pipeline + mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, + self._geo.y, + self._geo.z, + self._geo.faces, + scalars=mlab_plot, + figure=self._f) + if thresh is not None: + if array_plot.min() >= thresh: + warn("Data min is greater than threshold.") + else: + mesh = mlab.pipeline.threshold(mesh, low=thresh) - mlab.close(self._f) - #should we tear down other variables? + surf = mlab.pipeline.surface(mesh, colormap=colormap, + vmin=min, vmax=max, + opacity=float(alpha), figure=self._f) - def _get_display_range(self, scalar_data, min, max, sign): + # apply look up table if given + if lut is not None: + surf.module_manager.scalar_lut_manager.lut.table = lut - if scalar_data.min() >= 0: - sign = "pos" - elif scalar_data.max() <= 0: - sign = "neg" - self.sign = sign + # Get the original colormap table + orig_ctable = \ + surf.module_manager.scalar_lut_manager.lut.table.to_array().copy() - # Get data with a range that will make sense for automatic thresholding - if sign == "neg": - range_data = np.abs(scalar_data[np.where(scalar_data < 0)]) - elif sign == "pos": - range_data = scalar_data[np.where(scalar_data > 0)] + # Get the colorbar + if colorbar: + bar = mlab.scalarbar(surf) + self._format_cbar_text(bar) + bar.scalar_bar_representation.position2 = .8, 0.09 else: - range_data = np.abs(scalar_data) + bar = None - # Get the min and max from among various places - if min is None: - try: - min = config.getfloat("overlay", "min_thresh") - except ValueError: - min_str = config.get("overlay", "min_thresh") - if min_str == "robust_min": - min = stats.scoreatpercentile(range_data, 2) - elif min_str == "actual_min": - min = range_data.min() - else: - min = 2.0 - warn("The 'min_thresh' value in your config value must be " - "a float, 'robust_min', or 'actual_min', but it is %s. " - "I'm setting the overlay min to the config default of 2" % min) + return surf, orig_ctable, bar - if max is None: - try: - max = config.getfloat("overlay", "max_thresh") - except ValueError: - max_str = config.get("overlay", "max_thresh") - if max_str == "robust_max": - max = stats.scoreatpercentile(scalar_data, 98) - elif max_str == "actual_max": - max = range_data.max() - else: - max = stats.scoreatpercentile(range_data, 98) - warn("The 'max_thresh' value in your config value must be " - "a float, 'robust_min', or 'actual_min', but it is %s. " - "I'm setting the overlay min to the config default " - "of robust_max" % max) + def add_annotation(self, annot, ids, cmap): + """Add an annotation file""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab - return min, max + # Create an mlab surface to visualize the annot + mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, + self._geo.y, + self._geo.z, + self._geo.faces, + scalars=ids, + figure=self._f) + surf = mlab.pipeline.surface(mesh, name=annot, figure=self._f) - def _prepare_data(self, data): - """Ensure data is float64 and has proper endianness. + # Set the color table + surf.module_manager.scalar_lut_manager.lut.table = cmap - Note: this is largely aimed at working around a Mayavi bug. + # Set the brain attributes + annot = dict(surface=surf, name=annot, colormap=cmap) + return annot + + def add_label(self, label, label_name, color, alpha): + """Add an ROI label to the image""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, + self._geo.y, + self._geo.z, + self._geo.faces, + scalars=label, + figure=self._f) + surf = mlab.pipeline.surface(mesh, name=label_name, figure=self._f) + color = colorConverter.to_rgba(color, alpha) + cmap = np.array([(0, 0, 0, 0,), color]) * 255 + surf.module_manager.scalar_lut_manager.lut.table = cmap + return surf + + def add_morphometry(self, morph_data, colormap, measure, min, max): + """Add a morphometry overlay to the image""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, + self._geo.y, + self._geo.z, + self._geo.faces, + scalars=morph_data, + figure=self._f) + + surf = mlab.pipeline.surface(mesh, colormap=colormap, + vmin=min, vmax=max, + name=measure, figure=self._f) + + # Get the colorbar + bar = mlab.scalarbar(surf) + self._format_cbar_text(bar) + bar.scalar_bar_representation.position2 = .8, 0.09 + + # Fil in the morphometry dict + return dict(surface=surf, colorbar=bar, measure=measure) + + def add_foci(self, foci_coords, scale_factor, color, alpha, name): + """Add spherical foci, possibly mapping to displayed surf""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + # Create the visualization + points = mlab.points3d(foci_coords[:, 0], + foci_coords[:, 1], + foci_coords[:, 2], + np.ones(foci_coords.shape[0]), + scale_factor=(10. * scale_factor), + color=color, opacity=alpha, name=name, + figure=self._f) + return points + + def add_contour_overlay(self, scalar_data, min=None, max=None, + n_contours=7, line_width=1.5): + """Add a topographic contour overlay of the positive data""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + # Set up the pipeline + mesh = mlab.pipeline.triangular_mesh_source(self._geo.x, self._geo.y, + self._geo.z, + self._geo.faces, + scalars=scalar_data, + figure=self._f) + thresh = mlab.pipeline.threshold(mesh, low=min) + surf = mlab.pipeline.contour_surface(thresh, contours=n_contours, + line_width=line_width) + + # Set the colorbar and range correctly + bar = mlab.scalarbar(surf, + nb_colors=n_contours, + nb_labels=n_contours + 1) + bar.data_range = min, max + self._format_cbar_text(bar) + bar.scalar_bar_representation.position2 = .8, 0.09 + + # Set up a dict attribute with pointers at important things + return dict(surface=surf, colorbar=bar) + + def add_text(self, x, y, text, name, color=None, opacity=1.0): + """ Add a text to the visualization""" + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + + return mlab.text(x, y, text, name=name, color=color, + opacity=opacity, figure=self._f) + + def _orient_lights(self): + """Set lights to come from same direction relative to brain.""" + if self.hemi == "rh": + if self._f.scene is not None and \ + self._f.scene.light_manager is not None: + for light in self._f.scene.light_manager.lights: + light.azimuth *= -1 + + def _get_geo_colors(self, config_opts): + """Return an mlab colormap name, vmin, and vmax for binary curvature. + + Parameters + ---------- + config_opts : dict + dictionary of config file "visual" options + + Returns + ------- + colormap : string + mlab colormap name + vmin : float + curv colormap minimum + vmax : float + curv colormap maximum + reverse : boolean + boolean indicating whether the colormap should be reversed """ - data = data.copy() - data = data.astype(np.float64) - if data.dtype.byteorder == '>': - data.byteswap(True) - return data + colormap_map = dict(classic=("Greys", -1, 2, False), + high_contrast=("Greys", -.1, 1.3, False), + low_contrast=("Greys", -5, 5, False), + bone=("bone", -.2, 2, True)) - def _format_cbar_text(self, cbar): + try: + cortex_color = config_opts['cortex'] + except KeyError: + cortex_color = config.get("visual", "cortex") + try: + color_data = colormap_map[cortex_color] + except KeyError: + warn("" + "The 'cortex' setting in your config file must be one of " + "'classic', 'high_contrast', 'low_contrast', or 'bone', " + "but your value is '%s'. I'm setting the cortex colormap " + "to the 'classic' setting." % cortex_color) + color_data = colormap_map['classic'] + + return color_data - bg_color = self.scene_properties["bgcolor"] + def _format_cbar_text(self, cbar): + bg_color = self._bg_color if bg_color is None or sum(bg_color) < 2: text_color = (1., 1., 1.) else: @@ -1735,114 +2456,115 @@ def _format_cbar_text(self, cbar): cbar.label_text_property.color = text_color -class Overlay(object): - """Encapsulation of statistical neuroimaging overlay visualization.""" +class OverlayData(object): + """Encapsulation of statistical neuroimaging overlay viz data""" def __init__(self, scalar_data, geo, min, max, sign): - try: - from mayavi import mlab - except ImportError: - from enthought.mayavi import mlab - if scalar_data.min() >= 0: sign = "pos" elif scalar_data.max() <= 0: sign = "neg" - self.sign = sign - - # Byte swap copy; due to mayavi bug - mlab_data = scalar_data.copy() - mlab_data = mlab_data.astype(np.float64) - if scalar_data.dtype.byteorder == '>': - mlab_data.byteswap(True) + self.geo = geo if sign in ["abs", "pos"]: - pos_mesh = mlab.pipeline.triangular_mesh_source(geo.x, - geo.y, - geo.z, - geo.faces, - scalars=mlab_data) - # Figure out the correct threshold to avoid TraitErrors # This seems like not the cleanest way to do this - pos_data = scalar_data[np.where(scalar_data > 0)] - try: - pos_max = pos_data.max() - except ValueError: - pos_max = 0 + pos_max = np.max((0.0, np.max(scalar_data))) if pos_max < min: thresh_low = pos_max else: thresh_low = min - pos_thresh = mlab.pipeline.threshold(pos_mesh, low=thresh_low) - pos_surf = mlab.pipeline.surface(pos_thresh, colormap="YlOrRd", - vmin=min, vmax=max) - pos_bar = mlab.scalarbar(pos_surf, nb_labels=5) - pos_bar.reverse_lut = True - - self.pos = pos_surf - self.pos_bar = pos_bar + self.pos_lims = [thresh_low, min, max] + else: + self.pos_lims = None if sign in ["abs", "neg"]: - neg_mesh = mlab.pipeline.triangular_mesh_source(geo.x, - geo.y, - geo.z, - geo.faces, - scalars=mlab_data) - # Figure out the correct threshold to avoid TraitErrors # This seems even less clean due to negative convolutedness - neg_data = scalar_data[np.where(scalar_data < 0)] - try: - neg_min = neg_data.min() - except ValueError: - neg_min = 0 + neg_min = np.min((0.0, np.min(scalar_data))) if neg_min > -min: thresh_up = neg_min else: thresh_up = -min - neg_thresh = mlab.pipeline.threshold(neg_mesh, up=thresh_up) - neg_surf = mlab.pipeline.surface(neg_thresh, colormap="PuBu", - vmin=-max, vmax=-min) - neg_bar = mlab.scalarbar(neg_surf, nb_labels=5) + self.neg_lims = [thresh_up, -max, -min] + else: + self.neg_lims = None + # Byte swap copy; due to mayavi bug + self.mlab_data = _prepare_data(scalar_data) - self.neg = neg_surf - self.neg_bar = neg_bar +class OverlayDisplay(): + """Encapsulation of overlay viz plotting""" + + def __init__(self, ol, figure): + try: + from mayavi import mlab + assert mlab + except ImportError: + from enthought.mayavi import mlab + args = [ol.geo.x, ol.geo.y, ol.geo.z, ol.geo.faces] + kwargs = dict(scalars=ol.mlab_data, figure=figure) + if ol.pos_lims is not None: + pos_mesh = mlab.pipeline.triangular_mesh_source(*args, **kwargs) + pos_thresh = mlab.pipeline.threshold(pos_mesh, low=ol.pos_lims[0]) + self.pos = mlab.pipeline.surface(pos_thresh, colormap="YlOrRd", + vmin=ol.pos_lims[1], + vmax=ol.pos_lims[2], + figure=figure) + self.pos_bar = mlab.scalarbar(self.pos, nb_labels=5) + self.pos_bar.reverse_lut = True + else: + self.pos = None + + if ol.neg_lims is not None: + neg_mesh = mlab.pipeline.triangular_mesh_source(*args, **kwargs) + neg_thresh = mlab.pipeline.threshold(neg_mesh, + up=ol.neg_lims[0]) + self.neg = mlab.pipeline.surface(neg_thresh, colormap="PuBu", + vmin=ol.neg_lims[1], + vmax=ol.neg_lims[2], + figure=figure) + self.neg_bar = mlab.scalarbar(self.neg, nb_labels=5) + else: + self.neg = None self._format_colorbar() def remove(self): - - if self.sign in ["pos", "abs"]: + if self.pos is not None: self.pos.remove() self.pos_bar.visible = False - if self.sign in ["neg", "abs"]: + if self.neg is not None: self.neg.remove() self.neg_bar.visible = False def _format_colorbar(self): - - if self.sign in ["abs", "neg"]: - self.neg_bar.scalar_bar_representation.position = (0.05, 0.01) - self.neg_bar.scalar_bar_representation.position2 = (0.42, 0.09) - if self.sign in ["abs", "pos"]: + if self.pos is not None: self.pos_bar.scalar_bar_representation.position = (0.53, 0.01) self.pos_bar.scalar_bar_representation.position2 = (0.42, 0.09) + if self.neg is not None: + self.neg_bar.scalar_bar_representation.position = (0.05, 0.01) + self.neg_bar.scalar_bar_representation.position2 = (0.42, 0.09) class TimeViewer(HasTraits): - """ TimeViewer object providing a GUI for visualizing time series, such - as M/EEG inverse solutions, on Brain object(s) + """TimeViewer object providing a GUI for visualizing time series + + Useful for visualizing M/EEG inverse solutions on Brain object(s). + + Parameters + ---------- + brain : Brain (or list of Brain) + brain(s) to control """ - # Nested import of traisui for setup.py without X server + # Nested import of traisui for setup.py without X server try: - from traitsui.api import View, Item, VSplit, HSplit, Group + from traitsui.api import (View, Item, VSplit, HSplit, Group) except ImportError: try: - from traits.ui.api import View, Item, VSplit, HSplit, Group + from traits.ui.api import (View, Item, VSplit, HSplit, Group) except ImportError: - from enthought.traits.ui.api import View, Item, VSplit, HSplit, Group - + from enthought.traits.ui.api import (View, Item, VSplit, + HSplit, Group) min_time = Int(0) max_time = Int(1E9) current_time = Range(low="min_time", high="max_time", value=0) @@ -1862,24 +2584,16 @@ class TimeViewer(HasTraits): Group(HSplit(Item(name="fmin"), Item(name="fmid"), Item(name="fmax"), - Item(name="transparent"), - ), + Item(name="transparent") + ), label="Color scale", - show_border=True - ), - Item(name="smoothing_steps"), - Item(name="orientation") - ) + show_border=True), + Item(name="smoothing_steps"), + Item(name="orientation") + ) ) def __init__(self, brain): - """Initialize TimeViewer - - Parameters - ---------- - brain : Brain (or list of Brain) - brain(s) to control - """ super(TimeViewer, self).__init__() if isinstance(brain, (list, tuple)):