From 506993aae139fd1d46f859b5bbc0076a616d366f Mon Sep 17 00:00:00 2001 From: marklescroart Date: Thu, 16 May 2024 10:20:11 -0700 Subject: [PATCH] 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