function [Xl,Yl,Xr,Yr,ir,jr,X0r,Y0r,il,jl,X0l,Y0l] = gen_banks(Xsb,Ysb,B)
% This fucntion generates banks for a given centerline and constant width.
% INPUTS: Xsb,Ysb - centerline coordinate vectors
% B - channel half-width
% OUTPUTS: Xl,Yl and Xr,Yr - left and right bank coordinate vectors
% il,jl and ir,jr - left and right bank indices of where cutoffs
% occur
% X0l,Y0l and X0r,Y0r - left and right bank intersection points
% corresponding to cutoffs found in i,j
%
% There are two methods to generate banks.
% METHOD 1: Uses sin(pi-angle)*B to
% calculate the new banks. This method underestimates the bank width
% everywhere, moreso at higher curvatures. However it's much faster and
% does a decent job so it's currently implemented.
%
% METHOD 2: This code generates banks a distance B away from an input centerline
% given by Xsb,Ysb. It does so by creating lines parallel to the vectors
% defining the centerline that are B distance away in a perpendicular
% sense. The intersections of these parallel lines define the vetices of
% the banks. This method is currently pasted at the end of this file.
%
% Both methods fail when curvature is high (e.g.immediately after a cutoff)
% causing the bank to erroneously intersect itself near the apex of the
% bend, sometimes multiple times. Extra coding was required to distinguish
% between an erroneous cutoff and a legitimate one.
%
% The code returns the indices of actual cutoffs and the left and right
% bank coordinates after erroneous bank intersections have been removed.
% Xsb=Xs;Ysb=Ys; % for debugging/testing
% Reshape inputs
Xsb = Xsb(:);
Ysb = Ysb(:);
thresh = 7;
thresh2 = 10;
checkthresh = 5; % threshold for distinguishing between pairs of cutoff vertices. if a whole section of bank is removed (obviously) incorrectly, increase this threshold
min_ox_length = 20;
% n_gen_line = 10^5; % just needs to be large enough so that the line segments overlap
[angles] = cl_angles(Xsb,Ysb); %angles computed under assumption of constant valley direction found via linear regression!!!
Xl = Xsb-sin(pi-angles)*B;
Yl = Ysb-cos(pi-angles)*B;
Xr = Xsb+sin(pi-angles)*B;
Yr = Ysb+cos(pi-angles)*B;
% Loop to generate left and right banks
for lr = 1:2
if lr == 1; % Determines left or right bank
Xbank = Xr;
Ybank = Yr;
B = B; % right bank
else
Xbank = Xl;
Ybank = Yl;
B = -B; % left bank
end
% Now that the bank has been generated, need to take care of erroneous
% cases.
[X0fix,Y0fix,ifix,jfix] = intersections(Xbank,Ybank,1); % returns intersections of left bank with itself
% This if statment checks for different combinations of legitmate
% cutoffs mixed in with erroneous bank intersections
if isempty(ifix) == 1
elseif length(ifix) == 1 && (length(Xbank) - floor(jfix)) < 3 % case where the end of the river is intersecting (only one intersection point in this case, hence the isodd requirement)
ifix = [ifix ifix+1.5]; % adding 1.5 as a marker to the cutoff.m routine that this is an endpoint
jfix = [jfix-1.5 jfix]; % adding 1.5 as a marker to the cutoff.m routine that this is an endpoint
X0fix = [X0fix X0fix];
Y0fix = [Y0fix Y0fix];
elseif length(ifix) == 1 % case where there is only one bank intersection (must be a bank error in this case UNLESS the cutoff points intersect EXACTLY at the same node, extremely low likelihood)
cutoff = floor(ifix);
cuton = ceil(jfix);
Xbef = Xbank(1:cutoff);
Xaft = Xbank(cuton:end);
Xbank = [Xbef; X0fix; Xaft];
Ybef = Ybank(1:cutoff);
Yaft = Ybank(cuton:end);
Ybank = [Ybef; Y0fix; Yaft];
ifix = [];
jfix = [];
else % case where valid cutoffs are mixed in with SINGLE CROSS bank error intersection
keep = find(jfix-ifix > min_ox_length); % if the cutoff section is fewer than min_ox_length number of vertices, it's erroneous
ijXY = [ifix jfix X0fix Y0fix]; % puts all the intersections in matrix form for less coding
remmat = ijXY; % makes a copy of the intersection matrix
remmat(keep,:) = []; % leaves only the intersections that should be removed
keepmat = ijXY(keep,:); % a matrix of the intersections to keep; each row corresponds to an intersection. should be even number of rows!
ifix = keepmat(:,1); % leaves all the coupled intersection indices that are actual cutoffs
jfix = keepmat(:,2);
X0fix = keepmat(:,3);
Y0fix = keepmat(:,4);
icheck = floor(remmat(:,1));
jcheck = ceil(remmat(:,2));
% These loops are required to handle cases where the returned
% intersections are contained within one another; that is, an
% intersection may be at i=20 and j=30, but there's another one at
% i=22 and j=27 and another one at i = 19 and j = 29, for example.
for xx = 1:length(icheck)
for intno = 1:size(remmat,1)
if icheck(xx) > floor(remmat(intno,1)) && icheck(xx) < ceil(remmat(intno,2))
end
if jcheck(xx) > floor(remmat(intno,1)) && jcheck(xx) < ceil(remmat(intno,2))
end
end
end
for g = 1:length(icheck)-1
if icheck(g+1) == 0 && jcheck(g+1) == 0
elseif icheck(g+1) == 0
if jcheck(g+1) > jcheck(g) && abs(jcheck(g+1)-jcheck(g)) < checkthresh
jcheck(g) = jcheck(g+1);
end
jcheck(g+1) = 0;
elseif jcheck(g+1) == 0
if icheck(g+1) < icheck(g) && abs(icheck(g+1)-icheck(g)) < checkthresh
icheck(g) = icheck(g+1);
end
icheck(g+1) = 0;
end
end
ijcheck = [icheck jcheck];
ijcheck(icheck==0 | jcheck==0,:)=[];
icheck = ijcheck(:,1);
jcheck = ijcheck(:,2);
if length(icheck)~=length(jcheck)
disp('error in check vectors, gen banks')
end
removeindices = [];
for r = 1:length(icheck)
removeindices = [removeindices icheck(r):jcheck(r)];
end
removeindices = unique(removeindices);
% Now that the indices to remove have been found, remove them.
Xbank(removeindices) = [];
Ybank(removeindices) = [];
end
% Save the coordinates and intersections
if lr == 1
Xr = Xbank;
Yr = Ybank;
ir = ifix;
jr = jfix;
X0r = X0fix;
Y0r = Y0fix;
else
Xl = Xbank;
Yl = Ybank;
il = ifix;
jl = jfix;
X0l = X0fix;
Y0l = Y0fix;
end
clear Xbank Ybank res ifix jfix X0fix Y0fix cutout n_vertices_btwn
end % l/r end
end % function end
%% Method 2
%--- This code could probably be vectorized so that intersections isn't
% repeatedly called ---%
% % Generate parallel line vectors
% for k = 1:length(Xos)
% res(k,:) = parallelLine([Xos(k), Yos(k), dx(k), dy(k)], B);
% end
% % Save the intersections of the parallel lines as banks
% for b = 1:length(res)-1
% X1 = [res(b,1)-n_gen_line*res(b,3), res(b,1) + n_gen_line * res(b,3)];
% Y1 = [res(b,2)-n_gen_line*res(b,4), res(b,2) + n_gen_line * res(b,4)];
% X2 = [res(b+1,1)-n_gen_line*res(b+1,3), res(b+1,1) + n_gen_line * res(b+1,3)];
% Y2 = [res(b+1,2)-n_gen_line*res(b+1,4), res(b+1,2) + n_gen_line * res(b+1,4)];
% [X_int,Y_int,~,~] = intersections(X1,Y1,X2,Y2); % returns intersections of the overlapping parallel line segments
%
% Xbank(b,1) = X_int;
% Ybank(b,1) = Y_int;
% end
%% Attempt to remove erroneous cutoffs via curvature thresholding
% It worked but was dependent on thresholds and the resulting banks were
% non-realistic shortly following cutoffs. Also when a cutoff hasn't
% occurred recently, the curvature threshold mistakenly identifies low
% curvatures as outliers (because it's based on the standard deviation of
% all curvatures).
% min_c_thresh = 3; % number of standard deviations to identify a curvature associated with an improperly calculated bank
% [angles] = cl_angles(Xs,Ys); %angles computed under assumption of constant valley direction found via linear regression!!!
%
% Xl = Xs-sin(pi-angles)*B;
% Yl = Ys-cos(pi-angles)*B;
% Xr = Xs+sin(pi-angles)*B;
% Yr = Ys+cos(pi-angles)*B;
% clf; plot(Xs,Ys); hold on; plot(Xr,Yr,'r'); plot(Xl,Yl,'k');axis equal
% [Csr, ~, ~, ~] = curvars(Xr,Yr);
% [Csl, ~, ~, ~] = curvars(Xl,Yl);
%
% if sum(abs(Csl*B) > min_c_thresh*abs(std(Csl*B))) > 0 % if the threshold is exceeded for the left bank
% [Xl, Yl] = bankfix(Xl,Yl,Csl,B,min_c_thresh);
% end
% if sum(abs(Csr*B) > min_c_thresh*abs(std(Csr*B))) > 0 % if the threshold is exceeded for the left bank
% [Xr, Yr] = bankfix(Xr,Yr,Csr,B,min_c_thresh);
% end
%
%
% [XY] = interparc(200,testX,testY,'linear'); % s_pts is a function of the streamlength, see further in code. The initial ratio of pts/length is maintained.
%
%
% % Savitzky-Golay Filter
% nFilt = 1;
% order = 3;
% window = 5; % automatically select window size based on length of signal, must be odd integer
% tmp = [testX, testY];
% for iFilt=1:nFilt;
% tmp = sgolayfilt(tmp,order,window);
% end
% X = tmp(:,1);
% Y = tmp(:,2);
% Set threshold levels
% t1 = std(Cs*B)*min_c_thresh;
% t2 = std(Cs*B)*4;
% t3 = std(Cs*B)*6;
% Loop to find the number of points to remove on either side of a "greater
% than threshold" curvature.
% Left bank
% for ll = 1:length(Cs_reml)
% if abs(Cs_reml(ll)*B) > t3 % this line says "if the nondimensional curvature is 6 times higher than the nondimensional std of curvatures, remove 4 zeros on either side. the next lines are similar.
% num_zeros_l(ll) = 4;
% elseif abs(Cs_reml(ll)*B) > t2
% num_zeros_l(ll) = 3;
% elseif abs(Cs_reml(ll)*B) > t1
% num_zeros_l(ll) = 2;
% else
% disp('impossibru!')
% end
% end
%
% % Right bank
% for rr = 1:length(Cs_remr)
% if abs(Cs_remr(rr)*B) > t3 % this line says "if the nondimensional curvature is 6 times higher than the nondimensional std of curvatures, remove 4 zeros on either side. the next lines are similar.
% num_zeros_r(rr) = 4;
% elseif abs(Cs_remr(rr)*B) > t2
% num_zeros_r(rr) = 3;
% elseif abs(Cs_remr(rr)*B) > t1
% num_zeros_r(rr) = 2;
% else
% disp('impossibru!')
% end
% end
% % Bank generation method fails near high-curvature regions, leading to
% % false cutoff identification.
% [Cs, ~, ~, ~] = curvars(Xr,Yr);
% n_keepl = Cs>min_c_thresh*-std(Cs); % logical of points not exceeding threshold; for left bank, curvatures are negative
% n_keepr = Cs 0 % if there are erroneous bank intersections
% remcol = floor([remmat(:,1) [1:size(remmat,1)]'; remmat(:,2) [1:size(remmat,1)]']); % arrange the indices into a column with corresponding row marker. unnecessary but makes coding more intuitive.
% c=1; rowrow = [];
% for ii = 1:size(remcol,1)
% val = remcol(ii,1);
% rowmark = remcol(ii,2);
% for row = 1:size(remmat,1)
% below = floor(remmat(row,1));
% above = floor(remmat(row,2));
% if val > below && val < above
% rowrow(c,1) = rowmark; % the rowrow matrix is defined as [which row is there an in-between value, which row that value is an in-between of]
% rowrow(c,2) = row;
% c = c+1;
% end
% end
% end
% rowrow = unique(rowrow,'rows'); % still need to deal with [4,3; 3,4]
% removetheseinbtw = [];
% removethesesolo = [];
% if numel(rowrow) > 0
% for jj = 1:size(rowrow,1);
% ireminbtw(jj) = floor(min(remmat(rowrow(jj,1),1),remmat(rowrow(jj,2),1)));
% jreminbtw(jj) = ceil(max(remmat(rowrow(jj,1),2),remmat(rowrow(jj,2),2)));
% end
% for kk = 1:length(ireminbtw)
% removetheseinbtw = [removetheseinbtw, ireminbtw(kk):jreminbtw(kk)];
% end
% end
% remmat(rowrow,:) = []; % remove the intersections that are already accounted for
% iremsolo = floor(remmat(:,1));
% jremsolo = ceil(remmat(:,2));
% for mm = 1:length(iremsolo)
% removethesesolo = [removethesesolo, iremsolo(mm):jremsolo(mm)];
% end
% removethese = unique([removetheseinbtw; removethesesolo]);
%
% Xbank(removethese) = [];
% Ybank(removethese) = [];
% end
%
% ifix = keepmat(:,1); % leaves all the coupled intersection indices that are actual cutoffs
% jfix = keepmat(:,2);
% X0fix = keepmat(:,3);
% Y0fix = keepmat(:,4);
% % Now loop through the single, non-coupled intersections and fix
% % the banks there.
% for w = 1:length(irem)
% cutoff = floor(toremovei(w));
% cuton = floor(toremovej(w))+1;
% Xbef = Xbank(1:cutoff);
% Xaft = Xbank(cuton:end);
% Xbank = [Xbef; toremoveX0(w); Xaft];
% Ybef = Ybank(1:cutoff);
% Yaft = Ybank(cuton:end);
% Ybank = [Ybef; toremoveY0(w); Yaft];
% end
% end
%
% ifix = keepmat(:,1); % leaves all the coupled intersection indices
% jfix = keepmat(:,2);
% X0fix = keepmat(:,3);
% Y0fix = keepmat(:,4);
%
% % So at this point, there should be an even number of elements in the
% % i and j (ifix and jfix) vectors.
% if numel(ifix) == 0 % if there are no intersections, skip this block
% elseif mod(numel(ifix),2) == 1 % if it's odd, we have a problem
% disp('problem in gen_banks')
% else
% % Now distinguish the intersection pairs that are actual cutoffs
% % versus erroneous bank calcs by thresholding the number of
% % vertices between the first and second intersection
% ifixfirst = floor(ifix(1:2:end));
% jfixsec = ceil(jfix(2:2:end));
% X0fixfirst = X0fix(1:2:end);
% Y0fixfirst = Y0fix(1:2:end);
% n_vertices_btwn = jfixsec - ifixfirst;
% cutoff = ifixfirst(n_vertices_btwn < thresh2);
% cuton = jfixsec(n_vertices_btwn < thresh2);
% % Replaces the erroneous cutoff points with the intersection
% % point
% for cc = 1:length(cutoff)
% Xbank(cutoff(cc):cuton(cc)) = X0fixfirst(cc);
% Ybank(cutoff(cc):cuton(cc)) = Y0fixfirst(cc);
% end
% % Deletes the repeated inserted points
% [value, num] = mode(Xbank);
% if num > 1
% repeats = find(Xbank==value);
% repeats = repeats(2:end); % save one of the points, ie don't delete them all
% Xbank(repeats) = [];
% Ybank(repeats) = [];
% end
%
% for vv = 1:length(n_vertices_btwn)
% if n_vertices_btwn(vv) < thresh2
% cutout(2*vv-1) = 0;
% cutout(2*vv) = 0;
% else
% cutout(2*vv-1) = 1;
% cutout(2*vv) = 1;
% end
% end
% ifix = ifix(cutout==1);
% jfix = jfix(cutout==1);
% X0fix = X0fix(cutout==1);
% Y0fix = Y0fix(cutout==1);