function [Xs, Ys, cutidcs] = remove_cutoffs(Xs, Ys, cut_thresh_in_Bs, B) % This code works in a downstream direction along a set of centerline % coordinates to detect and remove cutoffs by finding non-adjacent nodes % that are some threshold distance (search_radius) near each other. % Jon Schwenk, 8/20/2014. jonschwenk@gmail.com % INPUTS: Xs, Ys - coordinates of pre-cut channel centerline % cut_thresh_in_Bs - minimum distance at which cutoff occurs, % given in channel half-widths (for purely % neck cutoffs, this should be channel width) % B - channel half-width % OUTPUTS: Xs, Ys - channel centerline coordinates with cutoffs removed % cutidcs - Mx2 array, where M = number of cutoffs identified. % cutidcs(M,1) = upstream cutoff node % cutidcs(M,2) = downstream cutoff node. len_thresh = 10*B; % two nodes must be separated by this along-stream distance to be considered a cutoff pair (must be smaller than minimum possible atom length). This threshold helps exclude adjacent centerline nodes. search_radius = cut_thresh_in_Bs * B; % ensure column vectors Xs = Xs(:); Ys = Ys(:); % compute spacing between nodes dS = sqrt(diff(Xs).^2+diff(Ys).^2); S = [0; cumsum(dS)]; p = [Xs Ys]; [idx, dists] = rangesearch(p, p, search_radius, 'Distance', 'euclidean', 'NSMethod', 'kdtree'); % returns for each node all neighbors within the search radius and their corresponding distances cutidcs = []; N = numel(Xs); search_idx = 1; % finds the cutoff indices % a while loop is used to speed up the cutoff search by skipping the nodes % of an identified cutoff while search_idx ~= N; % diffs = idx{search_idx} - search_idx; % how many indices away are the ones within the search radius? % thresh_remove = diffs > idx_thresh; diffsl = S(idx{search_idx}) - S(search_idx); thresh_remove = diffsl > len_thresh; if sum(thresh_remove) > 0 % if any nodes are possible cutoff ones poss_idcs = idx{search_idx}(thresh_remove); % only take those that are at least len_thresh away poss_dists = dists{search_idx}(thresh_remove); % only take those that are at least idx_thresh away cutidx = poss_idcs(poss_dists==min(poss_dists)); % take the closest node as the cutoff node % save the found cutoff indices cutidcs = [cutidcs; search_idx cutidx]; search_idx = cutidx; % skips searching the nodes along the cutoff if search_idx == N break end end search_idx = search_idx + 1; end % perform the cutoff (ie remove the nodes from the centerline) if isempty(cutidcs) == 0 remove = []; for cutno = 1:size(cutidcs,1); remove = [remove cutidcs(cutno,1):cutidcs(cutno,2)]; end Xs(remove) = []; Ys(remove) = []; end end % function