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')