Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Jump to bottom

[ENH] Neuropythy *_to_dva #606

Merged
merged 2 commits into from
Jan 29, 2024
Merged

[ENH] Neuropythy *_to_dva #606

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
237782c
add skeleton
jgranley Jan 27, 2024
ba8ec9e
add cortex_to_dva
jgranley Jan 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
120 changes: 117 additions & 3 deletions pulse2percept/topography/neuropythy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Open in desktop
Original file line number Diff line number Diff line change
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)