From 335dfa51ce9a4d8d342089e86b70fc98f0ec5b21 Mon Sep 17 00:00:00 2001 From: Matthew Shinkle Date: Thu, 3 Aug 2023 08:44:53 -0700 Subject: [PATCH 01/15] feat: add 'dartboard' masking and plotting --- cortex/dartboards.py | 1003 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1003 insertions(+) create mode 100755 cortex/dartboards.py diff --git a/cortex/dartboards.py b/cortex/dartboards.py new file mode 100755 index 00000000..f663ed36 --- /dev/null +++ b/cortex/dartboards.py @@ -0,0 +1,1003 @@ +# pylint: disable=missing-module-docstring +from matplotlib import pyplot as plt +from matplotlib import transforms as mtrans +from matplotlib.colors import Normalize +from sklearn.metrics import pairwise_distances +import matplotlib +import numpy as np +import os +from . import db +from . import options +from . import polyutils +from . import quickflat + +from .dataset.views import Vertex + +from scipy.spatial import ConvexHull # pylint: disable-msg=E0611 +from .utils import get_cmap + +# Utils + + +def flatten_list_of_lists(list_of_lists): + flat_list = [item for sublist in list_of_lists for item in sublist] + return flat_list + + +def sph2cart(r, az, elev): + """Convert spherical to cartesian coordinates + + Parameters + ---------- + r : scalar or array-like + radius + az : scalar or array-like + azimuth angle in degrees + elev : scalar or array-like + elevation angle in degrees + """ + z = r * np.sin(np.radians(elev)) + rcoselev = r * np.cos(np.radians(elev)) + x = rcoselev * np.cos(np.radians(az)) + y = rcoselev * np.sin(np.radians(az)) + return x, y, z + + +def cart2pol(x, y): + phi = np.arctan2(y, x) + rho = np.sqrt(x**2 + y**2) + return phi, rho + +def pol2cart(phi, rho): + x = rho * np.cos(phi) + y = rho * np.sin(phi) + return x, y + + +def nonzero_coords(coords): + return coords[np.any(coords != 0, axis=1)] + +def outline_to_polar(outline, centroid, theta_direction='ccw'): + phi, rho = cart2pol(*(outline.T-centroid[:, np.newaxis])) + if theta_direction not in (-1, 0, 'counter-clockwise', 'ccw'): + phi = np.pi - phi + return np.stack([phi, rho], axis=0) + +def interpolate_angles(angles, angle_sub_bins, period=2*np.pi): + def single_interpolate(angles): + angles = list(angles) + if angles[-1] % period != angles[0] % period: + angles.append(angles[0]) + for i in range(1, len(angles)): + if angles[i] < angles[i-1]: + angles[i] += period + interpolated_angles = [] + for i in range(len(angles)-1): + interpolated_angles += [a for a in np.linspace( + angles[i], angles[i+1], angle_sub_bins, endpoint=False)] + interpolated_angles.append(angles[-1]) + return np.array(interpolated_angles) % period + if np.array(angles).ndim == 2: + return np.array([single_interpolate(hemi_angles) for hemi_angles in angles]) + else: + return single_interpolate(angles) + +def warp_phis(phis, current_anchors, new_anchors=None, period=2*np.pi): + if isinstance(new_anchors, type(None)): + new_anchors = np.linspace(0,period,len(current_anchors),endpoint=False) + current_anchors = np.append(current_anchors, current_anchors[0]+period) + new_anchors = np.append(new_anchors, new_anchors[0]+period) + return np.interp(phis, current_anchors, new_anchors, period=round(period,14)) + + +def angle_between(start_angle, end_angle, check_angle, period=2*np.pi): + start_angle, end_angle, check_angle = start_angle % period, end_angle % period, check_angle % period + if start_angle > end_angle: + adjustment = period - start_angle + start_angle += adjustment + end_angle += adjustment + check_angle += adjustment + start_angle, end_angle, check_angle = start_angle % period, end_angle % period, check_angle % period + return (start_angle <= check_angle) * (check_angle < end_angle) + + +def vertical_flip_path(path, overlay): + if isinstance(path, matplotlib.path.Path): + path.vertices[:, 1] = overlay.svgshape[1] - path.vertices[:, 1] + elif isinstance(path, np.ndarray): + path[:, 1] = overlay.svgshape[1] - path[:, 1] + else: + raise ValueError( + "path should be either a matplotlib.path.Path or an np.ndarray") + return path + + +def determine_path_hemi(overlay, path, percent_cutoff=.4): + if isinstance(path, matplotlib.path.Path): + path = path.vertices + middle_x = overlay.svgshape[0]/2 + return np.mean(middle_x < path[:, 0]) >= percent_cutoff + + +def center_of_mass(X): + # calculate center of mass of a closed polygon + x = X[:, 0] + y = X[:, 1] + g = (x[:-1]*y[1:] - x[1:]*y[:-1]) + A = 0.5*g.sum() + cx = ((x[:-1] + x[1:])*g).sum() + cy = ((y[:-1] + y[1:])*g).sum() + return 1./(6*A)*np.array([cx, cy]) + + +def split_axis(ax, fig=None, split='horizontal', gap=0, **kwargs): + if fig is None: + fit = ax.figure + pos = ax.get_position() + ax.axis('off') + width = pos.x1 - pos.x0 + height = pos.y1 - pos.y0 + if split == 'horizontal': + subax0 = fig.add_axes([ + pos.x0, + pos.y0, + width/2 - width*gap/2, + height], + **kwargs) + subax1 = fig.add_axes([ + pos.x0 + width/2 + width*gap/2, + pos.y0, + width/2 - width*gap/2, + height], + **kwargs) + elif split == 'vertical': + subax0 = fig.add_axes([ + pos.x0, + pos.y0, + width, + height/2 - height*gap/2], + **kwargs) + subax1 = fig.add_axes([ + pos.x0, + pos.y0 + height/2 + height*gap/2, + width, + height/2 - height*gap/2], + **kwargs) + return subax0, subax1 + + +def overlay_axis(ax, fig=None, polar=False): + if fig is None: + fig = ax.figure + rect = ax.get_position() + overlay_ax = fig.add_axes(rect, polar=polar, frameon=False) + return overlay_ax + +def get_num_hemi_vertices(overlay, surface_type='flat'): + """Returns the number of vertices in the left and right hemispheres. + + Parameters + ---------- + overlay (cortex.svgoverlay.SVGOverlay): Pycortex overlay for subject + + Returns + ------- + list + List of two values: [n_vertices_left, n_vertices_right] + """ + subject_id = overlay.svgfile.split('/')[-2] + surfs_flat = [polyutils.Surface(*d) + for d in db.get_surf(subject_id, surface_type)] + return [len(s.pts) for s in surfs_flat] + +def colormap_2d( + cmap, + data=None, + vmin0=None, + vmax0=None, + vmin1=None, + vmax1=None, + map_to_uint8=False, +): + """Map values in two dimensions to color according to a 2D color map image + + Parameters + ---------- + data0 : array (1d) + First dimension of data to map + data1 : array (1d) + Second dimension of data to map + cmap : array (3d) + image of values to use for 2D color map + vmin0 : scalar + vmin for first dimension + vmin1 : scalar + vmin for second dimension + vmax0 : scalar + vmax for first dimension + vmax1 : scalar + vmax for second dimension + map_to_uint8 : bool + whether to map color values to uint8 values (0-255) + """ + if isinstance(cmap, str): + # load pycortex 2D colormap + cmapdir = options.config.get('webgl', 'colormaps') + colormaps = os.listdir(cmapdir) + colormaps = sorted([c for c in colormaps if '.png' in c]) + colormaps = dict((c[:-4], os.path.join(cmapdir, c)) for c in colormaps) + cmap = plt.imread(colormaps[cmap]) + + norm0 = Normalize(vmin0, vmax0) + norm1 = Normalize(vmin1, vmax1) + + def func(data): + data0, data1 = data + d0 = np.clip(norm0(data0), 0, 1) + d1 = np.clip(1 - norm1(data1), 0, 1) + dim0 = np.round(d0 * (cmap.shape[1] - 1)) + # Nans in data seemed to cause weird interaction with conversion to int + dim0 = np.nan_to_num(dim0).astype(int) + dim1 = np.round(d1 * (cmap.shape[0] - 1)) + dim1 = np.nan_to_num(dim1).astype(int) + + colored = cmap[dim1.ravel(), dim0.ravel()] + # May be useful to map r, g, b, a values between 0 and 255 + # to avoid problems with diff plotting functions...? + if map_to_uint8: + colored = (colored * 255).astype(np.uint8) + return colored + if isinstance(data, type(None)): + return func + else: + return func(data) + +def distances_from_vertex(overlay, start_idx, target_idxs=None, surf_type='fiducial'): + n_vertices_left, n_vertices_right = get_num_hemi_vertices(overlay) + subject_id = overlay.svgfile.split('/')[-2] + surfs = [polyutils.Surface(*d) for d in db.get_surf(subject_id, surf_type)] + hemi = start_idx >= n_vertices_left + if surf_type=='flat': + diffs = surfs[hemi].pts - surfs[hemi].pts[start_idx-n_vertices_left*hemi] + distances = np.linalg.norm(diffs[...,:2], axis=-1) + else: + distances = surfs[hemi].geodesic_distance(start_idx-hemi*n_vertices_left) + if not isinstance(target_idxs, type(None)): + hemi_target_idxs = [idx-hemi*n_vertices_left for idx in target_idxs if (hemi==0 and idx < n_vertices_left) or (hemi==1 and idx >= n_vertices_left)] + if len(hemi_target_idxs) > 0: + distances = distances[hemi_target_idxs] + return distances + +def angles_from_vertex(overlay, start_idx, target_idxs=None, period=2*np.pi, theta_direction='ccw'): + diff = overlay.coords - overlay.coords[start_idx] + angles = np.arctan2(diff[:,1], diff[:,0]) + period/4 + if theta_direction not in (-1, 0, 'counter-clockwise', 'ccw'): + angles = (-angles)%period + if not isinstance(target_idxs, type(None)): + angles = angles[target_idxs] + return (angles + period*3/4) % period + + +def get_roi_paths(overlay, roi, cleaned=True, filter_zeros=True, vertical_flip=True, overlay_file=None): + if roi in list(overlay.rois.shapes.keys()): + paths = overlay.rois[roi].splines + elif roi in list(overlay.sulci.shapes.keys()): + paths = overlay.sulci[roi].splines + else: + raise ValueError( + f"{roi} is not a valid ROI or sulcus in the overlay file.") + if cleaned: + paths = [path.cleaned() for path in paths] + if filter_zeros: + for path in paths: + path.vertices = nonzero_coords(path.vertices) + if vertical_flip: + paths = [vertical_flip_path(path, overlay) for path in paths] + hemi_paths = [[], []] + for path in paths: + hemi = determine_path_hemi(overlay, path) + hemi_paths[hemi].append(path) + return hemi_paths + +def get_roi_centroids(overlay, roi, return_indices=True): + centroids = [] + for hemi_paths in get_roi_paths(overlay, roi): + centroid = center_of_mass(np.concatenate([path.vertices for path in hemi_paths])) + if return_indices: + centroid = get_closest_vertices(overlay, centroid).item() + centroids.append(centroid) + return centroids + +def get_roi_outlines(overlay, roi, distance_metric='euclidean', return_indices=False): + paths = get_roi_paths(overlay, roi) + outlines=[[],[]] + for hemi in range(2): + for path in paths[hemi]: + path_vertices = path.vertices + if return_indices: + path_vertices = get_closest_vertices(overlay, path_vertices, distance_metric=distance_metric) + outlines[hemi].append(path_vertices) + return outlines + +def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidean'): + """ + Get vertex idxs closest to the points in a set of coords (e.g. from overlay.sulci or overlay.rois) + + Parameters + ---------- + overlay : cortex.svgoverlay.SVGOverlay + Pycortex overlay for subject + coords : matplotlib.path.Path or np.ndarray + Set of coords on flatmap + metric : str, optional + Distance metric to pass to sklearn for choosing closest vertices. + Defaults to 'euclidean'. + + Returns + ------- + np.ndarray + Array of idxs for the vertex closest to each point in coords. + """ + if isinstance(coords, matplotlib.path.Path): + coords = coords.vertices + elif not isinstance(coords, np.ndarray): + raise ValueError("coords should be a matplotlib.path.Path or an np.ndarray") + if coords.ndim == 1: + coords = coords[np.newaxis] + overlay_coords = overlay.coords.copy() + if hemi!='both': + n_vertices_left, n_vertices_right = get_num_hemi_vertices(overlay, 'flat') + if hemi in [0,'l','left']: + overlay_coords = overlay_coords[:n_vertices_left] + elif hemi in [1,'r','right']: + overlay_coords = overlay_coords[n_vertices_left:] + else: + raise ValueError("hemi must be 0, 1, 'left', 'right' or 'both'.") + min_idxs = np.nanargmin(pairwise_distances(coords,overlay_coords,metric=distance_metric), axis=1) + if hemi in [1,'r','right']: + min_idxs += n_vertices_left + return min_idxs + + +def get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, comparison_roi_full=False, distance_metric='euclidean', return_indices=True): + all_coords = overlay.coords + roi_verts = get_roi_verts(overlay, roi, full=roi_full) + comparison_roi_verts = get_roi_verts( + overlay, comparison_roi, full=comparison_roi_full) + nearest_vertices = [] + for hemi in range(2): + hemi_roi_coords = all_coords[roi_verts[hemi]] + hemi_comparison_roi_coords = all_coords[comparison_roi_verts[hemi]] + distances = pairwise_distances( + hemi_roi_coords, hemi_comparison_roi_coords, metric=distance_metric) + nearest_vertex = np.nanargmin(np.nanmin(distances, axis=1)) + if return_indices: + nearest_vertices.append(roi_verts[hemi][nearest_vertex]) + else: + nearest_vertices.append(hemi_roi_coords[:, nearest_vertex]) + return nearest_vertices + + +def get_roi_verts(overlay, roi, full=True, distance_metric='euclidean', overlay_file=None): + """Like the one in cx.utils, but works for sulci and gyri too. + + Parameters + ---------- + overlay (cortex.svgoverlay.SVGOverlay): Pycortex overlay for subject + roi : str, list or None, optional + ROIs to fetch. Can be ROI name (string), a list of ROI names, or + None, in which case all ROIs will be fetched. + mask : bool + if True, return a logical mask across vertices for the roi + if False, return a list of indices for the ROI + + Returns + ------- + roi_dict : dict + Dictionary of {roi name : roi verts}. ROI verts are for both + hemispheres, with right hemisphere vertex numbers sequential + after left hemisphere vertex numbers. + """ + # Get overlays + n_vertices_left, n_vertices_right = get_num_hemi_vertices(overlay) + + if full: + if roi in list(overlay.rois.shapes.keys()): + roi_verts = overlay.rois.get_mask(roi) + elif roi in list(overlay.sulci.shapes.keys()): + roi_verts = overlay.sulci.get_mask(roi) + else: + raise ValueError(f"ERROR: {roi} is not a valid ROI or sulcus.") + else: + paths = get_roi_paths(overlay, roi) + roi_verts = np.concatenate( + [get_closest_vertices(overlay, path, distance_metric=distance_metric) + for path in flatten_list_of_lists(paths)] + ) + hemi_roi_verts = [ + roi_verts[roi_verts < n_vertices_left], + roi_verts[roi_verts >= n_vertices_left], + ] + return hemi_roi_verts + +# Core functions + +def compute_eccentricity_angle_masks(overlay, centroids, eccentricities, angles): + """ + Compute dartboard-style masks for both hemispheres. + + Parameters + ---------- + overlay : str + pycortex overlay file for subject + centroids : list or tuple + Center vertices indices, one for each hemisphere. + eccentricities : list or tuple + Two lists of eccentricities for eccentricity bins, one for each hemisphere. + Should include lower and upper bounds, e.g. (0, 10, 20) will yield two bins from 0-10 and 10-20. + angles : list + Two lists of angles (degrees) for angle bins, one for each hemisphere. + Should include repeat of start angle, e.g., (0, 180, 270, 360) will yield three bins from + 0-180, 180-270, 270-360. + + Returns + ------- + np.ndarray + Vertex masks, one for each bin and hemisphere. + Shape will be (2, len(eccentricities), len(angles), len(vertices)). + """ + n_vertices_left, n_vertices_right = get_num_hemi_vertices(overlay) + vertex_distances = [distances_from_vertex(overlay, centroids[hemi]) for hemi in range(2)] + vertex_angles = [angles_from_vertex(overlay, centroids[hemi], theta_direction=hemi) for hemi in range(2)] + masks = np.zeros((2, np.array(eccentricities).shape[-1]-1, + np.array(angles).shape[-1]-1, n_vertices_left+n_vertices_right), dtype=bool) + for ecc_i in range(np.array(eccentricities).shape[-1]-1): + for angle_i in range(np.array(angles).shape[-1]-1): + angle_mask = angle_between(angles[0][angle_i], angles[0][angle_i+1], vertex_angles[0][:n_vertices_left], period=2*np.pi) + eccentricity_mask = (eccentricities[0][ecc_i] <= vertex_distances[0]) & (vertex_distances[0] < eccentricities[0][ecc_i+1]) + masks[0, ecc_i, angle_i, :n_vertices_left] = angle_mask & eccentricity_mask + angle_mask = angle_between(angles[1][angle_i], angles[1][angle_i+1], vertex_angles[1][n_vertices_left:], period=2*np.pi) + eccentricity_mask = (eccentricities[1][ecc_i] <= vertex_distances[1]) & (vertex_distances[1] < eccentricities[1][ecc_i+1]) + masks[1, ecc_i, angle_i, n_vertices_left:] = angle_mask & eccentricity_mask + return masks.astype(bool) + +def apply_masks(data, masks, mean_func=np.nanmean, cutoff=None): + """Given vertex data and eccentricity-angle masks, extracts and averages data for each mask. + + Args: + data (np.ndarray): Array of vertex data. Can have any number of dimensions, as long as + the last dimension matches that of the masks. + masks np.ndarray: Vertex masks for each bin and hemisphere, such as the output of + get_eccentricity_angle_masks. Shape should be (2, eccentricities, angles, vertices). + mean_func (function, optional): Function to average values of vertices within a bin. + Defaults to np.nanmean. + cutoff (int, float, optional): Cutoff value for minimum number of vertices for a bin. + If int, the threshold will be based on absolute number of vertices. If float, threshold + will be based on percentage of vertices included per bin. + Bins with fewer vertices than the cutoff will return np.nan. Defaults to None. + + Returns: + np.ndarray: Average vertex values per bin. Will be shape (data.shape[:-1], # vertices). + """ + values = np.zeros((*data.shape[:-1], *masks.shape[:-1])) + for hemi in range(masks.shape[0]): + for eccentricity in range(masks.shape[1]): + for angle in range(masks.shape[2]): + mask = masks[hemi, eccentricity, angle] + if isinstance(cutoff, type(None)): + set_to_nan = False + elif isinstance(cutoff, float): + print((~np.isnan(data[..., mask])).mean()) + set_to_nan = (~np.isnan(data[..., mask])).mean() < cutoff + elif isinstance(cutoff, int): + set_to_nan = (~np.isnan(data[..., mask])).sum() < cutoff + if set_to_nan: + values[..., hemi, eccentricity, angle] = np.nan + else: + values[..., hemi, eccentricity, angle] = mean_func( + data[..., mask], axis=-1) + return values + + +def show_dartboard(data, + data2=None, + axis=None, + image_resolution=500, + cmap=plt.cm.inferno, + vmin=None, + vmax=None, + vmin2=None, + vmax2=None, + theta_direction=-1, + show_grid=True, + grid_linewidth=0.5, + grid_linecolor='lightgray'): + """Given values masked by angle and eccentricity, shows them as a 'dartboard'-style + visualization. + + Args: + data (np.ndarray): np.ndarray: Average vertex values per bin. + Should be of shape (# eccentricities, # angles). + data2 (np.ndarray, optional): np.ndarray: Average vertex values per bin. + Should be of shape (# eccentricities, # angles). + axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. + If None, a new axis will be created. Defaults to None. + image_resolution (int, optional): Resolution of dartboard figure. Defaults to 500. + linewidth (float, optional): Linecolor for bin outlines. Defaults to .5. + cmap (matplotlib.colors.ListedColormap, str, optional): Colormap for data. Can be a colormap + object or a string, which will be loaded through matplotlib or pycortex. Can be 2d. + Defaults to plt.cm.inferno. + vmin (int, float, optional): vmin for first data dimension. Defaults to None. + vmax (int, float, optional): vmax for first data dimension. Defaults to None. + vmin2 (int, float, optional): vmin for optional second data dimension. Defaults to None. + vmax2 (int, float, optional): vmax for optional second data dimension. Defaults to None. + theta_direction (int, string, optional): Direction by which data wraps around dartboard. + 'clockwise' = 'cw' = 1 + 'counter-clockwise' = 'ccw' = -1 + Defaults to 1. + show_grid (bool, optional): Whether to show grid between bins. Defaults to True. + + Returns: + matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + """ + max_radius = 1 + data = np.array(data).astype(np.float) + if isinstance(data2, np.ndarray): + data = np.stack([data, data2], axis=0) + else: + data = data[np.newaxis] + if isinstance(cmap, str): + if data.shape[0] == 1: + try: + cmap = plt.get_cmap(cmap) + except ValueError: + cmap = get_cmap(cmap) + elif data.shape[0] == 2: + cmap = colormap_2d(cmap, vmin0=0, vmin1=0, vmax0=1, vmax1=1,) + if vmin is None: + vmin = np.nanmin(data[0]) + if vmax is None: + vmax = np.nanmax(data[0]) + nrm = Normalize(vmin=vmin, vmax=vmax) + data[0] = np.clip(nrm(data[0]), 0, 1) + if data.shape[0] == 2: + if vmin2 is None: + vmin2 = np.nanmin(data[1]) + if vmax2 is None: + vmax2 = np.nanmax(data[1]) + nrm = Normalize(vmin=vmin2, vmax=vmax2) + data[1] = np.clip(nrm(data[1]), 0, 1) + if axis is None: + _, axis = plt.subplots() + n_radii, n_angles = data[0].shape + angles = np.linspace(0, 2 * np.pi, n_angles, endpoint=False) + radii = np.linspace(0, max_radius, n_radii + 1) + dartboard = np.full((image_resolution, image_resolution, 4), np.nan) + # Define image + image_t = np.linspace(-max_radius, max_radius, image_resolution) + image_x, image_y = np.meshgrid(image_t, image_t) + radius = np.sqrt(image_x**2 + image_y**2) + theta = np.arctan2(image_x, image_y) + theta = np.pi * 2-np.mod(theta + np.pi * 5/2, np.pi*2) + # Define lines + x_lines, y_lines, _ = sph2cart( + np.ones_like(angles) * max_radius, np.degrees(angles), np.zeros_like(angles)) + # Loop to fill dartboard + for radius_i, (radius_small, radius_large) in enumerate(zip(radii[:-1], radii[1:])): + for angle_i, angle in enumerate(angles): + ang_i = ((theta-angle) >= 0) + rad_i = (radius <= radius_large) & (radius > radius_small) + patch_color = cmap(data[:, radius_i, angle_i]) + dartboard[ang_i & rad_i, :] = patch_color + if theta_direction in (-1, 'counter-clockwise', 'ccw'): + dartboard = np.flip(dartboard, axis=1) + axis.imshow(dartboard, extent=(-max_radius, + max_radius, -max_radius, max_radius)) + if show_grid: + n_lines = int(n_angles/2) + for i in range(n_lines): + if i % (n_angles/4) == 0: + axis.plot( + [x_lines[i], x_lines[i+n_lines] + ], [y_lines[i], y_lines[i+n_lines]], + grid_linecolor, lw=grid_linewidth*1.25 + ) + else: + axis.plot( + [x_lines[i], x_lines[i+n_lines] + ], [y_lines[i], y_lines[i+n_lines]], + 'grey', lw=grid_linewidth, #but why grey + ) + for radius in radii[:-1]: + circ = plt.Circle((0, 0), radius=radius, + edgecolor='grey', facecolor='none', lw=grid_linewidth) + axis.add_patch(circ) + circ = plt.Circle( + (0, 0), radius=radii[-1]*.99, edgecolor=grid_linecolor, facecolor='none', lw=grid_linewidth*2, linestyle='--') + axis.add_patch(circ) + axis.axis("off") + return axis + +def interpolate_outlines(phi, rho, resolution=50): + new_phi = np.linspace(0,2*np.pi,resolution) + new_rho = np.interp(new_phi, phi, rho, period=2*np.pi) + return new_phi, new_rho + +def get_interpolated_outlines(overlay, outline_roi, center_roi, anchor_angles, geodesic_distances=True, resolution=100): + raw_outlines = get_roi_outlines(overlay, outline_roi) + # Compute coordinates of overall center + center_coords = get_roi_centroids(overlay, center_roi, return_indices=False) + center_idxs = get_roi_centroids(overlay, center_roi, return_indices=True) + # Get outline center coordinates + outline_centroids_coords = get_roi_centroids(overlay, outline_roi, return_indices=False) + interpolated_outlines = [] + for hemi in range(2): + raw_outline = raw_outlines[hemi][0] + raw_outline_centroid = outline_centroids_coords[hemi] + # Subtract overall center and convert to polar based on center roi centroids + phis, rhos = cart2pol(*(raw_outline-center_coords[hemi]).T) + centroid_phi, centroid_rho = cart2pol(*(raw_outline_centroid-center_coords[hemi])) + # Flip right hemi + if hemi: + phis = np.pi - phis + centroid_phi = np.pi - centroid_phi + # Change rhos to geodesic distances + if geodesic_distances: + outline_centroid_idxs = get_closest_vertices(overlay, raw_outline_centroid) + centroid_rho = distances_from_vertex(overlay, center_idxs[hemi], outline_centroid_idxs) + raw_outline_idxs = get_closest_vertices(overlay, raw_outline) + rhos = distances_from_vertex(overlay, center_idxs[hemi], raw_outline_idxs) + # Warp angles based on anchor centroids + warped_phis = warp_phis(phis, anchor_angles[hemi]) + warped_centroid_phi = warp_phis(centroid_phi, anchor_angles[hemi]) + # Convert back to cartesian + x, y = pol2cart(warped_phis, rhos) + centroid_x, centroid_y = pol2cart(warped_centroid_phi, centroid_rho) + # Subtract own centroid + x -= centroid_x + y -= centroid_y + # Convert to polar relative to these centroids + relative_phi, relative_rho = cart2pol(x,y) + # Even angular interpolation + relative_phi_interp, relative_rho_interp = interpolate_outlines(relative_phi, relative_rho, resolution=resolution) + # Convert interpolated points back to cartesian + x_interp, y_interp = pol2cart(relative_phi_interp, relative_rho_interp) + # Add back own centroid + x_interp += centroid_x + y_interp += centroid_y + interpolated_outlines.append(np.stack(cart2pol(x_interp,y_interp), axis=1)) + return interpolated_outlines + +def show_outlines_on_dartboard( + outlines, axis=None, colors=None, polar=True, hemi=0, rmax=None, as_overlay=False, + **plot_kwargs): + """_summary_ + + Args: + outlines (_type_): _description_ + axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. + colors (_type_, optional): List of colors to iterate through when plotting outlines. + If None, will iterate through default pyplot color cycle. Defaults to None. + polar (bool, optional): Whether outlines are in polar coordinates. Defaults to True. + hemi (int, optional): Hemisphere of the brain from which outlines are drawn. + 0 = left, 1 = right. Defaults to 0. + rmax (int, float, optional): polar maximum radius. Defaults to None. + as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid + interfering with previously-plotted data. Defaults to False. + + Returns: + matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + """ + if axis is None: + fig, axis = plt.subplots( + subplot_kw={'projection': 'polar' if polar else 'rectilinear'}) + elif as_overlay: + fig = axis.figure + rect = axis.get_position() + axis = fig.add_axes(rect, polar=polar, frameon=False) + if colors is None: + colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] + elif isinstance(colors, str): + colors = [colors]*len(outlines) + for outline, color in zip(outlines, colors): + axis.plot(*outline.T, c=color, **plot_kwargs) + axis.grid(False) + if polar: + axis.spines['polar'].set_visible(False) + axis.set_xticklabels([]) + axis.set_yticklabels([]) + if polar: + if hemi in ['right', 'rh', 'r', 1]: + axis.set_theta_zero_location('W') + axis.set_theta_direction(-1) + if rmax is not None: + axis.set_rmax(rmax) + return axis + +def dartboards_table( + dataframe, dataframe_2=None, outlines_dataframe=None, show_titles=True, show_lines=True, + **dartboard_kwargs): + """Given dataframe(s) of dartboard values, plots them in dartboard figures arranged similarly + to the dataframe. + + Args: + dataframe (pandas.core.frame.DataFrame): Dataframe of masked dartboard values. + Should be arranged as subjects X experiments or experiments X subjects. + Elements should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and + (#hemis, # eccentricities, # angles) for two-hemi data. + NaN values will not be plotted. + dataframe_2 (pandas.core.frame.DataFrame, optional): + Dataframe of second data dimension. Elements should be of shape + (# eccentricities, # angles) for single-/cross-hemi data, and + (#hemis, # eccentricities, # angles) for two-hemi data. Defaults to None. + outlines_dataframe (pandas.core.frame.DataFrame, optional): Dataframe of ROI outlines in + polar coordinates, to plot over dartboards. Should be one of following shapes: + (2, # points) for single-/cross-hemisphere data + (# hemispheres, 2, # points) for separate hemispheres + Defaults to None. + show_titles (bool, optional): Whether to show column and index/row titles. Defaults to True. + show_lines (bool, optional): Whether to draw lines separating last column and row from rest + of data, assuming these are for cross-experiment/-subject plots. Defaults to True. + + Returns: + matplotlib.figure.Figure: Matplotlib figure in which data is plotted. + """ + cols_list = list(dataframe.columns.values) + rows_list = list(dataframe.index.values) + + fig, axes = plt.subplots( + len(rows_list), len(cols_list), figsize=[3*(len(cols_list)), 3*(len(rows_list))]) + for col_i, row_name in enumerate(rows_list): + for row_i, col_name in enumerate(cols_list): + axis = axes[col_i][row_i] + dartboard_values = dataframe[col_name][row_name] + if isinstance(dartboard_values, type(np.nan)): + axis.set_axis_off() + else: + if dataframe_2 is not None: + dartboard_values_2 = dataframe_2[col_name][row_name] + else: + dartboard_values_2 = None + if dartboard_values.ndim == 2: + show_dartboard(dartboard_values, dartboard_values_2, axis=axis, + **dartboard_kwargs) + if not isinstance(outlines_dataframe, type(None)): + show_outlines_on_dartboard( + outlines_dataframe[col_name][row_name][:,0], + axis, hemi=0, rmax=50, as_overlay=True, linewidth=1, colors='white') + elif dartboard_values.ndim == 3: + axis_left, axis_right = split_axis(axis, fig) + if dartboard_values_2 is None: + dartboard_values_2 = [None, None] + show_dartboard( + dartboard_values[0], dartboard_values_2[0], axis=axis_left, + theta_direction=-1, **dartboard_kwargs) + show_dartboard( + dartboard_values[1], dartboard_values_2[1], axis=axis_right, + theta_direction=1, **dartboard_kwargs) + if outlines_dataframe is not None: + if isinstance(outlines_dataframe[col_name][row_name], (list, tuple, np.ndarray)): + show_outlines_on_dartboard( + outlines_dataframe[col_name][row_name][:,0], + axis_left, hemi=0, rmax=50, as_overlay=True, linewidth=1, colors='white') + show_outlines_on_dartboard( + outlines_dataframe[col_name][row_name][:,1], + axis_right, hemi=1, rmax=50, as_overlay=True, linewidth=1, colors='white') + if show_titles: + for col_name, axis in zip(cols_list, axes[0, :]): + axis.set_title(col_name, fontdict={'size': 18, 'y': 1.05}) + + # if show_column_titles: + for row_name, axis in zip(rows_list, axes[:, 0]): + axis.annotate(row_name, xy=(-.1, .5), xytext=(0, 5), + xycoords='axes fraction', textcoords='offset points', + size=18, ha='center', va='center', rotation=90) + if show_lines: + renderer = fig.canvas.get_renderer() + + def get_bbox(axis): + return axis.get_tightbbox(renderer).transformed(fig.transFigure.inverted()) + bboxes = np.array(list(map(get_bbox, axes.flat)), + mtrans.Bbox).reshape(axes.shape) + + xmax = np.array(list(map(lambda b: b.x1, bboxes.flat)) + ).reshape(axes.shape).max(axis=0) + xmin = np.array(list(map(lambda b: b.x0, bboxes.flat)) + ).reshape(axes.shape).min(axis=0) + line_xs = np.c_[xmax[1:], xmin[:-1]].mean(axis=1) + + ymax = np.array(list(map(lambda b: b.y1, bboxes.flat)) + ).reshape(axes.shape).max(axis=1) + ymin = np.array(list(map(lambda b: b.y0, bboxes.flat)) + ).reshape(axes.shape).min(axis=1) + line_ys = np.c_[ymax[1:], ymin[:-1]].mean(axis=1) + + col_line = plt.Line2D( + [line_xs[-1], line_xs[-1]], [.125, .875], transform=fig.transFigure, color="black") + fig.add_artist(col_line) + + row_line = plt.Line2D([.1, .9], [line_ys[-1], line_ys[-1]], + transform=fig.transFigure, color="black") + fig.add_artist(row_line) + return fig + + +def draw_mask_outlines( + overlay, masks, axis=None, eccentricities=True, angles=True, + outer_line=True, as_overlay=True, **plot_kwargs): + """Given eccentricity-angle masks, fits and draws outlines of bins on the cortex. + + Args: + masks np.ndarray: Vertex masks for each bin and hemisphere, such as the output of + get_eccentricity_angle_masks. Shape should be (2, eccentricities, angles, vertices). + axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. + eccentricities (bool, optional): Whether to draw outlines for separate eccentricities. + Defaults to True. + angles (bool, optional): Whether to draw outlines for separate angles. + Defaults to True. + outer_line (bool, optional): Whether to draw outline of outermost eccentricity. + Defaults to True. + as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid + interfering with previously-plotted data. Defaults to False. + + Returns: + matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + """ + if axis is None: + _, axis = plt.subplots() + elif as_overlay: + axis = overlay_axis(axis) + + def fit_and_draw_hull(masks): + for mask in masks: + coords = overlay.coords[mask == 1] + hull = ConvexHull(coords) + for simplex in hull.simplices: + axis.plot(coords[simplex, 0], + coords[simplex, 1], **plot_kwargs) + if eccentricities: + fit_and_draw_hull(masks[:, :-1].sum(2).reshape(-1, masks.shape[-1])) + if angles: + fit_and_draw_hull(masks.sum(1).reshape(-1, masks.shape[-1])) + if outer_line: + fit_and_draw_hull(masks[:, -1].sum(1).reshape(-1, masks.shape[-1])) + axis.set_xticks([]) + axis.set_yticks([]) + axis.set_xlim(0, overlay.svgshape[0]) + axis.set_ylim(0, overlay.svgshape[1]) + axis.set_aspect(1) + return axis + + +def draw_mask_bins( + overlay, masks, axis=None, eccentricities=True, angles=True, + as_overlay=True, cmap='gray', alpha=.5, **plot_kwargs): + """Given eccentricity-angle masks, visualizes bins on the cortex as alternating colors on the + cortex. + + Args: + + masks (np.ndarray): Vertex masks, one for each bin and hemisphere. + Shape should will be (2, eccentricities, angles, vertices). + axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. + eccentricities (bool, optional): Whether to color bins for separate eccentricities. + Defaults to True. + angles (bool, optional): Whether to color bins for separate angles. + Defaults to True. + as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid + interfering with previously-plotted data. Defaults to False. + cmap (matplotlib.colors.ListedColormap, str, optional): Colormap to set colors for bins, + which will have values of either 0 or 1. Defaults to 'gray'. + alpha (float, optional): Transparency value, for overlaying on flatmaps. Defaults to .5. + + Returns: + matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + """ + if axis is None: + _, axis = plt.subplots() + elif as_overlay: + axis = overlay_axis(axis) + # Sum over hemispheres + masks = masks.sum(0) + if not eccentricities: + masks = masks.sum(0) + mod = 1 + else: + mod = masks.shape[1] + if not angles: + masks = masks.sum(-2) + masks = masks.reshape(-1, masks.shape[-1]) + bins = np.full(masks.shape[-1], np.nan) + for i, mask in enumerate(masks): + bins[mask == 1] = (i + i//mod) % 2 + bins_vx = Vertex(bins, overlay.svgfile.split('/')[-2], vmin=0, vmax=1) + img_grid, extent = quickflat.composite.make_flatmap_image(bins_vx) + axis.imshow(img_grid, extent=extent, alpha=alpha, cmap=cmap, **plot_kwargs) + axis.set_xticks([]) + axis.set_yticks([]) + axis.set_aspect(1) + return axis + + +def draw_anchor_lines( + overlay, center_idxs, anchor_idxs, axis=None, as_overlay=True, **plot_kwargs): + """Draw lines over the cortex stretching from the center of one ROI to the centers of other + 'anchor' ROIs. + + Args: + overlay (str): Subject identifier string, e.g. 'S1fs' + center_roi (str): ROI from which to draw lines. + anchor_rois (list, tuple): ROIs to which to draw lines from center_roi. + axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. + as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid + interfering with previously-plotted data. Defaults to False. + + Returns: + matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + """ + default_kwargs = { + 'color': 'lightgray', + 'linewidth': 1.5, + } + default_kwargs.update(plot_kwargs) + if axis is None: + _, axis = plt.subplots() + elif as_overlay: + axis = overlay_axis(axis) + for hemi in range(2): + center_coords = overlay.coords[center_idxs[hemi]] + for anchor_idx in anchor_idxs[hemi]: + anchor_coords = overlay.coords[anchor_idx] + axis.plot( + [center_coords[0], anchor_coords[0]], + [center_coords[1], anchor_coords[1]], + **default_kwargs) + axis.set_xticks([]) + axis.set_yticks([]) + axis.set_xlim(0, overlay.svgshape[0]) + axis.set_ylim(0, overlay.svgshape[1]) + axis.set_aspect(1) + return axis + + +def overlay_dartboards( + data, data_2=None, axis=None, fig=None, position_x=.5, position_y=.75, scale=.25, **dartboard_kwargs): + """Show dartboard data in a position overlaid on an existing axis. Primarily intended for + overlaying dartboards on corresponding flatmaps. + + Args: + data (np.ndarray): np.ndarray: Average vertex values per bin. + Should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and + (#hemis, # eccentricities, # angles) for two-hemi data. + data_2 (np.ndarray, optional): np.ndarray: Average vertex values per bin. + Should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and + (#hemis, # eccentricities, # angles) for two-hemi data. Defaults to None. + axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. + fig (matplotlib.figure.Figure, optional): Matplotlib figure in which data is plotted. + x (float, optional): X position in range of (0,1) for center of dartboards. Defaults to .5. + y (float, optional): Y position in range of (0,1) for center of dartboards. Defaults to .5. + scale (float, optional): Scale of dartboards relative to axis size. Defaults to .25. + + Returns: + matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + """ + overlaid_axis = axis.figure.add_axes( + mtrans.Bbox([ + [position_x-scale/2, position_y-scale/2], + [position_x+scale/2, position_y+scale/2]]), + frameon=False + ) + + if data.ndim == 2: + show_dartboard(data, data_2, axis=overlaid_axis, + **dartboard_kwargs) + return overlaid_axis + elif data.ndim == 3: + axis_left, axis_right = split_axis(overlaid_axis, fig) + show_dartboard( + data[0], data_2[0] if not isinstance(data_2, type(None)) else None, axis=axis_left, + theta_direction=-1, **dartboard_kwargs) + show_dartboard( + data[1], data_2[1] if not isinstance(data_2, type(None)) else None, axis=axis_right, + theta_direction=1, **dartboard_kwargs) + return axis_left, axis_right \ No newline at end of file From 6b82dfd7eeb6999eee693c8fcebe9da9a8787bfd Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Fri, 18 Aug 2023 14:01:22 -0700 Subject: [PATCH 02/15] ENH: added function to get dartboard data, cache masks --- cortex/dartboards.py | 391 ++++++++++++------ examples/dartboards/dartboard_basic.py | 21 + .../dartboards/dartboard_flatmap_overlay.py | 3 + 3 files changed, 278 insertions(+), 137 deletions(-) create mode 100644 examples/dartboards/dartboard_basic.py create mode 100644 examples/dartboards/dartboard_flatmap_overlay.py diff --git a/cortex/dartboards.py b/cortex/dartboards.py index f663ed36..49e54510 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -16,6 +16,10 @@ from scipy.spatial import ConvexHull # pylint: disable-msg=E0611 from .utils import get_cmap +import time # Temp + +CACHE_DIR = options.config.get('basic', 'cache') + # Utils @@ -64,6 +68,19 @@ def outline_to_polar(outline, centroid, theta_direction='ccw'): return np.stack([phi, rho], axis=0) def interpolate_angles(angles, angle_sub_bins, period=2*np.pi): + """interpolate between anchor angles + + CONSIDER UPDATING. + + Parameters + ---------- + angles : list + list of angles over which to interpolate. Should be ordered, in radians,... + angle_sub_bins : _type_ + _description_ + period : _type_, optional + _description_, by default 2*np.pi + """ def single_interpolate(angles): angles = list(angles) if angles[-1] % period != angles[0] % period: @@ -252,7 +269,25 @@ def func(data): else: return func(data) -def distances_from_vertex(overlay, start_idx, target_idxs=None, surf_type='fiducial'): +def _distances_from_vertex(overlay, start_idx, target_idxs=None, surf_type='fiducial'): + """Compute distance from vertex to XX + + Parameters + ---------- + overlay : SVGOverlay object + svg overlay object for a subject + start_idx : int + start index + target_idxs : list, optional + other indices idk, by default None + surf_type : str, optional + type of surface on which to compute distance, by default 'fiducial' + + Returns + ------- + _type_ + _description_ + """ """""" n_vertices_left, n_vertices_right = get_num_hemi_vertices(overlay) subject_id = overlay.svgfile.split('/')[-2] surfs = [polyutils.Surface(*d) for d in db.get_surf(subject_id, surf_type)] @@ -268,12 +303,32 @@ def distances_from_vertex(overlay, start_idx, target_idxs=None, surf_type='fiduc distances = distances[hemi_target_idxs] return distances -def angles_from_vertex(overlay, start_idx, target_idxs=None, period=2*np.pi, theta_direction='ccw'): +def _angles_from_vertex(overlay, start_idx, target_idxs=None, period=2*np.pi, theta_direction='ccw'): + """compute angle between two points + + Parameters + ---------- + overlay : SVGOverlay object + svg overlay object for a subject + start_idx : int + idk + target_idxs : list, optional + idk, by default None + period : scalar, optional + max angle, by default 2*np.pi + theta_direction : str, optional + direction for positive angles, by default 'ccw' + + Returns + ------- + _type_ + _description_ + """ diff = overlay.coords - overlay.coords[start_idx] angles = np.arctan2(diff[:,1], diff[:,0]) + period/4 if theta_direction not in (-1, 0, 'counter-clockwise', 'ccw'): angles = (-angles)%period - if not isinstance(target_idxs, type(None)): + if target_idxs is not None: angles = angles[target_idxs] return (angles + period*3/4) % period @@ -299,7 +354,24 @@ def get_roi_paths(overlay, roi, cleaned=True, filter_zeros=True, vertical_flip=T hemi_paths[hemi].append(path) return hemi_paths -def get_roi_centroids(overlay, roi, return_indices=True): +def get_roi_centroids(overlay, roi, return_indices=True, cache_dir=None): + """get centroids of a given ROI in both hemispheres + + Parameters + ---------- + overlay : SVGOverlay object + svg for a given subject + roi : str + name for region or sulcus to choose + return_indices : bool, optional + whether to return indices (True) or coordinates (False), + by default True + + Returns + ------- + centroids : list + coordinate or index for centroids + """ centroids = [] for hemi_paths in get_roi_paths(overlay, roi): centroid = center_of_mass(np.concatenate([path.vertices for path in hemi_paths])) @@ -308,7 +380,8 @@ def get_roi_centroids(overlay, roi, return_indices=True): centroids.append(centroid) return centroids -def get_roi_outlines(overlay, roi, distance_metric='euclidean', return_indices=False): + +def get_roi_outlines(overlay, roi, distance_metric='euclidean', return_indices=False, cache_dir=None): paths = get_roi_paths(overlay, roi) outlines=[[],[]] for hemi in range(2): @@ -319,9 +392,9 @@ def get_roi_outlines(overlay, roi, distance_metric='euclidean', return_indices=F outlines[hemi].append(path_vertices) return outlines -def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidean'): - """ - Get vertex idxs closest to the points in a set of coords (e.g. from overlay.sulci or overlay.rois) + +def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidean', cache_dir=None): + """Get vertex idxs nearest to the points in a set of coords (e.g. from overlay.sulci or overlay.rois) Parameters ---------- @@ -330,13 +403,13 @@ def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidea coords : matplotlib.path.Path or np.ndarray Set of coords on flatmap metric : str, optional - Distance metric to pass to sklearn for choosing closest vertices. + Distance metric to pass to sklearn for choosing nearest vertices. Defaults to 'euclidean'. Returns ------- np.ndarray - Array of idxs for the vertex closest to each point in coords. + Array of idxs for the vertex nearest to each point in coords. """ if isinstance(coords, matplotlib.path.Path): coords = coords.vertices @@ -359,10 +432,10 @@ def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidea return min_idxs -def get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, comparison_roi_full=False, distance_metric='euclidean', return_indices=True): +def _get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, comparison_roi_full=False, distance_metric='euclidean', return_indices=True, cache_dir=None): all_coords = overlay.coords - roi_verts = get_roi_verts(overlay, roi, full=roi_full) - comparison_roi_verts = get_roi_verts( + roi_verts = _get_roi_verts(overlay, roi, full=roi_full) + comparison_roi_verts = _get_roi_verts( overlay, comparison_roi, full=comparison_roi_full) nearest_vertices = [] for hemi in range(2): @@ -378,7 +451,7 @@ def get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, comp return nearest_vertices -def get_roi_verts(overlay, roi, full=True, distance_metric='euclidean', overlay_file=None): +def _get_roi_verts(overlay, roi, full=True, distance_metric='euclidean', overlay_file=None): """Like the one in cx.utils, but works for sulci and gyri too. Parameters @@ -447,8 +520,8 @@ def compute_eccentricity_angle_masks(overlay, centroids, eccentricities, angles) Shape will be (2, len(eccentricities), len(angles), len(vertices)). """ n_vertices_left, n_vertices_right = get_num_hemi_vertices(overlay) - vertex_distances = [distances_from_vertex(overlay, centroids[hemi]) for hemi in range(2)] - vertex_angles = [angles_from_vertex(overlay, centroids[hemi], theta_direction=hemi) for hemi in range(2)] + vertex_distances = [_distances_from_vertex(overlay, centroids[hemi]) for hemi in range(2)] + vertex_angles = [_angles_from_vertex(overlay, centroids[hemi], theta_direction=hemi) for hemi in range(2)] masks = np.zeros((2, np.array(eccentricities).shape[-1]-1, np.array(angles).shape[-1]-1, n_vertices_left+n_vertices_right), dtype=bool) for ecc_i in range(np.array(eccentricities).shape[-1]-1): @@ -644,9 +717,9 @@ def get_interpolated_outlines(overlay, outline_roi, center_roi, anchor_angles, g # Change rhos to geodesic distances if geodesic_distances: outline_centroid_idxs = get_closest_vertices(overlay, raw_outline_centroid) - centroid_rho = distances_from_vertex(overlay, center_idxs[hemi], outline_centroid_idxs) + centroid_rho = _distances_from_vertex(overlay, center_idxs[hemi], outline_centroid_idxs) raw_outline_idxs = get_closest_vertices(overlay, raw_outline) - rhos = distances_from_vertex(overlay, center_idxs[hemi], raw_outline_idxs) + rhos = _distances_from_vertex(overlay, center_idxs[hemi], raw_outline_idxs) # Warp angles based on anchor centroids warped_phis = warp_phis(phis, anchor_angles[hemi]) warped_centroid_phi = warp_phis(centroid_phi, anchor_angles[hemi]) @@ -671,22 +744,33 @@ def get_interpolated_outlines(overlay, outline_roi, center_roi, anchor_angles, g def show_outlines_on_dartboard( outlines, axis=None, colors=None, polar=True, hemi=0, rmax=None, as_overlay=False, **plot_kwargs): - """_summary_ + """ + Plot outlines of regions of interest on dartboard plots. - Args: - outlines (_type_): _description_ - axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. - colors (_type_, optional): List of colors to iterate through when plotting outlines. - If None, will iterate through default pyplot color cycle. Defaults to None. - polar (bool, optional): Whether outlines are in polar coordinates. Defaults to True. - hemi (int, optional): Hemisphere of the brain from which outlines are drawn. - 0 = left, 1 = right. Defaults to 0. - rmax (int, float, optional): polar maximum radius. Defaults to None. - as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid - interfering with previously-plotted data. Defaults to False. + Parameters + ---------- + outlines : _type_ + _description_ + axis : matplotlib.axes._subplots.AxesSubplot, optional + Matplotlib axis on which to plot. + colors : _type_, optional + List of colors to iterate through when plotting outlines. If None, will iterate through + the default pyplot color cycle. Defaults to None. + polar : bool, optional + Whether outlines are in polar coordinates. Defaults to True. + hemi : int, optional + Hemisphere of the brain from which outlines are drawn. + 0 = left, 1 = right. Defaults to 0. + rmax : int or float, optional + Polar maximum radius. Defaults to None. + as_overlay : bool, optional + Whether to create a new axis on top of the existing axis, to avoid interfering with + previously-plotted data. Defaults to False. - Returns: - matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + Returns + ------- + matplotlib.axes._subplots.AxesSubplot + Matplotlib axis in which data is plotted. """ if axis is None: fig, axis = plt.subplots( @@ -714,112 +798,22 @@ def show_outlines_on_dartboard( axis.set_rmax(rmax) return axis -def dartboards_table( - dataframe, dataframe_2=None, outlines_dataframe=None, show_titles=True, show_lines=True, - **dartboard_kwargs): - """Given dataframe(s) of dartboard values, plots them in dartboard figures arranged similarly - to the dataframe. +def dartboard_pair(dartboard_data, + axes=None, + show_grid=True, + grid_linecolor='gray', + grid_linewidth=0.5, + rois=None, + + **dartboard_kwargs, - Args: - dataframe (pandas.core.frame.DataFrame): Dataframe of masked dartboard values. - Should be arranged as subjects X experiments or experiments X subjects. - Elements should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and - (#hemis, # eccentricities, # angles) for two-hemi data. - NaN values will not be plotted. - dataframe_2 (pandas.core.frame.DataFrame, optional): - Dataframe of second data dimension. Elements should be of shape - (# eccentricities, # angles) for single-/cross-hemi data, and - (#hemis, # eccentricities, # angles) for two-hemi data. Defaults to None. - outlines_dataframe (pandas.core.frame.DataFrame, optional): Dataframe of ROI outlines in - polar coordinates, to plot over dartboards. Should be one of following shapes: - (2, # points) for single-/cross-hemisphere data - (# hemispheres, 2, # points) for separate hemispheres - Defaults to None. - show_titles (bool, optional): Whether to show column and index/row titles. Defaults to True. - show_lines (bool, optional): Whether to draw lines separating last column and row from rest - of data, assuming these are for cross-experiment/-subject plots. Defaults to True. + ): + # parse dartboard_data (array or vertex object) - Returns: - matplotlib.figure.Figure: Matplotlib figure in which data is plotted. - """ - cols_list = list(dataframe.columns.values) - rows_list = list(dataframe.index.values) - - fig, axes = plt.subplots( - len(rows_list), len(cols_list), figsize=[3*(len(cols_list)), 3*(len(rows_list))]) - for col_i, row_name in enumerate(rows_list): - for row_i, col_name in enumerate(cols_list): - axis = axes[col_i][row_i] - dartboard_values = dataframe[col_name][row_name] - if isinstance(dartboard_values, type(np.nan)): - axis.set_axis_off() - else: - if dataframe_2 is not None: - dartboard_values_2 = dataframe_2[col_name][row_name] - else: - dartboard_values_2 = None - if dartboard_values.ndim == 2: - show_dartboard(dartboard_values, dartboard_values_2, axis=axis, - **dartboard_kwargs) - if not isinstance(outlines_dataframe, type(None)): - show_outlines_on_dartboard( - outlines_dataframe[col_name][row_name][:,0], - axis, hemi=0, rmax=50, as_overlay=True, linewidth=1, colors='white') - elif dartboard_values.ndim == 3: - axis_left, axis_right = split_axis(axis, fig) - if dartboard_values_2 is None: - dartboard_values_2 = [None, None] - show_dartboard( - dartboard_values[0], dartboard_values_2[0], axis=axis_left, - theta_direction=-1, **dartboard_kwargs) - show_dartboard( - dartboard_values[1], dartboard_values_2[1], axis=axis_right, - theta_direction=1, **dartboard_kwargs) - if outlines_dataframe is not None: - if isinstance(outlines_dataframe[col_name][row_name], (list, tuple, np.ndarray)): - show_outlines_on_dartboard( - outlines_dataframe[col_name][row_name][:,0], - axis_left, hemi=0, rmax=50, as_overlay=True, linewidth=1, colors='white') - show_outlines_on_dartboard( - outlines_dataframe[col_name][row_name][:,1], - axis_right, hemi=1, rmax=50, as_overlay=True, linewidth=1, colors='white') - if show_titles: - for col_name, axis in zip(cols_list, axes[0, :]): - axis.set_title(col_name, fontdict={'size': 18, 'y': 1.05}) - - # if show_column_titles: - for row_name, axis in zip(rows_list, axes[:, 0]): - axis.annotate(row_name, xy=(-.1, .5), xytext=(0, 5), - xycoords='axes fraction', textcoords='offset points', - size=18, ha='center', va='center', rotation=90) - if show_lines: - renderer = fig.canvas.get_renderer() - - def get_bbox(axis): - return axis.get_tightbbox(renderer).transformed(fig.transFigure.inverted()) - bboxes = np.array(list(map(get_bbox, axes.flat)), - mtrans.Bbox).reshape(axes.shape) - - xmax = np.array(list(map(lambda b: b.x1, bboxes.flat)) - ).reshape(axes.shape).max(axis=0) - xmin = np.array(list(map(lambda b: b.x0, bboxes.flat)) - ).reshape(axes.shape).min(axis=0) - line_xs = np.c_[xmax[1:], xmin[:-1]].mean(axis=1) - - ymax = np.array(list(map(lambda b: b.y1, bboxes.flat)) - ).reshape(axes.shape).max(axis=1) - ymin = np.array(list(map(lambda b: b.y0, bboxes.flat)) - ).reshape(axes.shape).min(axis=1) - line_ys = np.c_[ymax[1:], ymin[:-1]].mean(axis=1) - - col_line = plt.Line2D( - [line_xs[-1], line_xs[-1]], [.125, .875], transform=fig.transFigure, color="black") - fig.add_artist(col_line) - - row_line = plt.Line2D([.1, .9], [line_ys[-1], line_ys[-1]], - transform=fig.transFigure, color="black") - fig.add_artist(row_line) - return fig + # loop over hemispheres + + # optionally plot ROIs + pass def draw_mask_outlines( @@ -1000,4 +994,127 @@ def overlay_dartboards( show_dartboard( data[1], data_2[1] if not isinstance(data_2, type(None)) else None, axis=axis_right, theta_direction=1, **dartboard_kwargs) - return axis_left, axis_right \ No newline at end of file + return axis_left, axis_right + + +def _get_anchor_string(anchors): + strs = [] + for a in anchors: + if isinstance(a, tuple): + anchor_name, anchor_type = a + else: + anchor_name, anchor_type = a, 'centroid' + strs.append('{}_{}'.format(anchor_name, anchor_type[0])) + anchor_str = 'a' + '-'.join(strs) + return anchor_str + +def get_dartboard_data(vertex_obj, + centroid, + anchors, + mean_func = np.nanmean, + eccentricities=np.linspace(0, 50, 8+1), + recache=False, + verbose=True, + ): + """retrieve dartboard data for a given vertex object and dartboard parameters + + Parameters + ---------- + vertex_obj : _type_ + _description_ + centroid : _type_ + _description_ + anchors : _type_ + _description_ + mean_func : function, optional + The function used to average all vertices in a given bin, by default np.nanmean + Consider what is being averaged to make an appropriate choice here (e.g. + correlations should be Fischer z transformed before averaging, which may + require a custom function) + eccentricities : _type_, optional + _description_, by default np.linspace(0, 50, 8+1) + + Returns + ------- + masks : array + Masks for each bin of dartboard histogram + data : array + One value for each bin of dartboard histogram. Array is + (hemispheres, eccentricities, angles) starting from 0 degrees + for first value (VERIFY ME FIX) + """ + # Subject + subject = vertex_obj.subject + # SVG overlay object + svg = db.get_overlay(subject) + # Check for cached files: + anchor_str = _get_anchor_string(anchors) + max_radii = (50, 50) # FIX ME + n_angles = 0 # FIX ME + n_eccentricities = 0 # FIX ME + surf_type = 'flat' # FIX ME + rad_str = 'mx%d_%d' % tuple(max_radii) + fname = f'{subject}_dartboard_mask_c{centroid}_{anchor_str}_{n_angles}ang_{n_eccentricities}ecc_{rad_str}_{surf_type}.npy' + fpath = os.path.expanduser(os.path.join(CACHE_DIR, subject, 'cache', fname)) + save_result = True + if os.path.exists(fpath) and not recache: + if verbose: + print(f'Loading masks from {fpath}') + masks = np.load(fpath) + else: + # Compute dartboard masks + if verbose: + print("Computing dartboard masks...") + # Compute centroids of each ROI + t0 = time.time() + # Compute & cache centroids + centroids = {} + for j, anchor in enumerate([centroid] + anchors): + if isinstance(anchor, tuple): + anchor_name, anchor_type = anchor + else: + anchor_name, anchor_type = anchor, 'centroid' + if j == 0: + center_name, center_type = anchor_name, anchor_type + if anchor_type == 'nearest': + centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name, cache_dir=cache_dir) + elif anchor_type == 'centroid': + centroids[anchor_name] = get_roi_centroids(svg, anchor_name, cache_dir=cache_dir) + else: + raise ValueError("unknown anchor type specified: %s\n(Must be 'nearest' or 'centroid')"%(anchor_type)) + t1 = time.time() + print('Time to get centroids:', t1 - t0) + + # Angles from center ROI to each of the anchors + anchor_angles_dict = {} + for anchor in anchors: + if isinstance(anchor, tuple): + anchor_name, anchor_type = anchor + else: + anchor_name, anchor_type = anchor, 'centroid' + anchor_angles_dict[anchor_name] = [_angles_from_vertex(svg, centroids[center_name][hemi], centroids[anchor_name][hemi], theta_direction=hemi) for hemi in range(2)] + t2 = time.time() + print('Time to compute angles:', t2 - t1) + + # Compute the masks, based on specified eccentricity bins, angle bins, and the previously-computed variables + masks = compute_eccentricity_angle_masks( + svg, centroids[center_name], + eccentricities=[ + eccentricities, + eccentricities, + ], + angles=[ + interpolate_angles(np.array([v for v in anchor_angles_dict.values()])[ + :, 0], angle_sub_bins=4), + interpolate_angles(np.array([v for v in anchor_angles_dict.values()])[ + :, 1], angle_sub_bins=4) + ], + ) + t3 = time.time() + print('Time to compute masks:', t3-t2) + if verbose: + print(f'Saving dartboard masks to {fpath}') + np.save(fpath, masks) + output = apply_masks(vertex_obj.data, masks, mean_func) + + return masks, output \ No newline at end of file diff --git a/examples/dartboards/dartboard_basic.py b/examples/dartboards/dartboard_basic.py new file mode 100644 index 00000000..79bc59a7 --- /dev/null +++ b/examples/dartboards/dartboard_basic.py @@ -0,0 +1,21 @@ +# Basic dartboard pair +import cortex + +subject = 'S1' +# Get curvature information for this subject +# Note that this is just an easy example of vertex data; +# Vertex data can always be created by creating a Volume +# and mapping it to vertices, like this: +# vol = cx.Volume() +# vx = vol.map() +vx = cortex.db.get_surfinfo(subject, type='curvature') + +center = 'M1H' + +anchors = [('S1H', 'centroid'), + ('M1F', 'centroid'), + ('FEF', 'centroid'), + ('M1M', 'centroid'), + ] + +out_m1h = cortex.dartboards.get_dartboard_data(vx, center, anchors) diff --git a/examples/dartboards/dartboard_flatmap_overlay.py b/examples/dartboards/dartboard_flatmap_overlay.py new file mode 100644 index 00000000..4b9ddc97 --- /dev/null +++ b/examples/dartboards/dartboard_flatmap_overlay.py @@ -0,0 +1,3 @@ +# Dartboard overlay on flatmap + +# WIP \ No newline at end of file From 9c98764272e0f945e7f5c7e2b50346be2a1ca507 Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Wed, 6 Sep 2023 10:54:00 -0700 Subject: [PATCH 03/15] ENH: dartboard outline and mask caching implemented. Still WIP on final dartboard functions and argument-passing --- cortex/dartboards.py | 530 +++++++++++++++++++++++++++++++++---------- 1 file changed, 409 insertions(+), 121 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index 49e54510..fbe7da62 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -1,4 +1,5 @@ # pylint: disable=missing-module-docstring +import cmath from matplotlib import pyplot as plt from matplotlib import transforms as mtrans from matplotlib.colors import Normalize @@ -12,7 +13,7 @@ from . import quickflat from .dataset.views import Vertex - +from scipy import interpolate from scipy.spatial import ConvexHull # pylint: disable-msg=E0611 from .utils import get_cmap @@ -137,6 +138,18 @@ def determine_path_hemi(overlay, path, percent_cutoff=.4): def center_of_mass(X): + """Calculate the center of mass of a closed polygon + + Parameters + ---------- + X : array + array of x, y coordinates that is (n_points x 2) + + Returns + ------- + center_of_mass + array of shape (2,) that contains (x, y) center of mass + """ # calculate center of mass of a closed polygon x = X[:, 0] y = X[:, 1] @@ -149,7 +162,7 @@ def center_of_mass(X): def split_axis(ax, fig=None, split='horizontal', gap=0, **kwargs): if fig is None: - fit = ax.figure + fig = ax.figure pos = ax.get_position() ax.axis('off') width = pos.x1 - pos.x0 @@ -354,7 +367,7 @@ def get_roi_paths(overlay, roi, cleaned=True, filter_zeros=True, vertical_flip=T hemi_paths[hemi].append(path) return hemi_paths -def get_roi_centroids(overlay, roi, return_indices=True, cache_dir=None): +def get_roi_centroids(overlay, roi, return_indices=True): """get centroids of a given ROI in both hemispheres Parameters @@ -381,7 +394,7 @@ def get_roi_centroids(overlay, roi, return_indices=True, cache_dir=None): return centroids -def get_roi_outlines(overlay, roi, distance_metric='euclidean', return_indices=False, cache_dir=None): +def get_roi_outlines(overlay, roi, distance_metric='euclidean', return_indices=False): paths = get_roi_paths(overlay, roi) outlines=[[],[]] for hemi in range(2): @@ -393,7 +406,7 @@ def get_roi_outlines(overlay, roi, distance_metric='euclidean', return_indices=F return outlines -def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidean', cache_dir=None): +def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidean'): """Get vertex idxs nearest to the points in a set of coords (e.g. from overlay.sulci or overlay.rois) Parameters @@ -432,7 +445,7 @@ def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidea return min_idxs -def _get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, comparison_roi_full=False, distance_metric='euclidean', return_indices=True, cache_dir=None): +def _get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, comparison_roi_full=False, distance_metric='euclidean', return_indices=True): all_coords = overlay.coords roi_verts = _get_roi_verts(overlay, roi, full=roi_full) comparison_roi_verts = _get_roi_verts( @@ -452,11 +465,12 @@ def _get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, com def _get_roi_verts(overlay, roi, full=True, distance_metric='euclidean', overlay_file=None): - """Like the one in cx.utils, but works for sulci and gyri too. + """Like the one in cortex.utils, but works for sulci and gyri too. Parameters ---------- - overlay (cortex.svgoverlay.SVGOverlay): Pycortex overlay for subject + overlay : SVGOverlay + Pycortex overlay for subject roi : str, list or None, optional ROIs to fetch. Can be ROI name (string), a list of ROI names, or None, in which case all ROIs will be fetched. @@ -585,33 +599,41 @@ def show_dartboard(data, show_grid=True, grid_linewidth=0.5, grid_linecolor='lightgray'): - """Given values masked by angle and eccentricity, shows them as a 'dartboard'-style - visualization. + """Given values masked by angle and eccentricity, shows them as a 'dartboard'-style visualization. - Args: - data (np.ndarray): np.ndarray: Average vertex values per bin. - Should be of shape (# eccentricities, # angles). - data2 (np.ndarray, optional): np.ndarray: Average vertex values per bin. - Should be of shape (# eccentricities, # angles). - axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. - If None, a new axis will be created. Defaults to None. - image_resolution (int, optional): Resolution of dartboard figure. Defaults to 500. - linewidth (float, optional): Linecolor for bin outlines. Defaults to .5. - cmap (matplotlib.colors.ListedColormap, str, optional): Colormap for data. Can be a colormap - object or a string, which will be loaded through matplotlib or pycortex. Can be 2d. - Defaults to plt.cm.inferno. - vmin (int, float, optional): vmin for first data dimension. Defaults to None. - vmax (int, float, optional): vmax for first data dimension. Defaults to None. - vmin2 (int, float, optional): vmin for optional second data dimension. Defaults to None. - vmax2 (int, float, optional): vmax for optional second data dimension. Defaults to None. - theta_direction (int, string, optional): Direction by which data wraps around dartboard. - 'clockwise' = 'cw' = 1 - 'counter-clockwise' = 'ccw' = -1 - Defaults to 1. - show_grid (bool, optional): Whether to show grid between bins. Defaults to True. + Parameters + ---------- + data : np.ndarray + Average vertex values per bin. Should be of shape (#eccentricities, #angles). + data2 : np.ndarray, optional + Average vertex values per bin. Should be of shape (#eccentricities, #angles). + axis : plt.Axes, optional + Matplotlib axis on which to plot. If None, a new axis will be created. Defaults to None. + image_resolution : int, optional + Resolution of dartboard figure. Defaults to 500. + linewidth : float, optional + Linecolor for bin outlines. Defaults to 0.5. + cmap : plt.colors.ListedColormap or str, optional + Colormap for data. Can be a colormap object or a string, which will be loaded through + matplotlib or pycortex. Can be 2d. Defaults to plt.cm.inferno. + vmin : float, optional + vmin for the first data dimension. Defaults to None. + vmax : float, optional + vmax for the first data dimension. Defaults to None. + vmin2 : float, optional + vmin for the optional second data dimension. Defaults to None. + vmax2 : float, optional + vmax for the optional second data dimension. Defaults to None. + theta_direction : int or str, optional + Direction by which data wraps around the dartboard. 'clockwise' = 'cw' = 1, + 'counter-clockwise' = 'ccw' = -1. Defaults to 1. + show_grid : bool, optional + Whether to show the grid between bins. Defaults to True. - Returns: - matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + Returns + ------- + plt.Axes + Matplotlib axis in which data is plotted. """ max_radius = 1 data = np.array(data).astype(np.float) @@ -696,13 +718,143 @@ def interpolate_outlines(phi, rho, resolution=50): new_rho = np.interp(new_phi, phi, rho, period=2*np.pi) return new_phi, new_rho -def get_interpolated_outlines(overlay, outline_roi, center_roi, anchor_angles, geodesic_distances=True, resolution=100): +def resample_roi_outline( + angles, + distances, + resolution=100, + every_n=5, + even_sampling_over='polar angle', +): + """resample roi outline to smooth outline + + Relies on angles being monotonically increasing + + Parameters + ---------- + angles : _type_ + _description_ + even_sampling_over : str, optional + _description_, by default 'polar angle' + resolution : int, optional + number of points to sample along the ROI border, by default 100 + + Returns + ------- + _type_ + _description_ + """ + delta = 1e-5 + # Periodic linear resample (upsample) to assure values from 0 to 2pi + if even_sampling_over == 'polar angle': + # Get evenly spaced (resampled) version of angles and distances + angles_linear_sampling = np.linspace(0, 2 * np.pi - delta, len(angles)) + dist_linear_sampling = np.interp( + angles_linear_sampling, angles, distances, period=2 * np.pi) + if every_n is None: + return angles_linear_sampling, dist_linear_sampling + # Smooth: Fit a new spline with points sampled at `every_n` + angle_input = angles_linear_sampling[::every_n] + angle_input[-1] = angles_linear_sampling[-1] + dist_input = dist_linear_sampling[::every_n] + dist_input[-1] = dist_linear_sampling[-1] + smooth = interpolate.interp1d( + angle_input, dist_input, kind='cubic', assume_sorted=True, ) + # Get evenly spaced target angles + angles_out = np.linspace(0, 2 * np.pi - delta, resolution) + dist_out = smooth(angles_out) + return angles_out, dist_out + + if even_sampling_over == 'path length': + ## Resample + xy = np.array(pol2cart(angles, distances)).T + path_dist = np.cumsum( + np.hstack([0, np.linalg.norm(np.diff(xy, axis=0), axis=1)])) + max_dist = np.max(path_dist) + interp = interpolate.interp1d( + path_dist, xy, axis=0, kind='linear', ) + # TO DO: sample evenly between anchor points? + regular_distances = np.linspace( + 0, max_dist, len(angles)) + xy_resampled = interp(regular_distances) + + if every_n is None: + # Map back to angles for output + angles_out, dist_out = cart2pol(*xy_resampled.T) + return angles_out, dist_out + ## Smooth + # Inputs + xy_input = xy_resampled[::every_n] + xy_input[-1] = xy_resampled[-1] + dist_input = regular_distances[::every_n] + dist_input[-1] = regular_distances[-1] + # Smooth x,y values + smooth = interpolate.interp1d(dist_input, + xy_input, + axis=0, kind='cubic', ) + path_dist_out = np.linspace(0, max_dist, resolution) + xy_out = smooth(path_dist_out) + # Convert back to angles + angles_out, dist_out = cart2pol(*xy_out.T) + return angles_out, dist_out + + +def _get_interpolated_outlines(overlay, + outline_roi, + geodesic_distances=True, + resolution=100, + even_sampling_over='polar angle', + every_n=5, + recache=False, + verbose=False, + **dartboard_spec): + """compute outline of roi for a given dartboard plot + + Parameters + ---------- + overlay : SVGOverlay object + svgoverlay for a given subject + outline_roi : str + name of ROI to plot (must sexist in overlays.svg) + center_roi : str + name of center ROI for dartboard + anchor_angles : array + array of angles ?? + geodesic_distances : bool, optional + Flag for whether to compute geodesic distances for eccentricity bins, + by default True + resolution : int, optional + number of points estimated for outline, by default 100 + + Returns + ------- + dict + dictionary of roi names (keys) & outlines (x, y coordinates of boundary) + """ + center_roi = dartboard_spec['center_roi'] + subject = overlay.svgfile.split('/')[-2] + dartboard_str = _get_dartboard_str(**dartboard_spec) + fname = f'{subject}_dartboard_{outline_roi}_outline_{dartboard_str}_resampleby_{even_sampling_over}_sm{every_n}_res{resolution}_geo{geodesic_distances}.npy' + fpath = os.path.expanduser(os.path.join(CACHE_DIR, subject, 'cache', fname)) + if os.path.exists(fpath) and not recache: + if verbose: + print(f'Loading {center_roi} outline from {fpath}') + interpolated_outlines = np.load(fpath) + return interpolated_outlines + else: + if verbose: + print(f'Computing dartboard ROI outline for {center_roi}') + # Compute interpolated outlines, no cache present raw_outlines = get_roi_outlines(overlay, outline_roi) # Compute coordinates of overall center center_coords = get_roi_centroids(overlay, center_roi, return_indices=False) center_idxs = get_roi_centroids(overlay, center_roi, return_indices=True) # Get outline center coordinates outline_centroids_coords = get_roi_centroids(overlay, outline_roi, return_indices=False) + # Get centroids and angles to anchors + center_name, centroids, anchor_angles_dict = _compute_centroids_angles_from_spec( + overlay, verbose=verbose, **dartboard_spec) + anchor_angles = np.array([anchor_angles_dict[roi[0]] for roi in dartboard_spec['anchors']]).T + # Above still has redundancy with getting ROI centroids... interpolated_outlines = [] for hemi in range(2): raw_outline = raw_outlines[hemi][0] @@ -732,13 +884,27 @@ def get_interpolated_outlines(overlay, outline_roi, center_roi, anchor_angles, g # Convert to polar relative to these centroids relative_phi, relative_rho = cart2pol(x,y) # Even angular interpolation - relative_phi_interp, relative_rho_interp = interpolate_outlines(relative_phi, relative_rho, resolution=resolution) + ## CHANGED HERE: May require re-ordering of angles, may rest on assumption + ## about what angle is first + #relative_phi_interp, relative_rho_interp = interpolate_outlines(relative_phi, relative_rho, resolution=resolution) + relative_phi_interp, relative_rho_interp = resample_roi_outline( + relative_phi, + relative_rho, + resolution=resolution, + every_n=every_n, + even_sampling_over=even_sampling_over, + ) # Convert interpolated points back to cartesian x_interp, y_interp = pol2cart(relative_phi_interp, relative_rho_interp) # Add back own centroid x_interp += centroid_x y_interp += centroid_y interpolated_outlines.append(np.stack(cart2pol(x_interp,y_interp), axis=1)) + # cache result + if verbose: + print(f'Saving {center_roi} outline to {fpath}') + np.save(fpath, interpolated_outlines) + return interpolated_outlines def show_outlines_on_dartboard( @@ -798,22 +964,50 @@ def show_outlines_on_dartboard( axis.set_rmax(rmax) return axis -def dartboard_pair(dartboard_data, - axes=None, - show_grid=True, - grid_linecolor='gray', - grid_linewidth=0.5, - rois=None, - - **dartboard_kwargs, - +# Need all inputs to spec dartboard. Make dartboard_spec a thing? a dict? +def show_dartboard_pair(dartboard_data, + dartboard_data2=None, + axes=None, + cmap=plt.cm.viridis, + show_grid=True, + grid_linecolor='gray', + grid_linewidth=0.5, + image_resolution=500, + rois=None, + figsize=(8,4), + **dartboard_spec, ): - # parse dartboard_data (array or vertex object) - - # loop over hemispheres - - # optionally plot ROIs - pass + if dartboard_data2 is None: + data_to_plot = (dartboard_data, (None, None)) + else: + data_to_plot = (dartboard_data, dartboard_data2) + if axes is None: + fig, axes = plt.subplots(1, 2, figsize=figsize) + # Loop over hemispheres + # Hemisphere index goes (0, 1) = (left, right) + directions = [-1, 1] + for hemi_index, data in enumerate(zip(*data_to_plot)): + axis = axes[hemi_index] + d0, d1 = data + if not isinstance(d0, np.ndarray): + d0 = get_dartboard_data(d0, **dartboard_spec) + if d1 is not None: + d1 = get_dartboard_data(d1, **dartboard_spec) + _ = show_dartboard(d0, data2=d1, + axis=axis, + theta_direction=directions[hemi_index], + image_resolution=image_resolution, + cmap=cmap, + show_grid=show_grid, + grid_linewidth=grid_linewidth, + grid_linecolor=grid_linecolor, + **dartboard_spec, + ) + # optionally plot ROIs + if rois is not None: + dartboard_str = _get_dartboard_str(dartboard_spec) + # _get_interpolated_outlines() + return axes def draw_mask_outlines( @@ -866,26 +1060,32 @@ def fit_and_draw_hull(masks): def draw_mask_bins( overlay, masks, axis=None, eccentricities=True, angles=True, as_overlay=True, cmap='gray', alpha=.5, **plot_kwargs): - """Given eccentricity-angle masks, visualizes bins on the cortex as alternating colors on the - cortex. - - Args: + """ + Given eccentricity-angle masks, visualizes bins on the cortex as alternating colors on the cortex. - masks (np.ndarray): Vertex masks, one for each bin and hemisphere. - Shape should will be (2, eccentricities, angles, vertices). - axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. - eccentricities (bool, optional): Whether to color bins for separate eccentricities. - Defaults to True. - angles (bool, optional): Whether to color bins for separate angles. - Defaults to True. - as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid - interfering with previously-plotted data. Defaults to False. - cmap (matplotlib.colors.ListedColormap, str, optional): Colormap to set colors for bins, - which will have values of either 0 or 1. Defaults to 'gray'. - alpha (float, optional): Transparency value, for overlaying on flatmaps. Defaults to .5. + Parameters + ---------- + masks : np.ndarray + Vertex masks, one for each bin and hemisphere. + Shape should be (2, eccentricities, angles, vertices). + axis : matplotlib.axes._subplots.AxesSubplot or None, optional + Matplotlib axis on which to plot. + eccentricities : bool, optional + Whether to color bins for separate eccentricities. Defaults to True. + angles : bool, optional + Whether to color bins for separate angles. Defaults to True. + as_overlay : bool, optional + Whether to create a new axis on top of the existing axis, to avoid interfering with previously-plotted data. + Defaults to False. + cmap : matplotlib.colors.ListedColormap or str, optional + Colormap to set colors for bins, which will have values of either 0 or 1. Defaults to 'gray'. + alpha : float, optional + Transparency value, for overlaying on flatmaps. Defaults to 0.5. - Returns: - matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + Returns + ------- + matplotlib.axes._subplots.AxesSubplot + Matplotlib axis in which data is plotted. """ if axis is None: _, axis = plt.subplots() @@ -954,6 +1154,62 @@ def draw_anchor_lines( return axis +def generate_dartboard_vertex_object(dartboard_mask, + subject, + type='grid', + bg_value=np.nan, + cmap='gray'): + """Make a vertex object containing index values for parts of the dartboard grid + + Parameters + ---------- + dartboard_mask : array + array of masks for each region of dartboard grid, of size + (2, n_eccentricites, n_angles, n_vertices). 2 is hemispheres. + subject : str + pycortex subject surface ID + type : str + one of: 'grid', 'solid', 'eccentricity', 'angle' + n_eccentricities : scalar int + number of eccentricity bins + n_angles : scalar int + number of radial bins + bg_value : float + value for non-dartboard locations + cmap : str + color map to use for pycortex vertex object returned + + TODO: + ---- + Optionally split hemispheres? + """ + assert dartboard_mask.ndim == 4, \ + 'dartboard mask must be (2 x eccentricities x angles x vertices)' + _, n_eccentricities, n_angles, n_vertices = dartboard_mask.shape + mask_img = np.zeros((n_vertices,), ) * bg_value + for ecc_i in range(n_eccentricities): + #ang_i = j % n_angles + #ecc_i = j // n_angles + for ang_i in range(n_angles): + if type == 'grid': + if ecc_i % 2: + value = ang_i % 2 + else: + value = 1 - (ang_i % 2) + elif type == 'solid': + value = 1 + elif type == 'eccentricity': + value = ecc_i + elif type == 'angle': + # Maybe revisit me + value = ang_i + #value = ang_i + for hem in [0, 1]: # left, right + mask_img[dartboard_mask[hem, ecc_i, ang_i, :]] = value + vx = Vertex(mask_img, subject, vmin=0, + vmax=np.nanmax(mask_img), cmap=cmap) + return vx + def overlay_dartboards( data, data_2=None, axis=None, fig=None, position_x=.5, position_y=.75, scale=.25, **dartboard_kwargs): """Show dartboard data in a position overlaid on an existing axis. Primarily intended for @@ -997,22 +1253,76 @@ def overlay_dartboards( return axis_left, axis_right -def _get_anchor_string(anchors): +def _get_dartboard_str(**dartboard_spec): + mx = dartboard_spec['max_radii'] + if not isinstance(mx, tuple): + mx = (mx, mx) strs = [] - for a in anchors: + for a in dartboard_spec['anchors']: if isinstance(a, tuple): anchor_name, anchor_type = a else: anchor_name, anchor_type = a, 'centroid' strs.append('{}_{}'.format(anchor_name, anchor_type[0])) - anchor_str = 'a' + '-'.join(strs) - return anchor_str + anchor_str = 'a' + '_'.join(strs) + fmt = dict( + center_roi=dartboard_spec['center_roi'], + max_radii=mx, + n_angles=dartboard_spec['n_angles'], + n_eccentricities=dartboard_spec['n_eccentricities'], + anchor_str=anchor_str, + rad_str='mx%d_%d' % tuple(mx), + ) + fname = '_center{center_roi}-anchors_{anchor_str}-{n_angles}ang-{n_eccentricities}ecc-{rad_str}' + return fname.format(**fmt) + +def _compute_centroids_angles_from_spec(svg, center_roi, anchors, verbose=False, **kwargs): + """kwargs catches extra inputs from dartboard_spec. + Perhaps not the most principled. + functional. + """ + # Compute centroids of each ROI + t0 = time.time() + # Compute centroids + centroids = {} + for j, anchor in enumerate([center_roi] + anchors): + if isinstance(anchor, tuple): + anchor_name, anchor_type = anchor + else: + anchor_name, anchor_type = anchor, 'centroid' + if j == 0: + center_name, center_type = anchor_name, anchor_type + if anchor_type == 'nearest': + centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name) + elif anchor_type == 'centroid': + centroids[anchor_name] = get_roi_centroids(svg, anchor_name) + else: + raise ValueError("unknown anchor type specified: %s\n(Must be 'nearest' or 'centroid')"%(anchor_type)) + t1 = time.time() + if verbose: + print('Time to get centroids:', t1 - t0) + + # Compute angles from center ROI to each of the anchors + anchor_angles_dict = {} + for anchor in anchors: + if isinstance(anchor, tuple): + anchor_name, anchor_type = anchor + else: + anchor_name, anchor_type = anchor, 'centroid' + anchor_angles_dict[anchor_name] = [_angles_from_vertex(svg, centroids[center_name][hemi], centroids[anchor_name][hemi], theta_direction=hemi) for hemi in range(2)] + t2 = time.time() + if verbose: + print('Time to compute angles:', t2 - t1) + return center_name, centroids, anchor_angles_dict def get_dartboard_data(vertex_obj, - centroid, + center_roi, anchors, - mean_func = np.nanmean, - eccentricities=np.linspace(0, 50, 8+1), + max_radii=(50, 50), + n_angles=16, + n_eccentricities=8, + eccentricities=None, + mean_func=np.nanmean, recache=False, verbose=True, ): @@ -1022,7 +1332,7 @@ def get_dartboard_data(vertex_obj, ---------- vertex_obj : _type_ _description_ - centroid : _type_ + center_roi : _type_ _description_ anchors : _type_ _description_ @@ -1032,7 +1342,7 @@ def get_dartboard_data(vertex_obj, correlations should be Fischer z transformed before averaging, which may require a custom function) eccentricities : _type_, optional - _description_, by default np.linspace(0, 50, 8+1) + _description_, by default np.linspace(0, max_radii, 8+1) Returns ------- @@ -1047,16 +1357,22 @@ def get_dartboard_data(vertex_obj, subject = vertex_obj.subject # SVG overlay object svg = db.get_overlay(subject) + # Handle max_radii (should be tuple) + if not isinstance(max_radii, tuple): + max_radii = (max_radii, max_radii) + # Allow manually specified eccentricities to override linear spacing + if eccentricities is None: + eccentricities = [np.linspace(0, mr, n_eccentricities + 1) for mr in max_radii] # Check for cached files: - anchor_str = _get_anchor_string(anchors) - max_radii = (50, 50) # FIX ME - n_angles = 0 # FIX ME - n_eccentricities = 0 # FIX ME - surf_type = 'flat' # FIX ME - rad_str = 'mx%d_%d' % tuple(max_radii) - fname = f'{subject}_dartboard_mask_c{centroid}_{anchor_str}_{n_angles}ang_{n_eccentricities}ecc_{rad_str}_{surf_type}.npy' + dartboard_spec = dict(center_roi=center_roi, + anchors=anchors, + n_angles=n_angles, + n_eccentricities=n_eccentricities, + eccentricities=eccentricities, + max_radii=max_radii,) + dartboard_str = _get_dartboard_str(**dartboard_spec) + fname = f'{subject}_dartboard_mask_{dartboard_str}.npy' fpath = os.path.expanduser(os.path.join(CACHE_DIR, subject, 'cache', fname)) - save_result = True if os.path.exists(fpath) and not recache: if verbose: print(f'Loading masks from {fpath}') @@ -1065,44 +1381,14 @@ def get_dartboard_data(vertex_obj, # Compute dartboard masks if verbose: print("Computing dartboard masks...") - # Compute centroids of each ROI - t0 = time.time() - # Compute & cache centroids - centroids = {} - for j, anchor in enumerate([centroid] + anchors): - if isinstance(anchor, tuple): - anchor_name, anchor_type = anchor - else: - anchor_name, anchor_type = anchor, 'centroid' - if j == 0: - center_name, center_type = anchor_name, anchor_type - if anchor_type == 'nearest': - centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name, cache_dir=cache_dir) - elif anchor_type == 'centroid': - centroids[anchor_name] = get_roi_centroids(svg, anchor_name, cache_dir=cache_dir) - else: - raise ValueError("unknown anchor type specified: %s\n(Must be 'nearest' or 'centroid')"%(anchor_type)) - t1 = time.time() - print('Time to get centroids:', t1 - t0) - - # Angles from center ROI to each of the anchors - anchor_angles_dict = {} - for anchor in anchors: - if isinstance(anchor, tuple): - anchor_name, anchor_type = anchor - else: - anchor_name, anchor_type = anchor, 'centroid' - anchor_angles_dict[anchor_name] = [_angles_from_vertex(svg, centroids[center_name][hemi], centroids[anchor_name][hemi], theta_direction=hemi) for hemi in range(2)] - t2 = time.time() - print('Time to compute angles:', t2 - t1) + center_name, centroids, anchor_angles_dict = _compute_centroids_angles_from_spec( + svg, verbose=verbose, **dartboard_spec) # Compute the masks, based on specified eccentricity bins, angle bins, and the previously-computed variables + t2 = time.time() masks = compute_eccentricity_angle_masks( svg, centroids[center_name], - eccentricities=[ - eccentricities, - eccentricities, - ], + eccentricities=eccentricities, angles=[ interpolate_angles(np.array([v for v in anchor_angles_dict.values()])[ :, 0], angle_sub_bins=4), @@ -1111,10 +1397,12 @@ def get_dartboard_data(vertex_obj, ], ) t3 = time.time() - print('Time to compute masks:', t3-t2) if verbose: + print('Time to compute masks:', t3-t2) print(f'Saving dartboard masks to {fpath}') np.save(fpath, masks) output = apply_masks(vertex_obj.data, masks, mean_func) - return masks, output \ No newline at end of file + return masks, output + + From 67775784b5d9b1a5c69c94f2ad45a64a13831143 Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Thu, 12 Oct 2023 15:55:29 -0700 Subject: [PATCH 04/15] ENH: Functioning dartboard pair plotting --- cortex/dartboards.py | 92 ++++++++++++++++++++++++++++++++++++-------- 1 file changed, 77 insertions(+), 15 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index fbe7da62..01117bef 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -7,6 +7,7 @@ import matplotlib import numpy as np import os + from . import db from . import options from . import polyutils @@ -967,32 +968,63 @@ def show_outlines_on_dartboard( # Need all inputs to spec dartboard. Make dartboard_spec a thing? a dict? def show_dartboard_pair(dartboard_data, dartboard_data2=None, + subject=None, axes=None, + rois=None, + # plot parameters cmap=plt.cm.viridis, + vmin=None, + vmax=None, + vmin2=None, + vmax2=None, + image_resolution=500, + figsize=(8, 4), + # grid line parameters show_grid=True, grid_linecolor='gray', grid_linewidth=0.5, - image_resolution=500, - rois=None, - figsize=(8,4), + # roi path line parameters + geodesic_distances=True, + path_resolution=100, + every_n=5, + even_sampling_over='path length', + roi_linewidth=1, + roi_color='r', + recache=False, + verbose=False, **dartboard_spec, ): + if subject is None: + if not isinstance(dartboard_data, Vertex): + raise ValueError("Must provide either `subject` argument or a Vertex object as `dartboard_data`") + subject = dartboard_data.subject + # WORKING HERE. need to re-figure-out what data_to_plot is. + if isinstance(dartboard_data, Vertex): + _, data0 = get_dartboard_data(dartboard_data, **dartboard_spec) + else: + data0 = dartboard_data.copy() if dartboard_data2 is None: - data_to_plot = (dartboard_data, (None, None)) + # None for left hem, None for right hem + data1 = (None, None) + elif isinstance(dartboard_data2, Vertex): + _, data1 = get_dartboard_data(dartboard_data2, **dartboard_spec) else: - data_to_plot = (dartboard_data, dartboard_data2) + data1 = dartboard_data2.copy() + + data_to_plot = (dartboard_data, dartboard_data2) + if axes is None: fig, axes = plt.subplots(1, 2, figsize=figsize) # Loop over hemispheres # Hemisphere index goes (0, 1) = (left, right) directions = [-1, 1] - for hemi_index, data in enumerate(zip(*data_to_plot)): + for hemi_index, data in enumerate(zip(data0, data1)): axis = axes[hemi_index] d0, d1 = data - if not isinstance(d0, np.ndarray): - d0 = get_dartboard_data(d0, **dartboard_spec) - if d1 is not None: - d1 = get_dartboard_data(d1, **dartboard_spec) + # if not isinstance(d0, np.ndarray): + # d0 = get_dartboard_data(d0, **dartboard_spec) + # if d1 is not None: + # d1 = get_dartboard_data(d1, **dartboard_spec) _ = show_dartboard(d0, data2=d1, axis=axis, theta_direction=directions[hemi_index], @@ -1001,12 +1033,42 @@ def show_dartboard_pair(dartboard_data, show_grid=show_grid, grid_linewidth=grid_linewidth, grid_linecolor=grid_linecolor, - **dartboard_spec, + vmin=vmin, + vmax=vmax, + vmin2=vmin2, + vmax2=vmax2, ) - # optionally plot ROIs - if rois is not None: - dartboard_str = _get_dartboard_str(dartboard_spec) - # _get_interpolated_outlines() + # optionally plot ROIs + if rois is not None: + overlay = db.get_overlay(subject) + outlines = {} + if not isinstance(roi_color, (list, tuple)): + roi_color = [roi_color] * len(rois) + for roi in rois: + outlines[roi] = _get_interpolated_outlines(overlay, + roi, + geodesic_distances=geodesic_distances, + resolution=path_resolution, + even_sampling_over=even_sampling_over, + every_n=every_n, + recache=recache, + verbose=verbose, + **dartboard_spec) + for hemi_index in range(2): + hemi_axis = axes[hemi_index] + outlines_to_plot = np.array([v for v in outlines.values()])[:, hemi_index] + if isinstance(dartboard_spec['max_radii'], tuple): + rmax = dartboard_spec['max_radii'][hemi_index] + else: + rmax = dartboard_spec['max_radii'] + show_outlines_on_dartboard( + outlines_to_plot, + hemi_axis, + hemi=hemi_index, + rmax=rmax, + as_overlay=True, # As overlay + linewidth=roi_linewidth, + colors=roi_color) return axes From 3c66a54b41c5356d4a5d44a393cab4277f377ab1 Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Tue, 17 Oct 2023 21:55:08 -0700 Subject: [PATCH 05/15] ENH: updated docstrings, copied old code version of dartboard plotted over flatmap, this needs update --- cortex/dartboards.py | 447 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 393 insertions(+), 54 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index 01117bef..ead53be3 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -1,11 +1,10 @@ -# pylint: disable=missing-module-docstring -import cmath from matplotlib import pyplot as plt from matplotlib import transforms as mtrans from matplotlib.colors import Normalize from sklearn.metrics import pairwise_distances import matplotlib import numpy as np +import time import os from . import db @@ -15,16 +14,15 @@ from .dataset.views import Vertex from scipy import interpolate -from scipy.spatial import ConvexHull # pylint: disable-msg=E0611 +from scipy.spatial import ConvexHull from .utils import get_cmap -import time # Temp - -CACHE_DIR = options.config.get('basic', 'cache') - -# Utils - +try: + CACHE_DIR = options.config.get('basic', 'cache') +except: + CACHE_DIR = options.config.get('basic','filestore') +## Utils def flatten_list_of_lists(list_of_lists): flat_list = [item for sublist in list_of_lists for item in sublist] return flat_list @@ -54,6 +52,7 @@ def cart2pol(x, y): rho = np.sqrt(x**2 + y**2) return phi, rho + def pol2cart(phi, rho): x = rho * np.cos(phi) y = rho * np.sin(phi) @@ -63,12 +62,14 @@ def pol2cart(phi, rho): def nonzero_coords(coords): return coords[np.any(coords != 0, axis=1)] + def outline_to_polar(outline, centroid, theta_direction='ccw'): phi, rho = cart2pol(*(outline.T-centroid[:, np.newaxis])) if theta_direction not in (-1, 0, 'counter-clockwise', 'ccw'): phi = np.pi - phi return np.stack([phi, rho], axis=0) + def interpolate_angles(angles, angle_sub_bins, period=2*np.pi): """interpolate between anchor angles @@ -101,6 +102,7 @@ def single_interpolate(angles): else: return single_interpolate(angles) + def warp_phis(phis, current_anchors, new_anchors=None, period=2*np.pi): if isinstance(new_anchors, type(None)): new_anchors = np.linspace(0,period,len(current_anchors),endpoint=False) @@ -204,6 +206,7 @@ def overlay_axis(ax, fig=None, polar=False): overlay_ax = fig.add_axes(rect, polar=polar, frameon=False) return overlay_ax + def get_num_hemi_vertices(overlay, surface_type='flat'): """Returns the number of vertices in the left and right hemispheres. @@ -221,6 +224,7 @@ def get_num_hemi_vertices(overlay, surface_type='flat'): for d in db.get_surf(subject_id, surface_type)] return [len(s.pts) for s in surfs_flat] + def colormap_2d( cmap, data=None, @@ -508,8 +512,7 @@ def _get_roi_verts(overlay, roi, full=True, distance_metric='euclidean', overlay ] return hemi_roi_verts -# Core functions - +## Core functions def compute_eccentricity_angle_masks(overlay, centroids, eccentricities, angles): """ Compute dartboard-style masks for both hemispheres. @@ -549,23 +552,25 @@ def compute_eccentricity_angle_masks(overlay, centroids, eccentricities, angles) masks[1, ecc_i, angle_i, n_vertices_left:] = angle_mask & eccentricity_mask return masks.astype(bool) + def apply_masks(data, masks, mean_func=np.nanmean, cutoff=None): """Given vertex data and eccentricity-angle masks, extracts and averages data for each mask. - Args: - data (np.ndarray): Array of vertex data. Can have any number of dimensions, as long as - the last dimension matches that of the masks. - masks np.ndarray: Vertex masks for each bin and hemisphere, such as the output of - get_eccentricity_angle_masks. Shape should be (2, eccentricities, angles, vertices). - mean_func (function, optional): Function to average values of vertices within a bin. - Defaults to np.nanmean. - cutoff (int, float, optional): Cutoff value for minimum number of vertices for a bin. - If int, the threshold will be based on absolute number of vertices. If float, threshold - will be based on percentage of vertices included per bin. - Bins with fewer vertices than the cutoff will return np.nan. Defaults to None. + Parameters + ---------- + data : np.ndarray + Array of vertex data. Can have any number of dimensions, as long as the last dimension matches that of the masks. + masks : np.ndarray + Vertex masks for each bin and hemisphere, such as the output of get_eccentricity_angle_masks. Shape should be (2, eccentricities, angles, vertices). + mean_func : function, optional + Function to average values of vertices within a bin. Defaults to np.nanmean. + cutoff : int or float, optional + Cutoff value for the minimum number of vertices for a bin. If an int, the threshold will be based on the absolute number of vertices. If a float, the threshold will be based on the percentage of vertices included per bin. Bins with fewer vertices than the cutoff will return np.nan. Defaults to None. - Returns: - np.ndarray: Average vertex values per bin. Will be shape (data.shape[:-1], # vertices). + Returns + ------- + np.ndarray + Average vertex values per bin. Will be of shape (data.shape[:-1], # vertices). """ values = np.zeros((*data.shape[:-1], *masks.shape[:-1])) for hemi in range(masks.shape[0]): @@ -714,11 +719,13 @@ def show_dartboard(data, axis.axis("off") return axis + def interpolate_outlines(phi, rho, resolution=50): new_phi = np.linspace(0,2*np.pi,resolution) new_rho = np.interp(new_phi, phi, rho, period=2*np.pi) return new_phi, new_rho + def resample_roi_outline( angles, distances, @@ -911,8 +918,7 @@ def _get_interpolated_outlines(overlay, def show_outlines_on_dartboard( outlines, axis=None, colors=None, polar=True, hemi=0, rmax=None, as_overlay=False, **plot_kwargs): - """ - Plot outlines of regions of interest on dartboard plots. + """Plot outlines of regions of interest on dartboard plots. Parameters ---------- @@ -994,6 +1000,78 @@ def show_dartboard_pair(dartboard_data, verbose=False, **dartboard_spec, ): + """Make two dartboard plots (one for each hemisphere) of dartboard data + + Parameters + ---------- + dartboard_data : array or pycortex Vertex object + if array, should be (); if vertex object, vertex data will be masked + according to `dartboard_spec` kwargs + dartboard_data2 : array or pycortex Vertex object, optional + optional second dimension of data for 2D colormap plots; + see `dartboard_data` for format, by default None + subject : str, optional + string identifier for pycortex subject; not necessary if a vertex object + is provided for dartboard_data, by default None + axes : matplotlib axis, optional + axis into which to plot; if None, a new plot is created, by default None + rois : list, optional + list of string names for ROIs to plot on dartboard space, by default None + cmap : str or cmap, optional + colormap or name of colormap to use to colormap data. string names for + pycortex colormaps, including 2D colormaps, are allowable, by default plt.cm.viridis + vmin : scalar, optional + vmin for first dimension of data to be plotted on dartboard, by default None + vmax : scalar, optional + vmax for first dimension of data to be plotted on dartboard, by default None + vmin2 : scalar, optional + vmin for second dimension of data to be plotted on dartboard, by default None + vmax2 : scalar, optional + vmax for second dimension of data to be plotted on dartboard, by default None + image_resolution : int, optional + resolution of dartboard images, by default 500 + figsize : tuple, optional + size in inches of resulting plot; should be twice as wide as it is tall, by default (8, 4) + show_grid : bool, optional + whether to show lines between dartboard bins, by default True + grid_linecolor : str, optional + colorspec for grid lines, by default 'gray' + grid_linewidth : float, optional + width of grid lines, by default 0.5 + geodesic_distances : bool, optional + whether to compute distances along folded cortical surface (True), or to + simply compute distances in flattened space (False). True is slower, by default True + path_resolution : int, optional + resolution of ROI paths plotted on dartboard in points, by default 100 + every_n : int, optional + sampling along ROI path for smoothing; ROI paths are warped into dartboard + space, a little smoothing usually helps aesthetically. 1 is no smoothing + , by default 5 + even_sampling_over : str, optional + How to resample ROI paths, by angle or along path length. 'angle' is perhaps slightly + more principled for convex ROIs, but 'path_length' gives better results for non-convex + ROIs, by default 'path length' + roi_linewidth : int, optional + width of lines for drawn ROIs, by default 1 + roi_color : list or str, optional + color of lines for drawn ROIs; if list is provided, ROIs are each colored + in order of the colors, by default 'r' + recache : bool, optional + recache mask and ROI outline data (data is cached in the pycortex cache + to speed up processing), by default False + verbose : bool, optional + verbose output, by default False + + Returns + ------- + list of matplotlib.axes + list of axes into which data was plotted + + Raises + ------ + ValueError + _description_ + """ """""" if subject is None: if not isinstance(dartboard_data, Vertex): raise ValueError("Must provide either `subject` argument or a Vertex object as `dartboard_data`") @@ -1177,19 +1255,26 @@ def draw_mask_bins( def draw_anchor_lines( overlay, center_idxs, anchor_idxs, axis=None, as_overlay=True, **plot_kwargs): - """Draw lines over the cortex stretching from the center of one ROI to the centers of other - 'anchor' ROIs. + """ + Draw lines over the cortex stretching from the center of one ROI to the centers of other 'anchor' ROIs. - Args: - overlay (str): Subject identifier string, e.g. 'S1fs' - center_roi (str): ROI from which to draw lines. - anchor_rois (list, tuple): ROIs to which to draw lines from center_roi. - axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. - as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid - interfering with previously-plotted data. Defaults to False. + Parameters + ---------- + overlay : str + Subject identifier string, e.g. 'S1fs'. + center_roi : str + ROI from which to draw lines. + anchor_rois : list or tuple + ROIs to which to draw lines from center_roi. + axis : matplotlib.axes._subplots.AxesSubplot, optional + Matplotlib axis on which to plot. + as_overlay : bool, optional + Whether to create a new axis on top of the existing axis to avoid interfering with previously-plotted data. Defaults to False. - Returns: - matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + Returns + ------- + matplotlib.axes._subplots.AxesSubplot + Matplotlib axis in which data is plotted. """ default_kwargs = { 'color': 'lightgray', @@ -1274,24 +1359,31 @@ def generate_dartboard_vertex_object(dartboard_mask, def overlay_dartboards( data, data_2=None, axis=None, fig=None, position_x=.5, position_y=.75, scale=.25, **dartboard_kwargs): - """Show dartboard data in a position overlaid on an existing axis. Primarily intended for + """ + Show dartboard data in a position overlaid on an existing axis. Primarily intended for overlaying dartboards on corresponding flatmaps. - Args: - data (np.ndarray): np.ndarray: Average vertex values per bin. - Should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and - (#hemis, # eccentricities, # angles) for two-hemi data. - data_2 (np.ndarray, optional): np.ndarray: Average vertex values per bin. - Should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and - (#hemis, # eccentricities, # angles) for two-hemi data. Defaults to None. - axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. - fig (matplotlib.figure.Figure, optional): Matplotlib figure in which data is plotted. - x (float, optional): X position in range of (0,1) for center of dartboards. Defaults to .5. - y (float, optional): Y position in range of (0,1) for center of dartboards. Defaults to .5. - scale (float, optional): Scale of dartboards relative to axis size. Defaults to .25. + Parameters + ---------- + data : np.ndarray + Average vertex values per bin. Should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and (#hemis, # eccentricities, # angles) for two-hemi data. + data_2 : np.ndarray, optional + Average vertex values per bin. Should be of shape (# eccentricities, # angles) for single-/cross-hemi data, and (#hemis, # eccentricities, # angles) for two-hemi data. Defaults to None. + axis : matplotlib.axes._subplots.AxesSubplot, optional + Matplotlib axis on which to plot. + fig : matplotlib.figure.Figure, optional + Matplotlib figure in which data is plotted. + x : float, optional + X position in the range of (0,1) for the center of dartboards. Defaults to 0.5. + y : float, optional + Y position in the range of (0,1) for the center of dartboards. Defaults to 0.5. + scale : float, optional + Scale of dartboards relative to the axis size. Defaults to 0.25. - Returns: - matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + Returns + ------- + matplotlib.axes._subplots.AxesSubplot + Matplotlib axis in which data is plotted. """ overlaid_axis = axis.figure.add_axes( mtrans.Bbox([ @@ -1396,8 +1488,25 @@ def get_dartboard_data(vertex_obj, _description_ center_roi : _type_ _description_ - anchors : _type_ - _description_ + anchors : list + list of strings or tuples specifying anchor points. If strings are provided, + code assumes strings to be names for regions of interest defined for this + pycortex subject (e.g. 'FFA'). Anchor points are computed as centroids of + the named ROIs. If tuples are provided, the first element of the tuple is the + name of the ROI, and the second element must be 'centroid' or 'nearest', + indicating whether the anchor point shoudl be the centroi of the named ROI, or + the nearest point in or on the named region to the center_roi. Sulci or other + labeled anatomical markers in other layers of overlays.svg for the subject may + be used, e.g. ('STS', 'nearest') for the nearest point to the center_roi that + falls on STS. + + Important: anchors must be specified counter clockwise from the RIGHT with + respect to the RIGHT hemisphere: + 2 + /\ + 3 <- center -> 1 + \/ + 4 mean_func : function, optional The function used to average all vertices in a given bin, by default np.nanmean Consider what is being averaged to make an appropriate choice here (e.g. @@ -1468,3 +1577,233 @@ def get_dartboard_data(vertex_obj, return masks, output +def dartboard_on_flatmap(vertex_data, + fn=np.nanmean, + vmin=None, + vmax=None, + cmap=None, + # Dartboard args + center_roi=None, + anchor_rois=None, + display_rois=None, + n_angles=16, + n_eccentricities=8, + max_radii=(50, 50), + surf_type='inflated', + # Plotting args + figsize=(12, 6), + dartboard_axes_dist_from_midline=0.15, + dartboard_axes_bottom=0.25, + dartboard_axes_width=0.1, + dartboard_axes_height=0.2, + dartboard_display_alpha=0.2, + quickflat_kw=None, + flatmap_line_linewidth=1.5, + flatmap_line_color='y', + flatmap_line_style=None, + show_anchor_lines=None, + show_dartboard_grid=True, + show_dartboard_edge=True, + #show_roi_outlines=True, + # ROI outline parameters + n_roi_border_points=64, + roi_outline_smooth_factor=5, # every 5th point kept, smoothed with cubic spline + roi_border_kw=None, + cache_dir=None, + verbose=False, + outline_kw=None, + **kwargs): + """Make a flatmap with overlaid dartboard plots + + Parameters + ---------- + vertex_data : cortex.Vertex (or cortex.Vertex2D, or cortex.VertexRGB*) instance + data to be plotted (VertexRGB still WIP). Note that vmin, vmax of dartboard + use vmin, vmax of this data. + masks : array-like, optional + masks for , by default None + fn : function, optional + function to use to collapse over vertices within each mask, by default np.nanmean + center_roi : str + center ROI + n_eccentricities : int, optional + number of radial bins, by default 5; should be consistent with number of radial + bins in `masks`, if `masks` is provided + n_angles : int, optional + number of polar angle bins, by default 8; should be consistent with number of + angular bins in `masks`, if `masks` is provided + dartboard_axes_dist_from_midline : float, optional + fraction of the figure by which dartboard axes are displaced from the figure midline, + by default 0.15 + dartboard_axes_bottom : float, optional + fraction of the figure by which dartboard axes are displaced from the figure bottom, + by default 0.25 + dartboard_axes_width : float, optional + width of dartboard axes as a fraction of the figure, by default 0.1 + dartboard_axes_height : float, optional + height of dartboard axes as a fraction of the figure, by default 0.2 + verbose : bool + whether or not to display verbose output + """ + if verbose: + print("Getting masks...") + masks = generate_dartboard_masks( + vertex_data.subject, + center_roi, + anchor_rois, + n_angles=n_angles, + n_eccentricities=n_eccentricities, + max_radii=max_radii, + cache_dir=cache_dir,) + # Manage inputs + if isinstance(vertex_data, cx.Vertex): + if vmin is None: + vmin = vertex_data.vmin + if vmax is None: + vmax = vertex_data.vmax + elif isinstance(vertex_data, cx.Vertex2D): + if vmin is None: + vmin = (vertex_data.vmin, vertex_data.vmin2) + if vmax is None: + vmax = (vertex_data.vmax, vertex_data.vmax2) + else: + raise NotImplementedError("No VertexRGB yet!") + if cmap is None: + cmap = vertex_data.cmap + #print(vmin, vmax) + if flatmap_line_style is None: + flatmap_line_style = '-' + if quickflat_kw is None: + quickflat_kw = {} + if roi_border_kw is None: + roi_border_kw = {} + if anchor_rois is None: + anchor_rois = [] + if outline_kw is None: + outline_kw = {} + # Get masked data + data_lh, data_rh = get_dartboard_data(vertex_data, masks, + n_angles=n_angles, n_eccentricities=n_eccentricities, + fn=fn) + # Vertex flatmap plot + fig, ax = plt.subplots(figsize=figsize) + _ = cx.quickflat.make_figure( + vertex_data, with_curvature=True, fig=ax, **quickflat_kw) + + # Augment plot + if show_dartboard_grid: + # Grid fill for dartboard area + vx_grid = generate_dartboard_vertex_object( + masks, vertex_data.subject, type='grid', bg_value=np.nan) + img_grid, extent = cx.quickflat.composite.make_flatmap_image(vx_grid) + ax.imshow(img_grid, extent=extent, + alpha=dartboard_display_alpha, cmap='gray') + if show_dartboard_edge: + # Solid fill for dartboard area + vx_fill = generate_dartboard_vertex_object( + masks, vertex_data.subject, type='solid', bg_value=0) + img_fill, extent = cx.quickflat.composite.make_flatmap_image(vx_fill) + xt = np.linspace(extent[0], extent[1], img_fill.shape[1]) + yt = np.linspace(extent[3], extent[2], img_fill.shape[0]) + xg, yg = np.meshgrid(xt, yt) + ax.contour(xg, yg, img_fill, [1], + linewidths=[flatmap_line_linewidth], + colors=[flatmap_line_color], + zorder=10, # (in front) + ) + if show_anchor_lines: + # Draw lines from center of dartboard to each anchor ROI center + if not isinstance(flatmap_line_style, (list, tuple)): + flatmap_line_style = [flatmap_line_style] * len(anchor_rois) + pts, _ = cx.db.get_surf(vertex_data.subject, + 'flat', merge=True, nudge=True) + roi_centers = {} + for roi in [center_roi] + anchor_rois: + vertex_indices = get_roi_centroids(vertex_data.subject, roi, + surf_type=surf_type, + cache_dir=cache_dir, + verbose=verbose) + roi_centers[roi] = np.array([pts[vertex_indices[0]][:2], + pts[vertex_indices[1]][:2]]) + if roi == center_roi: + center_roi_vertex = vertex_indices + for lr in [0, 1]: # left, right + for pt, ls in zip(anchor_rois, flatmap_line_style): + x = [roi_centers[center_roi][lr, 0], + roi_centers[pt][lr, 0]] + y = [roi_centers[center_roi][lr, 1], + roi_centers[pt][lr, 1]] + ax.plot(x, y, lw=flatmap_line_linewidth, + ls=ls, color=flatmap_line_color) + + # These are in unitless percentages of the figure size. (0,0 is bottom left) + left, bottom, width, height = [0.5 - dartboard_axes_dist_from_midline - dartboard_axes_width, + dartboard_axes_bottom, + dartboard_axes_width, + dartboard_axes_height] + ax_lhem = fig.add_axes([left, bottom, width, height]) + # lh_h = dartboard(data_lh, + # cmap=vertex_data.cmap, + # ax=ax_lhem, + # theta_direction='counter_clockwise',#1, + # vmin=vmin, + # vmax=vmax, + # max_radius=max_radii[0], + # **kwargs) + left, bottom, width, height = [0.5 + dartboard_axes_dist_from_midline, + dartboard_axes_bottom, + dartboard_axes_width, + dartboard_axes_height] + ax_rhem = fig.add_axes([left, bottom, width, height]) + # rh_h = dartboard(data_rh, + # cmap=vertex_data.cmap, + # ax=ax_rhem, + # theta_direction='clockwise', #-1, + # vmin=vmin, + # vmax=vmax, + # max_radius=max_radii[1], + # **kwargs, + # ) + print(cmap) + dartboard_pair((data_lh, data_rh), + vertex_data.subject, + center_roi, + anchor_rois, + display_rois=display_rois, + n_angles=n_angles, + n_eccentricities=n_eccentricities, + max_radii=max_radii, + fn=fn, + surf_type=surf_type, + axs=np.array([ax_lhem, ax_rhem]), + vmin=vmin, + vmax=vmax, + cmap=cmap, + outline_kw=outline_kw, + cache_dir=cache_dir, + verbose=verbose, + **kwargs, + ### + ) + if show_anchor_lines: + # TEMP: Will need changing if anchors are not meant to specify + # 90 degree ticks around dartboard. For that, should be something + # like defining regularly spaced angles: + # angles = np.linspace(0, 2 * np.pi, n_anchor_points) + # ... then computing max_radii[blah] * sin and cos of each angle for X and Y + # ... while following same line order as above. + ax_rhem.plot([-max_radii[1], max_radii[1]], [0, 0], '--', color=flatmap_line_color, + lw=flatmap_line_linewidth) + ax_rhem.plot([0, 0], [-max_radii[1], max_radii[1]], '-', color=flatmap_line_color, + lw=flatmap_line_linewidth) + ax_rhem.axis([-max_radii[1], max_radii[1], - + max_radii[1], max_radii[1]]) + + ax_lhem.plot([-max_radii[0], max_radii[0]], [0, 0], '--', color=flatmap_line_color, + lw=flatmap_line_linewidth) + ax_lhem.plot([0, 0], [-max_radii[0], max_radii[0]], '-', color=flatmap_line_color, + lw=flatmap_line_linewidth) + ax_lhem.axis([-max_radii[0], max_radii[0], - + max_radii[0], max_radii[0]]) + + return fig From e7d722d24547466a3288a5a665b0595f564e4546 Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Wed, 25 Oct 2023 09:02:30 -0700 Subject: [PATCH 06/15] ENH: moved error to make show_dartboard_pair more usable, misc docstring updates --- cortex/dartboards.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index ead53be3..6c009872 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -605,7 +605,7 @@ def show_dartboard(data, show_grid=True, grid_linewidth=0.5, grid_linecolor='lightgray'): - """Given values masked by angle and eccentricity, shows them as a 'dartboard'-style visualization. + """Given values masked by angle and eccentricity, shows them as a radial grid ('dartboard'-style visualization). Parameters ---------- @@ -1002,6 +1002,15 @@ def show_dartboard_pair(dartboard_data, ): """Make two dartboard plots (one for each hemisphere) of dartboard data + Includes many options for plotting, including for plotting of + ROI outlines over dartboard plot. There are several choices to + make here, since the ROI outlines must be warped appropriately. + + For lower-level functions, see: + `show_dartboard` # + `_get_interpolated_outlines` # get roi outlines + `show_outlines_on_dartboard` # plot roi outlines + Parameters ---------- dartboard_data : array or pycortex Vertex object @@ -1073,10 +1082,12 @@ def show_dartboard_pair(dartboard_data, _description_ """ """""" if subject is None: - if not isinstance(dartboard_data, Vertex): - raise ValueError("Must provide either `subject` argument or a Vertex object as `dartboard_data`") - subject = dartboard_data.subject - # WORKING HERE. need to re-figure-out what data_to_plot is. + if isinstance(dartboard_data, Vertex): + subject = dartboard_data.subject + else: + if rois is not None: + raise ValueError("Must provide either `subject` argument or a Vertex object as `dartboard_data` if you wish to display rois") + if isinstance(dartboard_data, Vertex): _, data0 = get_dartboard_data(dartboard_data, **dartboard_spec) else: @@ -1088,9 +1099,7 @@ def show_dartboard_pair(dartboard_data, _, data1 = get_dartboard_data(dartboard_data2, **dartboard_spec) else: data1 = dartboard_data2.copy() - - data_to_plot = (dartboard_data, dartboard_data2) - + if axes is None: fig, axes = plt.subplots(1, 2, figsize=figsize) # Loop over hemispheres @@ -1099,10 +1108,6 @@ def show_dartboard_pair(dartboard_data, for hemi_index, data in enumerate(zip(data0, data1)): axis = axes[hemi_index] d0, d1 = data - # if not isinstance(d0, np.ndarray): - # d0 = get_dartboard_data(d0, **dartboard_spec) - # if d1 is not None: - # d1 = get_dartboard_data(d1, **dartboard_spec) _ = show_dartboard(d0, data2=d1, axis=axis, theta_direction=directions[hemi_index], From d674d43a8d05c06593f433bcfa67e2be1b4b111e Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Wed, 25 Oct 2023 13:21:56 -0700 Subject: [PATCH 07/15] ENH: working code, needs some documentation improvements and possible exposing of other functions, but almost ready. --- cortex/__init__.py | 1 + cortex/dartboards.py | 193 ++++++++++++++++++++++--------------------- 2 files changed, 102 insertions(+), 92 deletions(-) diff --git a/cortex/__init__.py b/cortex/__init__.py index 9ed8226f..b88cd5d8 100644 --- a/cortex/__init__.py +++ b/cortex/__init__.py @@ -8,6 +8,7 @@ from cortex.volume import mosaic, unmask import cortex.export from cortex.version import __version__, __full_version__ +from cortex import dartboards try: from cortex import formats diff --git a/cortex/dartboards.py b/cortex/dartboards.py index 6c009872..e93ae119 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -12,7 +12,7 @@ from . import polyutils from . import quickflat -from .dataset.views import Vertex +from .dataset.views import Vertex, Vertex2D from scipy import interpolate from scipy.spatial import ConvexHull from .utils import get_cmap @@ -329,9 +329,9 @@ def _angles_from_vertex(overlay, start_idx, target_idxs=None, period=2*np.pi, th overlay : SVGOverlay object svg overlay object for a subject start_idx : int - idk + index for center vertex from which to compute angles target_idxs : list, optional - idk, by default None + index for anchor vertex to which to compute angles, by default None period : scalar, optional max angle, by default 2*np.pi theta_direction : str, optional @@ -465,7 +465,7 @@ def _get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, com if return_indices: nearest_vertices.append(roi_verts[hemi][nearest_vertex]) else: - nearest_vertices.append(hemi_roi_coords[:, nearest_vertex]) + nearest_vertices.append(hemi_roi_coords[nearest_vertex, :]) return nearest_vertices @@ -603,6 +603,7 @@ def show_dartboard(data, vmax2=None, theta_direction=-1, show_grid=True, + max_radius=None, grid_linewidth=0.5, grid_linecolor='lightgray'): """Given values masked by angle and eccentricity, shows them as a radial grid ('dartboard'-style visualization). @@ -641,7 +642,8 @@ def show_dartboard(data, plt.Axes Matplotlib axis in which data is plotted. """ - max_radius = 1 + if max_radius is None: + max_radius = 1 data = np.array(data).astype(np.float) if isinstance(data2, np.ndarray): data = np.stack([data, data2], axis=0) @@ -1105,7 +1107,11 @@ def show_dartboard_pair(dartboard_data, # Loop over hemispheres # Hemisphere index goes (0, 1) = (left, right) directions = [-1, 1] + max_radii = dartboard_spec['max_radii'] + if not isinstance(max_radii, (list, tuple)): + max_radii = [max_radii, max_radii] for hemi_index, data in enumerate(zip(data0, data1)): + max_radius = max_radii[hemi_index] axis = axes[hemi_index] d0, d1 = data _ = show_dartboard(d0, data2=d1, @@ -1116,6 +1122,7 @@ def show_dartboard_pair(dartboard_data, show_grid=show_grid, grid_linewidth=grid_linewidth, grid_linecolor=grid_linecolor, + max_radius=max_radius, vmin=vmin, vmax=vmax, vmin2=vmin2, @@ -1435,14 +1442,8 @@ def _get_dartboard_str(**dartboard_spec): fname = '_center{center_roi}-anchors_{anchor_str}-{n_angles}ang-{n_eccentricities}ecc-{rad_str}' return fname.format(**fmt) -def _compute_centroids_angles_from_spec(svg, center_roi, anchors, verbose=False, **kwargs): - """kwargs catches extra inputs from dartboard_spec. - Perhaps not the most principled. - functional. - """ - # Compute centroids of each ROI - t0 = time.time() - # Compute centroids + +def _get_anchor_points(svg, center_roi, anchors, return_indices=True): centroids = {} for j, anchor in enumerate([center_roi] + anchors): if isinstance(anchor, tuple): @@ -1452,15 +1453,30 @@ def _compute_centroids_angles_from_spec(svg, center_roi, anchors, verbose=False, if j == 0: center_name, center_type = anchor_name, anchor_type if anchor_type == 'nearest': - centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name) + centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name, return_indices=return_indices) elif anchor_type == 'centroid': - centroids[anchor_name] = get_roi_centroids(svg, anchor_name) + centroids[anchor_name] = get_roi_centroids(svg, anchor_name, return_indices=return_indices) else: - raise ValueError("unknown anchor type specified: %s\n(Must be 'nearest' or 'centroid')"%(anchor_type)) + raise ValueError("unknown anchor type specified: %s\n(Must be 'nearest' or 'centroid')"%(anchor_type)) + return centroids + + +def _compute_centroids_angles_from_spec(svg, center_roi, anchors, verbose=False, **kwargs): + """kwargs catches extra inputs from dartboard_spec. + Perhaps not the most principled. + functional. + """ + # Compute centroids of each ROI + t0 = time.time() + # Compute centroids + centroids = _get_anchor_points(svg, center_roi, anchors) t1 = time.time() if verbose: print('Time to get centroids:', t1 - t0) - + if isinstance(center_roi, tuple) and (len(center_roi) == 2): + center_name, _ = center_roi + else: + center_name = center_roi # Compute angles from center ROI to each of the anchors anchor_angles_dict = {} for anchor in anchors: @@ -1562,6 +1578,10 @@ def get_dartboard_data(vertex_obj, # Compute the masks, based on specified eccentricity bins, angle bins, and the previously-computed variables t2 = time.time() + # NOTE: magic number here, not great. 4 sub-bins only works specifically + # for 4 anchors and 16 total anglular bins. This could be re-computed with + # an assumption of even spacing, or fundamentally changed by specifying + # the desired angles of the anchor points masks = compute_eccentricity_angle_masks( svg, centroids[center_name], eccentricities=eccentricities, @@ -1589,12 +1609,13 @@ def dartboard_on_flatmap(vertex_data, cmap=None, # Dartboard args center_roi=None, - anchor_rois=None, - display_rois=None, + anchors=None, + rois=None, n_angles=16, n_eccentricities=8, max_radii=(50, 50), - surf_type='inflated', + eccentricities=None, + #surf_type='inflated', # Was for choosing which # Plotting args figsize=(12, 6), dartboard_axes_dist_from_midline=0.15, @@ -1605,7 +1626,7 @@ def dartboard_on_flatmap(vertex_data, quickflat_kw=None, flatmap_line_linewidth=1.5, flatmap_line_color='y', - flatmap_line_style=None, + flatmap_line_style=('--','-','--','-'), show_anchor_lines=None, show_dartboard_grid=True, show_dartboard_edge=True, @@ -1614,10 +1635,9 @@ def dartboard_on_flatmap(vertex_data, n_roi_border_points=64, roi_outline_smooth_factor=5, # every 5th point kept, smoothed with cubic spline roi_border_kw=None, - cache_dir=None, verbose=False, outline_kw=None, - **kwargs): + ): """Make a flatmap with overlaid dartboard plots Parameters @@ -1652,21 +1672,31 @@ def dartboard_on_flatmap(vertex_data, """ if verbose: print("Getting masks...") - masks = generate_dartboard_masks( - vertex_data.subject, - center_roi, - anchor_rois, - n_angles=n_angles, - n_eccentricities=n_eccentricities, - max_radii=max_radii, - cache_dir=cache_dir,) + # Allow manually specified eccentricities to override linear spacing + if eccentricities is None: + eccentricities = [np.linspace( + 0, mr, n_eccentricities + 1) for mr in max_radii] + dartboard_spec = dict(center_roi=center_roi, + anchors=anchors, + n_angles=n_angles, + n_eccentricities=n_eccentricities, + eccentricities=eccentricities, + max_radii=max_radii,) + # masks = generate_dartboard_masks( + # vertex_data.subject, + # center_roi, + # anchor_rois, + # n_angles=n_angles, + # n_eccentricities=n_eccentricities, + # max_radii=max_radii, + # cache_dir=cache_dir,) # Manage inputs - if isinstance(vertex_data, cx.Vertex): + if isinstance(vertex_data, Vertex): if vmin is None: vmin = vertex_data.vmin if vmax is None: vmax = vertex_data.vmax - elif isinstance(vertex_data, cx.Vertex2D): + elif isinstance(vertex_data, Vertex2D): if vmin is None: vmin = (vertex_data.vmin, vertex_data.vmin2) if vmax is None: @@ -1675,39 +1705,37 @@ def dartboard_on_flatmap(vertex_data, raise NotImplementedError("No VertexRGB yet!") if cmap is None: cmap = vertex_data.cmap - #print(vmin, vmax) if flatmap_line_style is None: flatmap_line_style = '-' if quickflat_kw is None: - quickflat_kw = {} + quickflat_kw = {'with_curvature' : True, } if roi_border_kw is None: roi_border_kw = {} - if anchor_rois is None: - anchor_rois = [] if outline_kw is None: outline_kw = {} - # Get masked data - data_lh, data_rh = get_dartboard_data(vertex_data, masks, - n_angles=n_angles, n_eccentricities=n_eccentricities, - fn=fn) + # # Get masked data + # data_lh, data_rh = get_dartboard_data(vertex_data, masks, + # n_angles=n_angles, n_eccentricities=n_eccentricities, + # fn=fn) + masks, to_plot = get_dartboard_data(vertex_data, **dartboard_spec, mean_func=fn) # Vertex flatmap plot fig, ax = plt.subplots(figsize=figsize) - _ = cx.quickflat.make_figure( - vertex_data, with_curvature=True, fig=ax, **quickflat_kw) + _ = quickflat.make_figure( + vertex_data, fig=ax, **quickflat_kw) # Augment plot if show_dartboard_grid: # Grid fill for dartboard area vx_grid = generate_dartboard_vertex_object( masks, vertex_data.subject, type='grid', bg_value=np.nan) - img_grid, extent = cx.quickflat.composite.make_flatmap_image(vx_grid) + img_grid, extent = quickflat.composite.make_flatmap_image(vx_grid) ax.imshow(img_grid, extent=extent, alpha=dartboard_display_alpha, cmap='gray') if show_dartboard_edge: # Solid fill for dartboard area vx_fill = generate_dartboard_vertex_object( masks, vertex_data.subject, type='solid', bg_value=0) - img_fill, extent = cx.quickflat.composite.make_flatmap_image(vx_fill) + img_fill, extent = quickflat.composite.make_flatmap_image(vx_fill) xt = np.linspace(extent[0], extent[1], img_fill.shape[1]) yt = np.linspace(extent[3], extent[2], img_fill.shape[0]) xg, yg = np.meshgrid(xt, yt) @@ -1717,29 +1745,33 @@ def dartboard_on_flatmap(vertex_data, zorder=10, # (in front) ) if show_anchor_lines: + overlay = db.get_overlay(vertex_data.subject) # Draw lines from center of dartboard to each anchor ROI center if not isinstance(flatmap_line_style, (list, tuple)): - flatmap_line_style = [flatmap_line_style] * len(anchor_rois) - pts, _ = cx.db.get_surf(vertex_data.subject, + flatmap_line_style = [flatmap_line_style] * len(anchors) + pts, _ = db.get_surf(vertex_data.subject, 'flat', merge=True, nudge=True) roi_centers = {} - for roi in [center_roi] + anchor_rois: - vertex_indices = get_roi_centroids(vertex_data.subject, roi, - surf_type=surf_type, - cache_dir=cache_dir, - verbose=verbose) - roi_centers[roi] = np.array([pts[vertex_indices[0]][:2], - pts[vertex_indices[1]][:2]]) - if roi == center_roi: - center_roi_vertex = vertex_indices + roi_center_indices = _get_anchor_points(overlay, center_roi, anchors) + for roi, c in roi_center_indices.items(): + roi_centers[roi] = np.array([pts[c[0]][:2], + pts[c[1]][:2]]) + center_pt = roi_centers.pop(center_roi) for lr in [0, 1]: # left, right - for pt, ls in zip(anchor_rois, flatmap_line_style): - x = [roi_centers[center_roi][lr, 0], + for pt, ls in zip(roi_centers.keys(), flatmap_line_style): + x = [center_pt[lr, 0], roi_centers[pt][lr, 0]] - y = [roi_centers[center_roi][lr, 1], + y = [center_pt[lr, 1], roi_centers[pt][lr, 1]] ax.plot(x, y, lw=flatmap_line_linewidth, ls=ls, color=flatmap_line_color) + # Replace below with...? + # overlaid_axis = fig.add_axes( + # mtrans.Bbox([ + # [position_x-scale/2, position_y-scale/2], + # [position_x+scale/2, position_y+scale/2]]), + # frameon=False + # ) # These are in unitless percentages of the figure size. (0,0 is bottom left) left, bottom, width, height = [0.5 - dartboard_axes_dist_from_midline - dartboard_axes_width, @@ -1747,51 +1779,28 @@ def dartboard_on_flatmap(vertex_data, dartboard_axes_width, dartboard_axes_height] ax_lhem = fig.add_axes([left, bottom, width, height]) - # lh_h = dartboard(data_lh, - # cmap=vertex_data.cmap, - # ax=ax_lhem, - # theta_direction='counter_clockwise',#1, - # vmin=vmin, - # vmax=vmax, - # max_radius=max_radii[0], - # **kwargs) left, bottom, width, height = [0.5 + dartboard_axes_dist_from_midline, dartboard_axes_bottom, dartboard_axes_width, dartboard_axes_height] ax_rhem = fig.add_axes([left, bottom, width, height]) - # rh_h = dartboard(data_rh, - # cmap=vertex_data.cmap, - # ax=ax_rhem, - # theta_direction='clockwise', #-1, - # vmin=vmin, - # vmax=vmax, - # max_radius=max_radii[1], - # **kwargs, - # ) - print(cmap) - dartboard_pair((data_lh, data_rh), - vertex_data.subject, - center_roi, - anchor_rois, - display_rois=display_rois, - n_angles=n_angles, - n_eccentricities=n_eccentricities, - max_radii=max_radii, - fn=fn, - surf_type=surf_type, - axs=np.array([ax_lhem, ax_rhem]), + + show_dartboard_pair(vertex_data, + **dartboard_spec, + rois=rois, + mean_func=fn, + axes=np.array([ax_lhem, ax_rhem]), vmin=vmin, vmax=vmax, cmap=cmap, - outline_kw=outline_kw, - cache_dir=cache_dir, + #outline_kw=outline_kw, verbose=verbose, - **kwargs, + #**kwargs, ### ) if show_anchor_lines: - # TEMP: Will need changing if anchors are not meant to specify + # Show vertical and horizontal lines on dartboard plots + # NOTE: Will need changing if anchors are not meant to specify # 90 degree ticks around dartboard. For that, should be something # like defining regularly spaced angles: # angles = np.linspace(0, 2 * np.pi, n_anchor_points) From a025e38881aea2ea1b6825719d608d5881307480 Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Wed, 1 Nov 2023 14:28:24 -0700 Subject: [PATCH 08/15] ENH: documentation added, examples modified --- cortex/dartboards.py | 64 ++++++++++++------- examples/dartboards/dartboard_basic.py | 32 +++++++--- .../dartboards/dartboard_flatmap_overlay.py | 33 +++++++++- .../dartboard_low_level_computations_plots.py | 31 +++++++++ requirements.txt | 1 + 5 files changed, 127 insertions(+), 34 deletions(-) create mode 100644 examples/dartboards/dartboard_low_level_computations_plots.py diff --git a/cortex/dartboards.py b/cortex/dartboards.py index e93ae119..fdfc4ed8 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -1,7 +1,6 @@ from matplotlib import pyplot as plt from matplotlib import transforms as mtrans from matplotlib.colors import Normalize -from sklearn.metrics import pairwise_distances import matplotlib import numpy as np import time @@ -14,7 +13,7 @@ from .dataset.views import Vertex, Vertex2D from scipy import interpolate -from scipy.spatial import ConvexHull +from scipy.spatial import ConvexHull, distance from .utils import get_cmap try: @@ -23,6 +22,12 @@ CACHE_DIR = options.config.get('basic','filestore') ## Utils + +def _pairwise_distances(array, metric='euclidean'): + """Simple wrapper for scipy distance function""" + d = distance.pdist(array, metric=metric) + return distance.squareform(d) + def flatten_list_of_lists(list_of_lists): flat_list = [item for sublist in list_of_lists for item in sublist] return flat_list @@ -444,7 +449,7 @@ def get_closest_vertices(overlay, coords, hemi='both', distance_metric='euclidea overlay_coords = overlay_coords[n_vertices_left:] else: raise ValueError("hemi must be 0, 1, 'left', 'right' or 'both'.") - min_idxs = np.nanargmin(pairwise_distances(coords,overlay_coords,metric=distance_metric), axis=1) + min_idxs = np.nanargmin(_pairwise_distances(coords,overlay_coords,metric=distance_metric), axis=1) if hemi in [1,'r','right']: min_idxs += n_vertices_left return min_idxs @@ -459,7 +464,7 @@ def _get_closest_vertex_to_roi(overlay, roi, comparison_roi, roi_full=False, com for hemi in range(2): hemi_roi_coords = all_coords[roi_verts[hemi]] hemi_comparison_roi_coords = all_coords[comparison_roi_verts[hemi]] - distances = pairwise_distances( + distances = _pairwise_distances( hemi_roi_coords, hemi_comparison_roi_coords, metric=distance_metric) nearest_vertex = np.nanargmin(np.nanmin(distances, axis=1)) if return_indices: @@ -1505,10 +1510,10 @@ def get_dartboard_data(vertex_obj, Parameters ---------- - vertex_obj : _type_ - _description_ - center_roi : _type_ - _description_ + vertex_obj : cortex.Vertex + vertex object containing data to be averaged over the defined radial bins + center_roi : str or tuple + name of roi (defined in overlays.svg for this subject) anchors : list list of strings or tuples specifying anchor points. If strings are provided, code assumes strings to be names for regions of interest defined for this @@ -1528,6 +1533,9 @@ def get_dartboard_data(vertex_obj, 3 <- center -> 1 \/ 4 + max_radii : int or tuple + maximum radius to extend from center ROI (in mm in fiducial space). Separate + values for each hemisphere can be provided as a tuple (Lhem, Rhem). mean_func : function, optional The function used to average all vertices in a given bin, by default np.nanmean Consider what is being averaged to make an appropriate choice here (e.g. @@ -1625,12 +1633,11 @@ def dartboard_on_flatmap(vertex_data, dartboard_display_alpha=0.2, quickflat_kw=None, flatmap_line_linewidth=1.5, - flatmap_line_color='y', + flatmap_line_color='c', flatmap_line_style=('--','-','--','-'), show_anchor_lines=None, show_dartboard_grid=True, show_dartboard_edge=True, - #show_roi_outlines=True, # ROI outline parameters n_roi_border_points=64, roi_outline_smooth_factor=5, # every 5th point kept, smoothed with cubic spline @@ -1651,12 +1658,31 @@ def dartboard_on_flatmap(vertex_data, function to use to collapse over vertices within each mask, by default np.nanmean center_roi : str center ROI + anchors : list + list of strings or tuples specifying anchor points. If strings are provided, + code assumes strings to be names for regions of interest defined for this + pycortex subject (e.g. 'FFA'). Anchor points are computed as centroids of + the named ROIs. If tuples are provided, the first element of the tuple is the + name of the ROI, and the second element must be 'centroid' or 'nearest', + indicating whether the anchor point shoudl be the centroi of the named ROI, or + the nearest point in or on the named region to the center_roi. Sulci or other + labeled anatomical markers in other layers of overlays.svg for the subject may + be used, e.g. ('STS', 'nearest') for the nearest point to the center_roi that + falls on STS. + + Important: anchors must be specified counter clockwise from the RIGHT with + respect to the RIGHT hemisphere: + 2 + /\ + 3 <- center -> 1 + \/ + 4 n_eccentricities : int, optional number of radial bins, by default 5; should be consistent with number of radial bins in `masks`, if `masks` is provided n_angles : int, optional - number of polar angle bins, by default 8; should be consistent with number of - angular bins in `masks`, if `masks` is provided + number of polar angle bins, by default 16; will currently break with anything but + 16 and 4 anchors... dartboard_axes_dist_from_midline : float, optional fraction of the figure by which dartboard axes are displaced from the figure midline, by default 0.15 @@ -1682,14 +1708,7 @@ def dartboard_on_flatmap(vertex_data, n_eccentricities=n_eccentricities, eccentricities=eccentricities, max_radii=max_radii,) - # masks = generate_dartboard_masks( - # vertex_data.subject, - # center_roi, - # anchor_rois, - # n_angles=n_angles, - # n_eccentricities=n_eccentricities, - # max_radii=max_radii, - # cache_dir=cache_dir,) + # Manage inputs if isinstance(vertex_data, Vertex): if vmin is None: @@ -1713,10 +1732,7 @@ def dartboard_on_flatmap(vertex_data, roi_border_kw = {} if outline_kw is None: outline_kw = {} - # # Get masked data - # data_lh, data_rh = get_dartboard_data(vertex_data, masks, - # n_angles=n_angles, n_eccentricities=n_eccentricities, - # fn=fn) + # Load dartboard masks & masked data masks, to_plot = get_dartboard_data(vertex_data, **dartboard_spec, mean_func=fn) # Vertex flatmap plot fig, ax = plt.subplots(figsize=figsize) diff --git a/examples/dartboards/dartboard_basic.py b/examples/dartboards/dartboard_basic.py index 79bc59a7..41473b96 100644 --- a/examples/dartboards/dartboard_basic.py +++ b/examples/dartboards/dartboard_basic.py @@ -1,4 +1,8 @@ -# Basic dartboard pair +"""Minimal demo of dartboard (radial grid bin) plot. + +Note that this type of radial grid averaging is best +suited for regions for which there is some hypothesis about +organization centered on a functionally or anatomically defined point in the cortex. This demo uses """ import cortex subject = 'S1' @@ -10,12 +14,22 @@ # vx = vol.map() vx = cortex.db.get_surfinfo(subject, type='curvature') -center = 'M1H' - -anchors = [('S1H', 'centroid'), - ('M1F', 'centroid'), - ('FEF', 'centroid'), - ('M1M', 'centroid'), - ] +# Create a dartboard centered on M1H (primary motor cortex hand representation), +# anchored to FEF (frontal eye fields above) anterior, +# M1F (primary motor cortex foot representation) superior, +# S1H (sensorymotor cortex hand representation) posterior, +# and M1M (primary motor cortex mouth representation) inferior. +# This should show the central sulcus going through the radial grid. +dartboard_spec = dict( + center = 'M1H', + anchors = [('FEF', 'centroid'), + ('M1F', 'centroid'), + ('S1H', 'centroid'), + ('M1M', 'centroid'), + ], + n_angles=16, + n_eccentricities=8, + max_radii=(50,50), +) -out_m1h = cortex.dartboards.get_dartboard_data(vx, center, anchors) +axs = cortex.dartboards.show_dartboard_pair(vx, **dartboard_spec) diff --git a/examples/dartboards/dartboard_flatmap_overlay.py b/examples/dartboards/dartboard_flatmap_overlay.py index 4b9ddc97..99755cff 100644 --- a/examples/dartboards/dartboard_flatmap_overlay.py +++ b/examples/dartboards/dartboard_flatmap_overlay.py @@ -1,3 +1,34 @@ # Dartboard overlay on flatmap -# WIP \ No newline at end of file +# Basic dartboard pair +import cortex + +subject = 'S1' +# Get curvature information for this subject +# Note that this is just an easy example of vertex data; +# Vertex data can always be created by creating a Volume +# and mapping it to vertices, like this: +# vol = cx.Volume() +# vx = vol.map() +vx = cortex.db.get_surfinfo(subject, type='curvature') + +# Create a dartboard centered on M1H (primary motor cortex hand representation), +# anchored to FEF (frontal eye fields above) anterior, +# M1F (primary motor cortex foot representation) superior, +# S1H (sensorymotor cortex hand representation) posterior, +# and M1M (primary motor cortex mouth representation) inferior. +# This should show the central sulcus going through the radial grid. +dartboard_spec = dict( + center='M1H', + anchors=[('FEF', 'centroid'), + ('M1F', 'centroid'), + ('S1H', 'centroid'), + ('M1M', 'centroid'), + ], + n_angles=16, + n_eccentricities=8, + max_radii=(50, 50), +) + +# Show dartboard plot overlaid on flatmap from which it was derived +axs = cortex.dartboards.dartboard_on_flatmap(vx, **dartboard_spec) diff --git a/examples/dartboards/dartboard_low_level_computations_plots.py b/examples/dartboards/dartboard_low_level_computations_plots.py new file mode 100644 index 00000000..8214038d --- /dev/null +++ b/examples/dartboards/dartboard_low_level_computations_plots.py @@ -0,0 +1,31 @@ +# This demo shows the steps of computation of dartboard plots (radial bin averaging and breakdown of plots with low-level functions) +import cortex + +subject = 'S1' +# Get curvature information for this subject +# Note that this is just an easy example of vertex data; +# Vertex data can always be created by creating a Volume +# and mapping it to vertices, like this: +# vol = cx.Volume() +# vx = vol.map() +vx = cortex.db.get_surfinfo(subject, type='curvature') + +# Create a dartboard centered on M1H (primary motor cortex hand representation), +# anchored to FEF (frontal eye fields above) anterior, +# M1F (primary motor cortex foot representation) superior, +# S1H (sensorymotor cortex hand representation) posterior, +# and M1M (primary motor cortex mouth representation) inferior. +# This should show the central sulcus going through the radial grid. +dartboard_spec = dict( + center='M1H', + anchors=[('FEF', 'centroid'), + ('M1F', 'centroid'), + ('S1H', 'centroid'), + ('M1M', 'centroid'), + ], + n_angles=16, + n_eccentricities=8, + max_radii=(50, 50), +) + +masks, data = cortex.dartboards.get_dartboard_data(vx, **dartboard_spec) diff --git a/requirements.txt b/requirements.txt index b30db96c..acfa10d6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -14,3 +14,4 @@ pillow nibabel>=2.1 networkx>=2.1 imageio +wget From 650e0bc30d0487f4d1bca408db70efe71108d52d Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Fri, 10 Nov 2023 11:09:03 -0800 Subject: [PATCH 09/15] ENH: wip changes to demos, pairwise_distances function switch to scipy --- cortex/dartboards.py | 10 +++++----- .../dartboard_low_level_computations_plots.py | 9 +++++++++ 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index fdfc4ed8..40e53cf6 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -22,11 +22,11 @@ CACHE_DIR = options.config.get('basic','filestore') ## Utils - -def _pairwise_distances(array, metric='euclidean'): - """Simple wrapper for scipy distance function""" - d = distance.pdist(array, metric=metric) - return distance.squareform(d) +from sklearn.metrics import pairwise_distances as _pairwise_distances +# def _pairwise_distances(array, metric='euclidean'): +# """Simple wrapper for scipy distance function""" +# d = distance.pdist(array, metric=metric) +# return distance.squareform(d) def flatten_list_of_lists(list_of_lists): flat_list = [item for sublist in list_of_lists for item in sublist] diff --git a/examples/dartboards/dartboard_low_level_computations_plots.py b/examples/dartboards/dartboard_low_level_computations_plots.py index 8214038d..83682b75 100644 --- a/examples/dartboards/dartboard_low_level_computations_plots.py +++ b/examples/dartboards/dartboard_low_level_computations_plots.py @@ -1,5 +1,6 @@ # This demo shows the steps of computation of dartboard plots (radial bin averaging and breakdown of plots with low-level functions) import cortex +import numpy as np subject = 'S1' # Get curvature information for this subject @@ -28,4 +29,12 @@ max_radii=(50, 50), ) +# define vertex-averaging function +def mean_nonan(x, axis=None, threshold=0.8): + if np.mean(np.isnan(x), axis=axis) > threshold: + return np.nan + else: + return np.nanmean(x) + +# Get vertex-wise masks for each dartboard bin (`masks`) and masked data masks, data = cortex.dartboards.get_dartboard_data(vx, **dartboard_spec) From 99065c6829a44763e6bb09b661da8c1cc9651449 Mon Sep 17 00:00:00 2001 From: Arnab Biswas Date: Tue, 14 Nov 2023 12:18:19 -0800 Subject: [PATCH 10/15] ENH: add param cutoff to get_dartboard_data Adding cutoff param to select n_vertices or % vertices below which to set dartboard bin to nan. Functionality already exists in apply_masks function. Just added the param to get_dartboard_data. --- cortex/dartboards.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index 40e53cf6..b0acf297 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -1503,6 +1503,7 @@ def get_dartboard_data(vertex_obj, n_eccentricities=8, eccentricities=None, mean_func=np.nanmean, + cutoff = None, recache=False, verbose=True, ): @@ -1543,6 +1544,11 @@ def get_dartboard_data(vertex_obj, require a custom function) eccentricities : _type_, optional _description_, by default np.linspace(0, max_radii, 8+1) + cutoff : int or float, optional + Cutoff value for the minimum number of vertices for a bin. If an int, the threshold + will be based on the absolute number of vertices. If a float, the threshold will be + based on the percentage of vertices included per bin. Bins with fewer vertices than + the cutoff will return np.nan. Defaults to None. Returns ------- @@ -1605,7 +1611,7 @@ def get_dartboard_data(vertex_obj, print('Time to compute masks:', t3-t2) print(f'Saving dartboard masks to {fpath}') np.save(fpath, masks) - output = apply_masks(vertex_obj.data, masks, mean_func) + output = apply_masks(vertex_obj.data, masks, mean_func, cutoff) return masks, output From a0796f500e1f59089d174ab18cd4d4f39b041166 Mon Sep 17 00:00:00 2001 From: MShinkle <49564869+MShinkle@users.noreply.github.com> Date: Mon, 20 Nov 2023 17:29:16 -0800 Subject: [PATCH 11/15] feat: replace sklearn depend with scipy alternative --- cortex/dartboards.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index b0acf297..2564f2cb 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -21,12 +21,11 @@ except: CACHE_DIR = options.config.get('basic','filestore') -## Utils -from sklearn.metrics import pairwise_distances as _pairwise_distances -# def _pairwise_distances(array, metric='euclidean'): -# """Simple wrapper for scipy distance function""" -# d = distance.pdist(array, metric=metric) -# return distance.squareform(d) +def pairwise_distances(X, Y=None, metric='euclidean', **kwargs): + if Y is None: + Y = X + D = distance.cdist(X, Y, metric=metric, **kwargs) + return D def flatten_list_of_lists(list_of_lists): flat_list = [item for sublist in list_of_lists for item in sublist] From 6b3abfec64467387195cc18f1b9376a8a13b1d36 Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Mon, 4 Dec 2023 10:42:14 -0800 Subject: [PATCH 12/15] FIX: added _ to _pairwise_distances; code should run. --- cortex/dartboards.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index 2564f2cb..8aeb79ef 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -21,7 +21,7 @@ except: CACHE_DIR = options.config.get('basic','filestore') -def pairwise_distances(X, Y=None, metric='euclidean', **kwargs): +def _pairwise_distances(X, Y=None, metric='euclidean', **kwargs): if Y is None: Y = X D = distance.cdist(X, Y, metric=metric, **kwargs) From 96c0bb5bfd95ff1be65c2dc57a2e408777bf6b8d Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Mon, 4 Dec 2023 10:45:10 -0800 Subject: [PATCH 13/15] FIX: fixed spelling error in docstring --- cortex/dartboards.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index 8aeb79ef..1db44670 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -1520,7 +1520,7 @@ def get_dartboard_data(vertex_obj, pycortex subject (e.g. 'FFA'). Anchor points are computed as centroids of the named ROIs. If tuples are provided, the first element of the tuple is the name of the ROI, and the second element must be 'centroid' or 'nearest', - indicating whether the anchor point shoudl be the centroi of the named ROI, or + indicating whether the anchor point should be the centroi of the named ROI, or the nearest point in or on the named region to the center_roi. Sulci or other labeled anatomical markers in other layers of overlays.svg for the subject may be used, e.g. ('STS', 'nearest') for the nearest point to the center_roi that @@ -1669,7 +1669,7 @@ def dartboard_on_flatmap(vertex_data, pycortex subject (e.g. 'FFA'). Anchor points are computed as centroids of the named ROIs. If tuples are provided, the first element of the tuple is the name of the ROI, and the second element must be 'centroid' or 'nearest', - indicating whether the anchor point shoudl be the centroi of the named ROI, or + indicating whether the anchor point should be the centroi of the named ROI, or the nearest point in or on the named region to the center_roi. Sulci or other labeled anatomical markers in other layers of overlays.svg for the subject may be used, e.g. ('STS', 'nearest') for the nearest point to the center_roi that From 1ab0be7785bd0b1e4ed8cfb5bc3fffd47f3ba3d7 Mon Sep 17 00:00:00 2001 From: Mark Lescroart Date: Wed, 17 Jan 2024 12:30:29 -0800 Subject: [PATCH 14/15] FIX: updated roi outline code, WIP --- cortex/dartboards.py | 198 ++++++++++++++---- .../dartboard_low_level_computations_plots.py | 9 + 2 files changed, 171 insertions(+), 36 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index 1db44670..f19d22a3 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -63,6 +63,28 @@ def pol2cart(phi, rho): return x, y +def circ_dist(a, b, degrees=True): + """Compute angle between two angles w/ circular statistics + + Parameters + ---------- + a : scalar or array + angle(s) TO WHICH to compute angular (aka rotational) distance + b : scalar or array + angle(s) FROM WHICH to compute angular (aka rotational) distance + degrees : bool + Whether a and b are in degrees (defaults to True; False means they are in radians) + """ + if degrees: + a = np.radians(a) + b = np.radians(b) + phase = np.e**(1j*a) / np.e**(1j*b) + dist = np.arctan2(phase.imag, phase.real) + if degrees: + dist = np.degrees(dist) + return dist + + def nonzero_coords(coords): return coords[np.any(coords != 0, axis=1)] @@ -737,25 +759,38 @@ def resample_roi_outline( distances, resolution=100, every_n=5, + start_angle=np.pi/2, + fixed_sampling_between=(0, 90, 180, 270), # TO DO even_sampling_over='polar angle', ): """resample roi outline to smooth outline + + Also useful for creating average ROI outlines for different participants + that can be averaged in dartboard space to create an average ROI. - Relies on angles being monotonically increasing + Relies on angles being monotonically increasing. Parameters ---------- angles : _type_ - _description_ - even_sampling_over : str, optional - _description_, by default 'polar angle' + angle to each ROI outline point. This is theta of (rho, theta) + [i.e. (radius, polar angle)] in polar coordinates. + even_sampling_over : str, or None, optional + method by which to sample along roi border, either 'path length', + 'polar angle', or None; by default 'path length'. Sampling over + polar angle from the center of the ROI is a sensible option if the + ROI is approximately circular, without drastic elongation or non-convex + regions along the border. For ROIs with non-convex regions, + sampling over path length (at even distances along the ROI border) + will give more sensible results. Given that the failure case for + sampling resolution : int, optional number of points to sample along the ROI border, by default 100 Returns ------- - _type_ - _description_ + angles, distances + returns angles and distances in polar coordinates. """ delta = 1e-5 # Periodic linear resample (upsample) to assure values from 0 to 2pi @@ -777,10 +812,12 @@ def resample_roi_outline( angles_out = np.linspace(0, 2 * np.pi - delta, resolution) dist_out = smooth(angles_out) return angles_out, dist_out - - if even_sampling_over == 'path length': + elif even_sampling_over == 'path length': ## Resample xy = np.array(pol2cart(angles, distances)).T + vertex_order = sort_roi_border(xy, angles, start_angle=start_angle, max_deg_from_zero=5, verbose=False) + xy = xy[vertex_order] + angles = angles[vertex_order] path_dist = np.cumsum( np.hstack([0, np.linalg.norm(np.diff(xy, axis=0), axis=1)])) max_dist = np.max(path_dist) @@ -810,9 +847,69 @@ def resample_roi_outline( # Convert back to angles angles_out, dist_out = cart2pol(*xy_out.T) return angles_out, dist_out + elif even_sampling_over is None: + return angles, distances + +def sort_roi_border(xy, angles, start_angle=0, start_index=None, max_deg_from_zero=5, verbose=False): + """start_index overrides start_angle""" + dsts = distance.squareform(distance.pdist(xy)) + if start_index is None: + # Choose max distance vertex that is approximately straight up + # from whatever center the angles were computed from. + # Best would be maybe to compute max distance from actual center + # that is the center from which angles were computed. + deg_from_zero = 2 + while deg_from_zero <= max_deg_from_zero: + max_radians_from_up = np.radians(max_deg_from_zero) + close_to_straight_up, = np.nonzero(np.abs(circ_dist(angles, start_angle, degrees=False)) <= max_radians_from_up) + if len(close_to_straight_up) > 0: + # pick farthest away from approx. center that is straight up + dist_from_median = np.linalg.norm(xy - np.median(xy, axis=0), axis=1) + jj = np.argmax(dist_from_median[close_to_straight_up]) + start_index = close_to_straight_up[jj] + break + else: + deg_from_zero += 1 + if start_index is None: + start_index = np.argmin( + np.abs(circ_dist(angles, start_angle, degrees=False))) + vertex_order = [start_index] + ct = 0 + for i_vert in range(len(xy)): + # Special case for first one: go clockwise + if i_vert == 0: + new_verts = np.argsort(dsts[start_index])[:5] + new_verts = np.array([x for x in new_verts if x not in vertex_order]) + angles_from_origin = circ_dist(angles[new_verts], angles[start_index], degrees=False) + is_clockwise = angles_from_origin > 0 + #if ~any(is_clockwise): + # plt.plot(*xy.T) + # plt.plot() + min_clockwise_angle = np.min(angles_from_origin[is_clockwise]) + j = new_verts[angles_from_origin == min_clockwise_angle][0] + vertex_order.append(j) + else: + success = False + n = 3 + while (not success) and (n<15): + new_verts = np.argsort(dsts[vertex_order[-1]])[:n] + new_verts = np.array([x for x in new_verts if x not in vertex_order]) + success = len(new_verts==1) + n += 1 + #if n > 4: + # print(n) + if len(new_verts) > 0: + vertex_order.append(new_verts[-1]) + else: + ct += 1 + if verbose: + print(f"Skipped a vertex ({ct} skipped)") + #1/0 + vertex_order.append(vertex_order[0]) + return np.array(vertex_order) -def _get_interpolated_outlines(overlay, +def get_interpolated_outlines(overlay, outline_roi, geodesic_distances=True, resolution=100, @@ -898,9 +995,6 @@ def _get_interpolated_outlines(overlay, # Convert to polar relative to these centroids relative_phi, relative_rho = cart2pol(x,y) # Even angular interpolation - ## CHANGED HERE: May require re-ordering of angles, may rest on assumption - ## about what angle is first - #relative_phi_interp, relative_rho_interp = interpolate_outlines(relative_phi, relative_rho, resolution=resolution) relative_phi_interp, relative_rho_interp = resample_roi_outline( relative_phi, relative_rho, @@ -913,6 +1007,7 @@ def _get_interpolated_outlines(overlay, # Add back own centroid x_interp += centroid_x y_interp += centroid_y + # Convert back to polar coordinates interpolated_outlines.append(np.stack(cart2pol(x_interp,y_interp), axis=1)) # cache result if verbose: @@ -926,13 +1021,18 @@ def show_outlines_on_dartboard( **plot_kwargs): """Plot outlines of regions of interest on dartboard plots. + Plots all ROIs for one hemisphere to one axis. + Parameters ---------- - outlines : _type_ - _description_ + outlines : array + array of size (n_rois, n_points_along_boundary, coordinates), containing polar + coordinates for outline of each ROI. Polar coordinates are expected to be + (rho, theta), i.e. (radius, polar angle) + axis : matplotlib.axes._subplots.AxesSubplot, optional Matplotlib axis on which to plot. - colors : _type_, optional + colors : matplotlib colorspec, optional List of colors to iterate through when plotting outlines. If None, will iterate through the default pyplot color cycle. Defaults to None. polar : bool, optional @@ -1014,7 +1114,7 @@ def show_dartboard_pair(dartboard_data, For lower-level functions, see: `show_dartboard` # - `_get_interpolated_outlines` # get roi outlines + `get_interpolated_outlines` # get roi outlines `show_outlines_on_dartboard` # plot roi outlines Parameters @@ -1139,7 +1239,7 @@ def show_dartboard_pair(dartboard_data, if not isinstance(roi_color, (list, tuple)): roi_color = [roi_color] * len(rois) for roi in rois: - outlines[roi] = _get_interpolated_outlines(overlay, + outlines[roi] = get_interpolated_outlines(overlay, roi, geodesic_distances=geodesic_distances, resolution=path_resolution, @@ -1150,6 +1250,10 @@ def show_dartboard_pair(dartboard_data, **dartboard_spec) for hemi_index in range(2): hemi_axis = axes[hemi_index] + # This may be somewhat confusing. Seems sensible to make the function plot + # either one ROI at a time or to take a dict as input to obviate the need + # for mapping to an array (and a totally different format than any other + # function provides) here. outlines_to_plot = np.array([v for v in outlines.values()])[:, hemi_index] if isinstance(dartboard_spec['max_radii'], tuple): rmax = dartboard_spec['max_radii'][hemi_index] @@ -1171,21 +1275,27 @@ def draw_mask_outlines( outer_line=True, as_overlay=True, **plot_kwargs): """Given eccentricity-angle masks, fits and draws outlines of bins on the cortex. - Args: - masks np.ndarray: Vertex masks for each bin and hemisphere, such as the output of - get_eccentricity_angle_masks. Shape should be (2, eccentricities, angles, vertices). - axis (matplotlib.axes._subplots.AxesSubplot, optional): Matplotlb axis on which to plot. - eccentricities (bool, optional): Whether to draw outlines for separate eccentricities. - Defaults to True. - angles (bool, optional): Whether to draw outlines for separate angles. - Defaults to True. - outer_line (bool, optional): Whether to draw outline of outermost eccentricity. - Defaults to True. - as_overlay (bool, optional): Whether to create new axis on top of existing axis, to avoid - interfering with previously-plotted data. Defaults to False. - - Returns: - matplotlib.axes._subplots.AxesSubplot: Matplotlib axis in which data is plotted. + Parameters + ---------- + masks : np.ndarray + Vertex masks for each bin and hemisphere, such as the output of get_eccentricity_angle_masks. + Shape should be (2, eccentricities, angles, vertices). + axis : matplotlib.axes._subplots.AxesSubplot, optional + Matplotlib axis on which to plot. + eccentricities : bool, optional + Whether to draw outlines for separate eccentricities. Defaults to True. + angles : bool, optional + Whether to draw outlines for separate angles. Defaults to True. + outer_line : bool, optional + Whether to draw the outline of the outermost eccentricity. Defaults to True. + as_overlay : bool, optional + Whether to create a new axis on top of the existing axis to avoid interfering with + previously-plotted data. Defaults to False. + + Returns + ------- + matplotlib.axes._subplots.AxesSubplot + Matplotlib axis in which data is plotted. """ if axis is None: _, axis = plt.subplots() @@ -1526,8 +1636,9 @@ def get_dartboard_data(vertex_obj, be used, e.g. ('STS', 'nearest') for the nearest point to the center_roi that falls on STS. - Important: anchors must be specified counter clockwise from the RIGHT with - respect to the RIGHT hemisphere: + Important: anchors must be specified counter clockwise from the RIGHT (per mathematical + convention that rightward is zero degrees) with respect to the LEFT (alphabetically + first) hemisphere: 2 /\ 3 <- center -> 1 @@ -1675,8 +1786,9 @@ def dartboard_on_flatmap(vertex_data, be used, e.g. ('STS', 'nearest') for the nearest point to the center_roi that falls on STS. - Important: anchors must be specified counter clockwise from the RIGHT with - respect to the RIGHT hemisphere: + Important: anchors must be specified counter clockwise from the RIGHT (per mathematical + convention that rightward is zero degrees) with respect to the LEFT (alphabetically + first) hemisphere: 2 /\ 3 <- center -> 1 @@ -1842,3 +1954,17 @@ def dartboard_on_flatmap(vertex_data, max_radii[0], max_radii[0]]) return fig + + +def average_outlines(outlines): + outlines = np.array(outlines) + for roi_i in range(outlines.shape[1]): + for hemi in range(2): + outlines[:, roi_i, hemi] = np.array( + pol2cart(*outlines[:, roi_i, hemi].T)).T + averaged_outlines = np.nanmean(outlines, axis=0) + for roi_i in range(averaged_outlines.shape[0]): + for hemi in range(2): + averaged_outlines[roi_i, hemi] = np.array( + cart2pol(*averaged_outlines[roi_i, hemi].T)).T + return averaged_outlines diff --git a/examples/dartboards/dartboard_low_level_computations_plots.py b/examples/dartboards/dartboard_low_level_computations_plots.py index 83682b75..1775c133 100644 --- a/examples/dartboards/dartboard_low_level_computations_plots.py +++ b/examples/dartboards/dartboard_low_level_computations_plots.py @@ -38,3 +38,12 @@ def mean_nonan(x, axis=None, threshold=0.8): # Get vertex-wise masks for each dartboard bin (`masks`) and masked data masks, data = cortex.dartboards.get_dartboard_data(vx, **dartboard_spec) + +# Get region of interest outlines. These will take some time to compute +# the first time, and will load from a cache thereafter (and thus be quick) +rois_to_plot = ['FEF','M1F', 'S1H', 'M1M','M1H'] +roi_outlines = {} +for roi in rois_to_plot: + roi_outlines[roi] = 0 + +# Plot masked data in dartboard plot From 506993aae139fd1d46f859b5bbc0076a616d366f Mon Sep 17 00:00:00 2001 From: marklescroart Date: Thu, 16 May 2024 10:20:11 -0700 Subject: [PATCH 15/15] ENH: updated anchor point finding --- cortex/dartboards.py | 137 +++++++++++++++++++++++++++++-------------- 1 file changed, 93 insertions(+), 44 deletions(-) diff --git a/cortex/dartboards.py b/cortex/dartboards.py index f19d22a3..2ff3d35b 100755 --- a/cortex/dartboards.py +++ b/cortex/dartboards.py @@ -160,6 +160,7 @@ def vertical_flip_path(path, overlay): def determine_path_hemi(overlay, path, percent_cutoff=.4): + """Returns 0 (False) for left, 1 (True) for right hemisphere""" if isinstance(path, matplotlib.path.Path): path = path.vertices middle_x = overlay.svgshape[0]/2 @@ -178,7 +179,9 @@ def center_of_mass(X): ------- center_of_mass array of shape (2,) that contains (x, y) center of mass - """ + """ + if not np.all(np.isclose(X[0], X[-1])): + X = np.vstack([X, X[0]]) # calculate center of mass of a closed polygon x = X[:, 0] y = X[:, 1] @@ -378,6 +381,7 @@ def _angles_from_vertex(overlay, start_idx, target_idxs=None, period=2*np.pi, th def get_roi_paths(overlay, roi, cleaned=True, filter_zeros=True, vertical_flip=True, overlay_file=None): + """returns left, right hemisphere paths""" if roi in list(overlay.rois.shapes.keys()): paths = overlay.rois[roi].splines elif roi in list(overlay.sulci.shapes.keys()): @@ -670,7 +674,7 @@ def show_dartboard(data, """ if max_radius is None: max_radius = 1 - data = np.array(data).astype(np.float) + data = np.array(data).astype(float) if isinstance(data2, np.ndarray): data = np.stack([data, data2], axis=0) else: @@ -748,6 +752,27 @@ def show_dartboard(data, return axis +def clockwise_test(pts): + """Test for clockwise ordering of points + + Parameters + ---------- + pts : array + 2D array (n_pts, xy) defining polygon points to be evaluated for clockwise rotation + + Returns + ------- + is_clockwise : bool + True for clockwise, False for counter-clockwise + """ + if not np.all(pts[0] == pts[-1]): + pts = np.vstack([pts, pts[0]]) + sum = 0 + for c0, c1 in zip(pts[:-1], pts[1:]): + sum += (c1[0] - c0[0]) * (c1[1] + c0[1]) + return sum > 0 + + def interpolate_outlines(phi, rho, resolution=50): new_phi = np.linspace(0,2*np.pi,resolution) new_rho = np.interp(new_phi, phi, rho, period=2*np.pi) @@ -873,40 +898,54 @@ def sort_roi_border(xy, angles, start_angle=0, start_index=None, max_deg_from_ze if start_index is None: start_index = np.argmin( np.abs(circ_dist(angles, start_angle, degrees=False))) - vertex_order = [start_index] - ct = 0 - for i_vert in range(len(xy)): - # Special case for first one: go clockwise - if i_vert == 0: - new_verts = np.argsort(dsts[start_index])[:5] - new_verts = np.array([x for x in new_verts if x not in vertex_order]) - angles_from_origin = circ_dist(angles[new_verts], angles[start_index], degrees=False) - is_clockwise = angles_from_origin > 0 - #if ~any(is_clockwise): - # plt.plot(*xy.T) - # plt.plot() - min_clockwise_angle = np.min(angles_from_origin[is_clockwise]) - j = new_verts[angles_from_origin == min_clockwise_angle][0] - vertex_order.append(j) - else: - success = False - n = 3 - while (not success) and (n<15): - new_verts = np.argsort(dsts[vertex_order[-1]])[:n] - new_verts = np.array([x for x in new_verts if x not in vertex_order]) - success = len(new_verts==1) - n += 1 - #if n > 4: - # print(n) - if len(new_verts) > 0: - vertex_order.append(new_verts[-1]) - else: - ct += 1 - if verbose: - print(f"Skipped a vertex ({ct} skipped)") - #1/0 - vertex_order.append(vertex_order[0]) - return np.array(vertex_order) + is_clockwise = clockwise_test(xy) + n = len(xy) + if is_clockwise: + print('Detected CLOCKWISE order') + vertex_order = np.hstack([np.arange(start_index, n), + np.arange(0, start_index), + start_index]) + else: + print('Detected COUNTER-CLOCKWISE order') + vertex_order = np.hstack([np.arange(start_index, 0, -1), + np.arange(n-1, start_index, -1), + start_index]) + + # vertex_order = [start_index] + # ct = 0 + # for i_vert in range(len(xy)): + # # Special case for first one: go clockwise + # if i_vert == 0: + # new_verts = np.argsort(dsts[start_index])[:5] + # new_verts = np.array([x for x in new_verts if x not in vertex_order]) + # angles_from_origin = circ_dist(angles[new_verts], angles[start_index], degrees=False) + # is_clockwise = angles_from_origin > 0 + # #if ~any(is_clockwise): + # # plt.plot(*xy.T) + # # plt.plot() + # min_clockwise_angle = np.min(angles_from_origin[is_clockwise]) + # j = new_verts[angles_from_origin == min_clockwise_angle][0] + # vertex_order.append(j) + # else: + # success = False + # n = 3 + # while (not success) and (n<15): + # new_verts = np.argsort(dsts[vertex_order[-1]])[:n] + # new_verts = np.array([x for x in new_verts if x not in vertex_order]) + # success = len(new_verts==1) + # n += 1 + # #if n > 4: + # # print(n) + # if len(new_verts) > 0: + # vertex_order.append(new_verts[-1]) + # else: + # ct += 1 + # if verbose: + # print(f"Skipped a vertex ({ct} skipped)") + # #1/0 + # vertex_order.append(vertex_order[0]) + # return np.array(vertex_order) + return vertex_order def get_interpolated_outlines(overlay, @@ -926,15 +965,18 @@ def get_interpolated_outlines(overlay, svgoverlay for a given subject outline_roi : str name of ROI to plot (must sexist in overlays.svg) - center_roi : str - name of center ROI for dartboard - anchor_angles : array - array of angles ?? geodesic_distances : bool, optional Flag for whether to compute geodesic distances for eccentricity bins, by default True resolution : int, optional number of points estimated for outline, by default 100 + every_n : int, optional + smooth path by fitting a cubic spline to `every_n` vertices on border + recache : bool, optional + whether to force re-computation (recache=True) or allow loading from cached + file (recache=False) + verbose : bool, optional + verbose output or not Returns ------- @@ -1163,8 +1205,8 @@ def show_dartboard_pair(dartboard_data, space, a little smoothing usually helps aesthetically. 1 is no smoothing , by default 5 even_sampling_over : str, optional - How to resample ROI paths, by angle or along path length. 'angle' is perhaps slightly - more principled for convex ROIs, but 'path_length' gives better results for non-convex + How to resample ROI paths, by angle or along path length. 'polar angle' is slightly + more principled for convex ROIs, but 'path length' gives better results for non-convex ROIs, by default 'path length' roi_linewidth : int, optional width of lines for drawn ROIs, by default 1 @@ -1566,10 +1608,17 @@ def _get_anchor_points(svg, center_roi, anchors, return_indices=True): anchor_name, anchor_type = anchor, 'centroid' if j == 0: center_name, center_type = anchor_name, anchor_type - if anchor_type == 'nearest': - centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name, return_indices=return_indices) + if 'nearest' in anchor_type: + if '_to_' in anchor_type: + _, _, roi_to = anchor_type.split('_') + centroids[f'{anchor_name}-{roi_to}'] = _get_closest_vertex_to_roi(svg, anchor_name, roi_to, return_indices=return_indices) + else: + centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name, return_indices=return_indices) elif anchor_type == 'centroid': centroids[anchor_name] = get_roi_centroids(svg, anchor_name, return_indices=return_indices) + elif anchor_type == 'superior': + pass + #centroids[anchor_name] = _get_closest_vertex_to_roi(svg, anchor_name, center_name, return_indices=return_indices) else: raise ValueError("unknown anchor type specified: %s\n(Must be 'nearest' or 'centroid')"%(anchor_type)) return centroids