Skip to content

Commit

Permalink
[ENH] Neuropythy *_to_dva (#606)
Browse files Browse the repository at this point in the history 
* add skeleton

* add cortex_to_dva
  • Loading branch information
@jgranley
jgranley committed Jan 29, 2024
1 parent 762244f commit 6e8471a
Show file tree
Hide file tree
Showing 2 changed files with 224 additions and 3 deletions.
120 changes: 117 additions & 3 deletions pulse2percept/topography/neuropythy.py
View file
@@ -1,5 +1,6 @@
import numpy as np
import os
from scipy.spatial import cKDTree

from .cortex import CorticalMap

Expand Down Expand Up @@ -56,6 +57,8 @@ def get_default_params(self):
'jitter_thresh' : 0.3,
# no split map
'left_offset' : None,
# max nearest neighbor distance for v*_to_dva
'cort_nn_thresh' : 1000, # 1mm
}
return {**super().get_default_params(), **params}

Expand Down Expand Up @@ -100,6 +103,10 @@ def load_meshes(self, subject):
left, right = ny.vision.predict_retinotopy(subject, sym_angle=False)

self.predicted_retinotopy = (left, right)
cortex_pts = [] # (npoints, 3)
addr_idxs = [] # (npoints, 2)
region_idxs = [] # (npoints, 2)
hemi_idxs = [] # (npoints, 2)
region_meshes = {}
for region in self.regions:
region_lbl = int(region[-1])
Expand All @@ -109,13 +116,25 @@ def load_meshes(self, subject):
(ang, ecc) = (retinotopy['angle'], retinotopy['eccen'])
ang = np.pi/180 * (90 - ang)
(x, y) = (ecc*np.cos(ang), ecc*np.sin(ang))
# doesn't matter what surface is used here for dva_to_*
# but use pial so plotting is easier
submesh = subject.hemis[hemi].pial_surface.submesh(ii)
# doesn't matter what surface is used here for dva_to_v*
# but grab all the points for a ckdtree for v*_to_dva
for surface in ['midgray', 'white', 'pial']:
submesh = subject.hemis[hemi].surface(surface).submesh(ii)
cortex_pts.append(submesh.coordinates.T.astype('float32'))
addr_idxs.append(np.arange(submesh.coordinates.shape[1]))
region_idxs.append([region for _ in range(submesh.coordinates.shape[1])])
hemi_idxs.append([1 if hemi=='rh' else 0 for _ in range(submesh.coordinates.shape[1])])
ii = submesh.labels
submesh = submesh.copy(coordinates=[x[ii], y[ii]])
vfmeshes.append(submesh)
region_meshes[region] = tuple(vfmeshes)

self.addr_idxs = {'addr' : np.concatenate(addr_idxs),
'region' : np.concatenate(region_idxs),
'hemi' : np.concatenate(hemi_idxs)}
cortex_pts = np.concatenate(cortex_pts)
self.cortex_tree = cKDTree(cortex_pts)

return region_meshes


Expand Down Expand Up @@ -278,6 +297,101 @@ def dva_to_v3(self, x, y, surface='midgray'):
ret[:, ~idx] = l
ret = ret.reshape((3, *shape))
return ret[0], ret[1], ret[2]


def cortex_to_dva(self, xc, yc, zc):
"""
Gives the visual field position(s) of the cortex point(s) `(xc,yc,zc)`.
Parameters
----------
xc, yc, zc : float or array_like
The x, y, and z-coordinate(s) of the cortex point(s) to look up (in mm).
Returns
-------
x, y : array_like
The x and y-coordinate(s) of the visual field point(s) (in dva).
"""
xc = np.array(xc, dtype='float32')
yc = np.array(yc, dtype='float32')
zc = np.array(zc, dtype='float32')
if np.shape(xc) != np.shape(yc) or np.shape(xc) != np.shape(zc):
raise ValueError("x, y, and z must have the same shape")
id_nan = np.isnan(xc) | np.isnan(yc) | np.isnan(zc)
query = np.stack([np.ravel(xc[~id_nan]), np.ravel(yc[~id_nan]), np.ravel(zc[~id_nan])], axis=-1) / 1000 # convert to mm
if np.size(query) == 0:
return np.ones((*np.shape(xc), 2)) * np.nan
dist, idx = self.cortex_tree.query(query, k=5)#, distance_upper_bound=self.cort_nn_thresh / 1000)
idx_nan = np.all(dist > self.cort_nn_thresh / 1000, axis=-1)
dist[dist > self.cort_nn_thresh / 1000] = 999999999 # make high so it doesn't contribute to geometric mean
neighbors = np.array([[self.region_meshes[self.addr_idxs['region'][i]][self.addr_idxs['hemi'][i]].coordinates[:, self.addr_idxs['addr'][i]]
for i in nb_pts]
for nb_pts in idx])
# use geometric mean based on distance
pts = np.sum(neighbors * 1/dist[..., None], axis=1) / np.sum(1/dist, axis=1)[:, None]
pts[idx_nan] = [np.nan, np.nan]
out = np.ones((*np.shape(xc), 2)) * np.nan
out[~id_nan] = pts
return pts[:, 0], pts[:, 1]


def v1_to_dva(self, xv1, yv1, zv1):
"""
Convert points in v1 to dva. Uses the mean of the 5 nearest neighbors
in the cortical mesh, weighted by 1/distance, to interpolate dva.
Any points that are more than self.cort_nn_thresh um from the
nearest neighbor will be set to nan.
Parameters
----------
xv1, yv1, zv1 : float or array_like
The x, y, and z-coordinate(s) of the v1 point(s) to look up (in mm).
Returns
-------
x, y : array_like
The x and y-coordinate(s) of the visual field point(s) (in dva).
"""
return self.cortex_to_dva(xv1, yv1, zv1)

def v2_to_dva(self, xv2, yv2, zv2):
"""
Convert points in v2 to dva. Uses the mean of the 5 nearest neighbors
in the cortical mesh, weighted by 1/distance, to interpolate dva.
Any points that are more than self.cort_nn_thresh um from the
nearest neighbor will be set to nan.
Parameters
----------
xv2, yv2, zv2 : float or array_like
The x, y, and z-coordinate(s) of the v2 point(s) to look up (in mm).
Returns
-------
x, y : array_like
The x and y-coordinate(s) of the visual field point(s) (in dva).
"""
return self.cortex_to_dva(xv2, yv2, zv2)

def v3_to_dva(self, xv3, yv3, zv3):
"""
Convert points in v3 to dva. Uses the mean of the 5 nearest neighbors
in the cortical mesh, weighted by 1/distance, to interpolate dva.
Any points that are more than self.cort_nn_thresh um from the
nearest neighbor will be set to nan.
Parameters
----------
xv3, yv3, zv3 : float or array_like
The x, y, and z-coordinate(s) of the v3 point(s) to look up (in mm).
Returns
-------
x, y : array_like
The x and y-coordinate(s) of the visual field point(s) (in dva).
"""
return self.cortex_to_dva(xv3, yv3, zv3)



107 changes: 107 additions & 0 deletions pulse2percept/topography/tests/test_neuropythy.py
View file
Expand Up @@ -221,5 +221,112 @@ def test_neuropythy_scoreboard():



@pytest.mark.slow()
@pytest.mark.parametrize('regions', [['v1'], ['v1', 'v3'], ['v1', 'v2', 'v3']])
def test_cortex_to_dva(regions):
nmap = NeuropythyMap('fsaverage', regions=regions, jitter_boundary=True)
npt.assert_equal(nmap.predicted_retinotopy is not None, True)
npt.assert_equal(nmap.region_meshes is not None, True)
if 'v1' in regions:
npt.assert_equal(nmap.region_meshes['v1'] is not None, True)
if 'v2' in regions:
npt.assert_equal(nmap.region_meshes['v2'] is not None, True)
if 'v3' in regions:
npt.assert_equal(nmap.region_meshes['v3'] is not None, True)


npt.assert_equal(list(nmap.region_meshes.keys()), regions)

if 'v1' in regions:
# should work with all shapes
npt.assert_equal(nmap.v1_to_dva(0, 0, 0)[0], np.array([np.nan]))
nmap.v1_to_dva([100, 200, 300], [100, 200, 300], [100, 200, 300])
nmap.v1_to_dva(np.eye(3), np.eye(3), np.eye(3))

x = np.array([-10035.355, -13315.073, 12075.739, 13630.971])
y = np.array([ -96637.12, -102852.29, -95358.4, -101546.41])
z = np.array([-10769.129, -3861.491, -7168.826, 924.938])


xdva, ydva = nmap.v1_to_dva(x, y, z)
npt.assert_equal(x.shape, (4,))
npt.assert_equal(y.shape, (4,))
npt.assert_almost_equal(xdva, np.array([1, 1, -1, -1]), decimal=1)
npt.assert_almost_equal(ydva, np.array([1, -1, 1, -1]), decimal=1)

x = np.arange(-10, -1, .1)
y = np.arange(-10, -1, .1)
x1, y2 = nmap.v1_to_dva(*nmap.dva_to_v1(x, y))
npt.assert_allclose(x, x1, rtol=.05, atol=0.1)
npt.assert_allclose(y, y2, rtol=.05, atol=0.1)


# test cort_nn_thresh
idx = np.argmax(nmap.subject.hemis['rh'].surface('midgray').coordinates[0])
x = np.array([nmap.subject.hemis['rh'].surface('midgray').coordinates[0][idx]])
y = np.array([nmap.subject.hemis['rh'].surface('midgray').coordinates[1][idx]])
z = np.array([nmap.subject.hemis['rh'].surface('midgray').coordinates[2][idx]])
xdva, ydva = nmap.v1_to_dva(x, y, z)
npt.assert_equal(xdva != np.array([np.nan]), True)
npt.assert_equal(ydva != np.array([np.nan]), True)
x1 = x + 999
xdva, ydva = nmap.v1_to_dva(x1, y, z)
npt.assert_equal(xdva != np.array([np.nan]), True)
npt.assert_equal(ydva != np.array([np.nan]), True)
x1 = x +1001
xdva, ydva = nmap.v1_to_dva(x1, y, z)
npt.assert_equal(xdva, np.array([np.nan]))
npt.assert_equal(ydva, np.array([np.nan]))




if 'v2' in regions:
npt.assert_equal(nmap.v2_to_dva(0, 0, 0)[0], np.array([np.nan]))
nmap.v2_to_dva([100, 200, 300], [100, 200, 300], [100, 200, 300])
nmap.v2_to_dva(np.eye(3), np.eye(3), np.eye(3))


x = np.array([-11731.504, -20458.03, 22283.799] )
y = np.array([ -93461.92, -100803.35, -99334.945])
z = np.array([-11246.644, 1673.845, 7011.859])

xdva, ydva = nmap.v2_to_dva(x, y, z)
npt.assert_equal(xdva.shape, (3,))
npt.assert_equal(ydva.shape, (3,))
npt.assert_allclose(xdva, np.array([1, 1, -1]), rtol=.05, atol=0.1)
npt.assert_allclose(ydva, np.array([1, -1, -1]), rtol=.05, atol=0.1)

x = np.arange(-10, -1, .1)
y = np.arange(-10, -1, .1)
x1, y2 = nmap.v2_to_dva(*nmap.dva_to_v2(x, y))
npt.assert_allclose(x, x1, rtol=.05, atol=0.1)
npt.assert_allclose(y, y2, rtol=.05, atol=0.1)


if 'v3' in regions:
npt.assert_equal(nmap.v3_to_dva(0, 0, 0)[0], np.array([np.nan]))
nmap.v3_to_dva([100, 200, 300], [100, 200, 300], [100, 200, 300])
nmap.v3_to_dva(np.eye(3), np.eye(3), np.eye(3))


x = np.array([-23812.113, -23514.828, 28547.275])
y = np.array([-84409.51, -93015.07, -93238.63])
z = np.array([-15261.302, 4050.124, 8467.487])

xdva, ydva = nmap.v3_to_dva(x, y, z)

npt.assert_equal(xdva.shape, (3,))
npt.assert_equal(ydva.shape, (3,))
npt.assert_allclose(xdva, np.array([1, 1, -1]), rtol=.05, atol=0.1)
npt.assert_allclose(ydva, np.array([1, -1, -1]), rtol=.05, atol=0.1)

x = np.arange(-10, -1, .1)
y = np.arange(-10, -1, .1)
x1, y2 = nmap.v3_to_dva(*nmap.dva_to_v3(x, y))
npt.assert_allclose(x, x1, rtol=.05, atol=0.1)
npt.assert_allclose(y, y2, rtol=.05, atol=0.1)




0 comments on commit 6e8471a

Please sign in to comment.