[DOC] Neuropythy example (#602)

* add neuropythy example

* add neuropythy example

* update doc

* restructure cortical.rst

doc/topics/cortical.rst

@@ -18,19 +18,17 @@ left and right.  The visual cortex is divided into multiple regions, including
 v1, v2, and v3, with each region performing a different function required
 to process visual information.
 
-Each region processes an aspect (such as color or motion) of the entire visual
-field.  Within a region, different parts of the visual field are processed by
+Within a region, different parts of the visual field are processed by
 different neurons.  We can define a mapping between locations in the visual field
-and locations in the cortex.  This mapping is called a visual field map, or
-topography.
+and locations in the cortex.  This mapping is called a visual field map, also
+called retinotopy or visuotopy.
 
 Model Plotting
 ^^^^^^^^^^^^^^
-One way to visualize the mapping between the visual field and the cortex is
-to plot a spatial model.  A spatial model consists of a set of points in the
+One way to visualize the mapping between the visual field and the cortex in pulse2percept
+is to plot a model. A model simulates a set of points in the
 visual field and the corresponding points in the cortex (using a visual field
-map).  The plot of a model shows all of these points, either in the visual
-field or on the cortex depending on the parameters used to create the plot.
+map).
 
 The first step is to create a model, for example
 :py:class:`~pulse2percept.models.cortex.ScoreboardModel`.  We can create the
@@ -59,21 +57,17 @@ spaced, and are represented by `+` symbols.
 
 If we don't set `use_dva=True`, then the visual field mapping will be applied
 to the points in the visual field, and the points on the cortex will be
-plotted instead.  Because we created a model with three regions, each point
-in the visual field will be transformed to three new points: one in v1, one
-in v2, and one in v3.
+plotted instead. ScoreboardModel by default uses the 
+:py:class:`~pulse2percept.topography.Polimeni2006Map` visual field map, but 
+it can also use :py:class:`~pulse2percept.topography.NeuropythyMap` for
+3D patient-specific MRI based retinotopy.
 
 The cortex is split into left and right hemispheres, with each side being
 responsible for processing information from one eye.  In reality, the left
-and right hemispheres of our brain are disconnected, but to simplify the
-code, pulse2percept represents them as one continuous space.  The origin
-of both hemispheres corresponds to the center of our visual field, as defined
-by the visual field map.  However, because we represent both hemispheres as
-one continuous space, this would result in both hemispheres overlapping
-around the origin.
-
-To avoid this, the left hemisphere are offset by 20mm, meaning the origin
-of the left hemisphere is (-20, 0).  In addition, cortical visual field maps
+and right hemispheres of our brain are disconnected, but to simplify 
+the Polimeni map represents them as one continuous space. 
+The left hemisphere is offset by 20mm, meaning the origin
+of the left hemisphere is (-20000, 0).  In addition, cortical visual field maps
 have a `split_map` attribute set to `True`, which means that no current will
 be allowed to cross between the hemispheres.
 
@@ -89,7 +83,7 @@ center of our visual field is represented by a larger area on the cortex than
 equally sized area at the periphery of our visual field, an effect called
 cortical magnification.
 
-Another style option for the plot is `"hull"`` (the default):
+Another style option for the plot is `"hull"`:
 
 .. ipython:: python
     :okwarning:
@@ -186,13 +180,53 @@ and "counter" electrodes.
 
 .. _topics-ensemble-implant:
 
+Neuralink
+^^^^^^^^^
+:py:class:`~pulse2percept.implants.cortex.Neuralink` is an implant 
+consisting of multiple Neuralink threads. Currently the only thread implemented
+is the :py:class:`~pulse2percept.implants.cortex.LinearEdgeThread` which 
+consists of 32 electrodes. 
+
+.. ipython:: python
+
+    from pulse2percept.implants.cortex import LinearEdgeThread
+    thread = LinearEdgeThread()
+    thread.plot3D()
+    @savefig neuralink_thread.png align=center
+    plt.axis('equal')
+
+
+Neuralink works well with the :py:class:`~pulse2percept.topography.NeuropythyMap`,
+which is a 3D patient-specific MRI based retinotopy. You can easily create
+a Neuralink implant with multiple threads using the NeuropythyMap as follows:
+
+.. ipython:: python
+    :okwarning:
+
+    from pulse2percept.implants.cortex import Neuralink
+    from pulse2percept.topography import NeuropythyMap
+    from pulse2percept.models.cortex import ScoreboardModel
+    map = NeuropythyMap('fsaverage', regions=['v1'])
+    model = ScoreboardModel(vfmap=map, xrange=(-4, 0), yrange=(-4, 4), xystep=.25).build()
+    neuralink = Neuralink.from_neuropythy(map, xrange=model.xrange, yrange=model.yrange, xystep=1, rand_insertion_angle=0)
+    fig = plt.figure(figsize=(10, 5))
+    ax1 = fig.add_subplot(121, projection='3d')
+    neuralink.plot3D(ax=ax1)
+    model.plot3D(style='cell', ax=ax1)
+    ax2 = fig.add_subplot(122)
+    neuralink.plot(ax=ax2)
+    model.plot(style='cell', ax=ax2)
+    @savefig neuralink.png align=center
+    plt.show()
+
+
 Ensemble Implants
 -----------------
 
 :py:class:`~pulse2percept.implants.EnsembleImplant` is a new class which
 allows the user to use multiple implants in tandem. It can be used with any 
 implant type, but was made for use with small implants meant to be used together,
-such as :py:class:`~pulse2percept.implants.cortex.Cortivis`. This tutorial will 
+such as :py:class:`~pulse2percept.implants.cortex.ICVP`. This tutorial will 
 demonstrate how to create an :py:class:`~pulse2percept.implants.EnsembleImplant`,
 to combine multiple :py:class:`~pulse2percept.implants.cortex.Cortivis` objects.
 
@@ -224,6 +258,7 @@ is the index of the implant in the constructor list. Implants can also be passed
 a dictionary, in which case the naming pattern is `key-electrode` where `key` is the
 electrode's dictionary key.
 
+
 .. _topics-cortical-models:
 
 Models
@@ -364,6 +399,27 @@ see that the predicted percept is now larger due to cortical magnification:
     percept.plot()
 
 
+Pulse2percept currently has 2 cortical models, :py:class:`~pulse2percept.models.cortex.ScoreboardModel` 
+and :py:class:`~pulse2percept.models.cortex.DynaphosModel`. The ScoreboardModel 
+is a simple model that assumes that each electrode creates a circular patch of 
+brightness. The DynaphosModel is a more complex model that takes into account
+both spatial current spread and temporal effects such as charge accumulation. 
+
+.. ipython:: python
+
+    from pulse2percept.models.cortex import DynaphosModel
+    from pulse2percept.stimuli import BiphasicPulseTrain
+    from pulse2percept.implants.cortex import Orion
+
+    model = DynaphosModel().build()
+    implant = Orion()
+    implant.stim = {e : BiphasicPulseTrain(20, 200, .45) for e in implant.electrode_names}
+    percept = model.predict_percept(implant)
+    @savefig model_dynaphos.png align=center
+    percept.plot()
+
+You can also play the percept as a video with `percept.play()`.
+
 .. _topics-cortical-developers:
 
 For Developers
@@ -373,16 +429,6 @@ In this section we will discuss some of the changes made under the hood
 accomadate cortical features, as well as some important notes for developers
 to keep in mind.
 
-Grid Plotting
-^^^^^^^^^^^^^
-Previously, the grid's plot function took in an optional `transform`
-function that would be applied to all of the points.  This parameter has been
-removed, and instead the plot function will automatically apply all of the
-transforms in the grid's `retinotopy` attribute.  This is a dictionary of
-different transforms that can be applied to the grid, such as dva to retinal
-coordinates or dva to cortical coordinates.  If you want to plot the grid
-without transformed points, you can pass in `use_dva=True`.
-
 Units
 ^^^^^
 Keep in mind that pulse2percept uses units of microns for length, microamps

doc/users/references.rst

@@ -19,6 +19,9 @@ Studies referenced throughout the Documentation:
 .. [Al-Atabany2010] WT Al-Atabany, T Tong, PA Degenaar (2010). Improved 
                     content aware scene retargeting for retinitis pigmentosa
                     patients. *BioMedical Engineering OnLine*, 9(1), 52. 
+.. [Benson2018] Benson NC, Winawer J (2018) Bayesian Analysis of Retinotopic Maps. 
+                eLife 2018, doi:`10.1101/325597
+                <https://www.biorxiv.org/content/10.1101/325597v4>`_.
 .. [Beyeler2019] M Beyeler, D Nanduri, JD Weiland, A Rokem, GM Boynton, I Fine
                  (2019). A model of ganglion axon pathways accounts for
                  percepts elicited by retinal implants. *Scientific Reports*

examples/implants/plot_implants_cortical.py

@@ -12,6 +12,9 @@
 
 import matplotlib.pyplot as plt
 from pulse2percept.implants.cortex import *
+import pulse2percept as p2p
+import warnings
+warnings.filterwarnings("ignore") # ignore matplotlib warnings
 
 orion = Orion()
 orion.plot(annotate=True)
@@ -38,4 +41,51 @@
 
 icvp = ICVP()
 icvp.plot(annotate=True)
-plt.show()
+plt.show()
+
+###############################################################################
+# Neuralink Threads (Neuralink)
+# --------------------------------------
+#
+# :py:class:`~pulse2percept.implants.cortex.Neuralink` is an implant consisting 
+# of one or more :py:class:`~pulse2percept.implants.cortex.NeuralinkThread`s. 
+# The simplest one is :py:class:`~pulse2percept.implants.cortex.LinearEdgeThread`,
+# consisting of 32 electrodes.
+#
+# These threads are designed to be 3D, and can be plotted in both 3D and 2D space,
+# which will show the x-y projection of the points. Here are some examples
+# showing the 1) a close up of the thread, 2) a 3D plot of two threads, and 3) the
+# 2D projection of the two threads.
+fig = plt.figure(figsize=(10, 4))
+thread1 = LinearEdgeThread()
+thread2 = LinearEdgeThread(500, 500, orient=[1, 1, 1])
+ax0 = fig.add_subplot(131, projection='3d')
+thread1.plot3D(ax=ax0)
+ax0.set_xlim(-100, 100)
+ax0.set_ylim(-100, 100)
+ax0.set_zlim(10, 200)
+ax1 = fig.add_subplot(132, projection='3d')
+thread1.plot3D(ax=ax1)
+thread2.plot3D(ax=ax1)
+ax2 = fig.add_subplot(133)
+thread1.plot(ax=ax2)
+thread2.plot(ax=ax2)
+plt.axis('equal')
+plt.show()
+
+###############################################################################
+# Neuralink implants can be easily created from a NeuropythyMap
+# enabling easy placement of implants across the cortical surface.
+# See the neuropythy example for more details.
+nmap = p2p.topography.NeuropythyMap(subject='fsaverage', regions=['v1'])
+model = p2p.models.cortex.ScoreboardModel(rho=500, xrange=(-6, 0), yrange=(-5, 5), xystep=.25, vfmap=nmap).build()
+nlink = Neuralink.from_neuropythy(nmap, xrange=model.xrange, yrange=model.yrange, 
+                                                      xystep=2, rand_insertion_angle=0)
+print(len(nlink.implants), " total threads")
+fig = plt.figure(figsize=(10, 4))
+ax1 = fig.add_subplot(121, projection='3d')
+model.plot3D(ax=ax1, style='cell')
+nlink.plot3D(ax=ax1)
+ax2 = fig.add_subplot(122)
+model.plot(style='cell', ax=ax2)
+nlink.plot(ax=ax2)

examples/models/plot_neuropythy.py

@@ -0,0 +1,139 @@
+# -*- coding: utf-8 -*-
+"""
+===============================================================================
+Neuropythy and Neuralink: Patient specific visual field maps based on MRI
+===============================================================================
+
+Neuropythy [Benson2018]_ is a python package that predicts patient-specific 
+visuotopies based on MRI scans of the human visual cortex. It is possible to use 
+Neuropythy within pulse2percept as the visual field map for cortical models.
+
+First, make sure neuropythy is installed, which can be done with
+`pip install neuropythy`.
+
+
+Creating a Neuropythy visual field map
+--------------------------------------
+:py:class:`~pulse2percept.topography.NeuropythyMap` requires a subject parameter,
+which can be either:
+1) 'fsaverage' (the average of freesurfer subjects)
+2) a string with a subject from the Benson Winawer 2018 dataset ('S1201'-'S1208')
+3) the path to a freesurfer subject directory in a format neuropythy can load
+4) the name of a subject in any directory specified by `ny.config['freesurfer_subject_paths']` or
+5) an instantiated `neuropythy.mri.core.Subject`
+
+If you do not have it downloaded already, the Benson Winawer 2018 dataset will
+be automatically downloaded and cached to the provided `cache_dir` parameter,
+which defaults to `~/.neuropythy_p2p`. If you have already downloaded the dataset,
+elsewhere, just add the subjects folder of the dataset to `ny.config['freesurfer_subject_paths']` and p2p
+will be able to load subjects from your predownloaded directory.
+
+"""
+# sphinx_gallery_thumbnail_number = 5
+
+import pulse2percept as p2p
+import matplotlib.pyplot as plt
+import warnings
+warnings.filterwarnings("ignore") # ignore matplotlib warnings
+
+# this could take a while if you haven't already downloaded the dataset
+nmap = p2p.topography.NeuropythyMap(subject='fsaverage', regions=['v1'])
+nmap
+
+##################################################################################
+# NeuropythyMap provides a number of methods to transform visual field coordinates
+# into cortical coordinates:
+#
+# * :py:meth:`~pulse2percept.topography.NeuropythyMap.dva_to_v1`
+# * :py:meth:`~pulse2percept.topography.NeuropythyMap.dva_to_v2`
+# * :py:meth:`~pulse2percept.topography.NeuropythyMap.dva_to_v3`
+#
+# In contrast to other visual field maps, you may also specify a cortical surface
+# that the visual field coordinates should be mapped to. By default, the cortical
+# surface is set to 'midgray', but you can also specify 'white' or 'pial', or other
+# surfaces accepted by neuropythy.
+#
+#
+# Lets use our map in a model and visualize the transform
+
+model = p2p.models.cortex.ScoreboardModel(vfmap=nmap, regions=['v1'])
+model.build()
+fig, axes = plt.subplots(ncols=2, figsize=(10, 4))
+model.plot(style='cell', ax=axes[0])
+model.plot(style='scatter', ax=axes[1])
+
+##################################################################################
+# Note, if we call the usual plot() method, it will be plotting a 2D
+# projection of the points. This can look quite strange sometimes, since the
+# points are inherently 3D. To plot the points in 3D, we can use the
+# plot3D() method.
+
+fig = plt.figure(figsize=(10, 4))
+ax1 = fig.add_subplot(121, projection='3d')
+model.plot3D(ax=ax1, style='cell')
+ax2 = fig.add_subplot(122, projection='3d')
+model.plot3D(ax=ax2, style='scatter')
+
+##################################################################################
+# Neuropythy can use up to 'v1', 'v2', 'v3'. If you want to use all three regions,
+# you can specify them as a list. 
+#
+# Lets use all three regions and plot the result (note it can get a little messy):
+nmap = p2p.topography.NeuropythyMap(subject='fsaverage', regions=['v1', 'v2', 'v3'])
+model = p2p.models.cortex.ScoreboardModel(xrange=(-20, 20), yrange=(-20, 20), xystep=0.5,
+                                          vfmap=nmap, regions=['v1', 'v2', 'v3'])
+fig = plt.figure(figsize=(10, 4))
+ax1 = fig.add_subplot(121, projection='3d')
+model.plot3D(ax=ax1, style='cell')
+ax2 = fig.add_subplot(122, projection='3d')
+model.plot3D(ax=ax2, style='scatter')
+
+
+##################################################################################
+# 
+# Note that the regions are not necessarily
+# contiguous, and the map will likely be discontinuous at the boundaries between
+# visual areas. To combat this, you can set the `jitter_boundary` to True
+# to add a small amount of noise to the boundary to make it less noticeable.
+# This is turned off by default.
+#
+# Placing implants with neuropythy maps
+# -------------------------------------
+# When using a cortical map, it is important to place your implant in the correct
+# 3D location; the default z value of 0 will likely not be on a cortical surface.
+# 
+# 
+# For the Neuralink implant, there exists a helper function, 
+# :py:meth:`~pulse2percept.implants.NeuralinkImplant.from_neuropythy`, which will
+# automatically place neuralink threads at specified visual field locations across 
+# the cortical surface. Threads will be inserted perpendicular to the cortical surface
+# up to `rand_insertion_angle`. Similar helper functions are under development
+# for other cortical implants.
+#
+# Lets place a Neuralink implant across the right hemisphere of the cortex:
+nmap = p2p.topography.NeuropythyMap(subject='fsaverage', regions=['v1'])
+model = p2p.models.cortex.ScoreboardModel(vfmap=nmap, xrange=(-4, 0), yrange=(-4, 4), xystep=.25).build()
+nlink = p2p.implants.cortex.Neuralink.from_neuropythy(nmap, xrange=model.xrange, yrange=model.yrange, 
+                                                      xystep=1,rand_insertion_angle=0)
+print(len(nlink.implants), " total threads")
+fig = plt.figure(figsize=(10, 4))
+ax1 = fig.add_subplot(121, projection='3d')
+model.plot3D(ax=ax1, style='cell')
+nlink.plot3D(ax=ax1)
+ax2 = fig.add_subplot(122)
+model.plot(style='cell', ax=ax2)
+nlink.plot(ax=ax2)
+
+##################################################################################
+# Finally, lets predict what a percept would look like if we stimulated one 
+# electrode on each thread, using the scoreboard model:
+nmap = p2p.topography.NeuropythyMap(subject='fsaverage', regions=['v1'], jitter_boundary=True)
+model = p2p.models.cortex.ScoreboardModel(rho=800, xrange=(-15, 15), yrange=(-15, 15), xystep=.25, vfmap=nmap).build()
+nlink = p2p.implants.cortex.Neuralink.from_neuropythy(nmap, xrange=model.xrange, yrange=model.yrange, 
+                                                      xystep=3,rand_insertion_angle=0)
+nlink.stim = {e : 1 for e in [i for i in nlink.electrode_names if i[-1] == '1']}
+percept = model.predict_percept(nlink)
+percept.plot()
+
+
+

examples/models/plot_visual_field_maps.py

@@ -27,6 +27,10 @@
    provides a mapping between visual field coordinates and cortical coordinates 
    using a wedge-dipole model for regions V1, V2, and V3. See Appendix B of 
    [Polimeni2006]_ for details.
+* :py:class:`~pulse2percept.topography.NeuropythyMap`: Neuropythy is a python
+   package that predicts patient-specific visuotopies based on MRI scans of the
+   human visual cortex. It provides a mapping between visual field coordinates
+   and 3D cortical coordinates for regions V1, V2, V3. See [Benson2018]_ for details.
 
 All of these visual field maps follow the templates in either 
 :py:class:`~pulse2percept.topography.RetinalMap` or