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