function [Xs, Ys, cutidcs] = smooth_after_cutoffs(Xs, Ys, cutidcs, Npre)
% Npre is the number of points in the centerline before any cutoffs were
% removed.
numsmooth = 15; % number of nodes before and after cutoff to use in smoothing
SGnfilt = 1; % number of times to filter interpolated coordinates
SGorder = 5; % order of the filter for interpolated coordinates -> lower order suppresses curvature more
add_pts = 2; % number of points to add to interpolated section across cutoff
% keyboard
ncuti = 0;
smoothidcs = cutidcs(:,1)-1;
for i = 1:numel(smoothidcs) % this loop finds the index of the node immediately upstream of each cutoff; indices must be adjusted for upstream cutoffs
if i == 1
else
ncuti = ncuti + (cutidcs(i-1,2)-cutidcs(i-1,1)+1);
smoothidcs(i) = smoothidcs(i) - ncuti;
end
end
M = numel(smoothidcs);
for cutno = M:-1:1 % work upstream to downstream
if cutidcs(cutno,1) == 1 % if the cutoff occurs at the beginning of the centerline
elseif cutidcs(cutno,2) == Npre % if the cutoff occurs at the end of the centerline
else
cutoff = smoothidcs(cutno);
cuton = cutoff + 1;
% interpolate two nodes linearly across the cutoff
x1 = Xs(cutoff); x2 = Xs(cuton);
y1 = Ys(cutoff); y2 = Ys(cuton);
dx = x2 - x1;
dy = y2 - y1;
xi1 = x1 + (1/3)*dx;
xi2 = x1 + (2/3)*dx;
yi1 = y1 + (1/3)*dy;
yi2 = y1 + (2/3)*dy;
numsmooth_beg = min(numsmooth, cutoff-1); % in case the cutoff occurs near the beginning of the river
numsmooth_end = min(numsmooth, length(Xs)-cuton-1); % in case the cutoff occurs near the end of the river
Xsmooth = double([Xs(cutoff-numsmooth_beg:cutoff); xi1; xi2; Xs(cuton:cuton+numsmooth_end)]); % have to convert to double to suppress a warning about mixed data types from ode45 in interparc.
Ysmooth = double([Ys(cutoff-numsmooth_beg:cutoff); yi1; yi2; Ys(cuton:cuton+numsmooth_end)]);
% Interpolate across the cutoff
% [XY] = interparc(numsmooth*2+add_pts,Xsforsmooth,Ysforsmooth,'spline'); % s_pts is a function of the streamlength, see further in code. The initial ratio of pts/length is maintained.
% Xsmooth = single(XY(:,1)); % convert back to a single
% Ysmooth = single(XY(:,2));
% Use Savitzky-Golay filter to smooth interpolated coordinates around cutoff
window = max(2.*round((SGorder+1+1)/2)-1,2.*round((length(Xsmooth)/4+1)/2)-1); % window must be larger than order of filter, hence the max(). must also be odd, hence the rounding.
tmp = [Xsmooth, Ysmooth];
for iFilt=1:SGnfilt;
tmp = sgolayfilt(tmp,SGorder,window);
end
Xsmooth = tmp(:,1);
Ysmooth = tmp(:,2);
Xs = [Xs(1:cutoff-numsmooth-1); Xsmooth; Xs(cuton+numsmooth+1:end)];
Ys = [Ys(1:cutoff-numsmooth-1); Ysmooth; Ys(cuton+numsmooth+1:end)];
end
cutidcs(cutno,2) = cutidcs(cutno,2) + add_pts; % account for added nodes
if cutidcs(cutno,2) > numel(Xs) % account for cutoffs near the end but not exactly at end of centerline
cutidcs(cutno,2) = numel(Xs);
end
end % for
end % function
% clf
% plot(Xs,Ys,'.');axis equal; hold on
% plot(Xs(cutoff),Ys(cutoff),'ro')