function [Xreturn,Yreturn,numcuts,Xl,Yl,Xr,Yr,cutindices, ncut] = cutoff(Xs,Ys,B,parallel) % This function checks for cutoffs by seeing where the banks intersect % themselves, and then performs the cut by removing the section of % centerline. The point where the % banks first intersect (looking downstream) is inserted into the % centerline. % % INPUTS: Xs,Ys - centerline coordinate vectors % B - centerline half-width (scalar) % % OUTPUTS: Xreturn,Yreturn - new centerline coordinate vectors after % cutoffs % numcuts - number of cutoffs performed in the time step % Xl,Yl and Xr,Yr - left and right bank coordinate vectors for % post-cutoff centerline %Xs=X;Ys=Y; % B=15; % Due to high curvatures (esp after a cutoff), the centerline can errantly intersect itself. If % this happens, cut off the loop. Eventually the curvatures will smooth % out. cutindices = []; ncut = []; [~,~,icl,jcl] = intersections(Xs,Ys,1); % returns intersections of centerline with itself while numel([icl; jcl]) > 0 cutoff = floor(icl(1)); % assumes that only one i,j pair is returned per intersection cuton = ceil(jcl(1)); Xs(cutoff:cuton) = []; Ys(cutoff:cuton) = []; cutindices = [cutindices, cutoff, cuton]; % for statistics module and node tracking ncut = [ncut, cuton-cutoff+1]; [~,~,icl,jcl] = intersections(Xs,Ys,1); end numsmooth = 10; % number of nodes before and after cutoff to use in smoothing spline SGnfilt = 5; % number of times to filter interpolated coordinates SGorder = 3; % order of the filter for interpolated coordinates -> lower order suppresses curvature more add_pts = 15; % number of points to add to interpolated section across cutoff % Reshape inputs Xs = Xs(:); Ys = Ys(:); % Return bank coordinates and intersection indices if parallel [Xl,Yl,Xr,Yr,ir,jr,X0r,Y0r,il,jl,X0l,Y0l] = gen_banks_par(Xs,Ys,B); else [Xl,Yl,Xr,Yr,ir,jr,X0r,Y0r,il,jl,X0l,Y0l] = gen_banks(Xs,Ys,B); end % Perform cutoffs numcuts = 0; while numel([ir;il]) > 0 preXs = length(Xs); numcuts = numcuts + 1; % Cutoffs are performed moving downstream. This can make a difference % if there are multiple cutoffs in a time step. if numel(ir) == 0 Xforcut = Xl; Yforcut = Yl; LR = 0; % LR is just a marker that lets the code know later if it's the right or left bank that was cut; LR=0 ==> Left bank if numel(il)== 1 LorRia = il(1); LorRib = NaN; LorRj = jl(1)-1; % the -1 is to make up for the +1 added to cuton later (exceeds length of stream if the intersection happens at the last node) else LorRia = il(1); LorRib = il(2); LorRj = jl(2); end elseif numel(il) == 0 Xforcut = Xr; Yforcut = Yr; LR = 1; if numel(ir) == 1 LorRia = ir(1); LorRib = NaN; LorRj = jr(1)-1; % the -1 is to make up for the +1 added to cuton later (exceeds length of stream if the intersection happens at the last node) else LorRia = ir(1); LorRib = ir(2); LorRj = jr(2); end elseif min(ir) < min(il) LorRia = ir(1); LorRib = ir(2); LorRj = jr(2); Xforcut = Xr; Yforcut = Yr; LR = 1; else LorRia = il(1); LorRib = il(2); LorRj = jl(2); Xforcut = Xl; Yforcut = Yl; LR = 0; end if LorRib - LorRia == 1.5 % this lets the code know that it's dealing with an end-of-the-river intersection, pretty rare case but can happen. no smoothing necessary cutoff = floor(LorRia); cutindices = [cutindices, cutoff, length(Xs)]; % for statistics module and node tracking, the cuton is added later ncut = [ncut, length(Xs)-cutoff]; Xs = Xs(1:cutoff); Ys = Ys(1:cutoff); else cutoff = floor(LorRia); % finds index at which banks intersect first cuton = ceil(LorRj)+1; %finds index at which banks intersect last cutindices = [cutindices, cutoff, cuton]; % for statistics module and node tracking, the cuton is added later % Extract the centerline around the cutoff for interpolation and % smoothing 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 Xsforsmooth = double([Xs(cutoff-numsmooth_beg:cutoff); Xs(cuton:cuton+numsmooth_end)]); % have to convert to double to suppress a warning about mixed data types from ode45 in interparc. Ysforsmooth = double([Ys(cutoff-numsmooth_beg:cutoff); 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(Xsforsmooth)/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)]; % Get new banklines and cutoff indices after performing cutoff finalXs = length(Xs); ncut = [ncut, preXs - finalXs]; % for statistics module and node tracking end [Xl,Yl,Xr,Yr,ir,jr,X0r,Y0r,il,jl,X0l,Y0l] = gen_banks(Xs,Ys, B); end Xreturn = single(Xs); % use single precision to save memory. in spot-checks it made absolutely no detectable difference but cuts file sizes (and therefore memory use) in half Yreturn = single(Ys); end