From 6e8471a786851204fe6bfa38778d3c575612d602 Mon Sep 17 00:00:00 2001 From: Jacob Granley Date: Mon, 29 Jan 2024 00:11:57 -0800 Subject: [PATCH] [ENH] Neuropythy *_to_dva (#606) * add skeleton * add cortex_to_dva --- pulse2percept/topography/neuropythy.py | 120 +++++++++++++++++- .../topography/tests/test_neuropythy.py | 107 ++++++++++++++++ 2 files changed, 224 insertions(+), 3 deletions(-) diff --git a/pulse2percept/topography/neuropythy.py b/pulse2percept/topography/neuropythy.py index dfb57983..6e255003 100644 --- a/pulse2percept/topography/neuropythy.py +++ b/pulse2percept/topography/neuropythy.py @@ -1,5 +1,6 @@ import numpy as np import os +from scipy.spatial import cKDTree from .cortex import CorticalMap @@ -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} @@ -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]) @@ -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 @@ -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) diff --git a/pulse2percept/topography/tests/test_neuropythy.py b/pulse2percept/topography/tests/test_neuropythy.py index 574011fe..bb3e2871 100644 --- a/pulse2percept/topography/tests/test_neuropythy.py +++ b/pulse2percept/topography/tests/test_neuropythy.py @@ -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) + + \ No newline at end of file