%parse_Transect_1_Channel_Width_Measurements_Second_Pass_Points_for_LPSC_2014

%************************SET CONTROL PARAMETERS ***************************
WhichPass = 3 %Refers to the first, second, third passes of Transect 1 picks (widths and wavelengths)
distancecutoff = 500; %do not evaluate line segments more distant than this (for channel widths)
dottingstrategy = [0 1]; %if =1, dot uniformly along channel banks; if =0, just find width from vertices.
dbd = 2.5; %distance between dots (m)

%************************SET GRAPHICS AND HOUSEKEEPING PARAMETERS ***************************
picturesque = 0;
output_type = 'terse'; %'none', 'terse', or 'complete'
%************************IMPORT DATA FROM CSV FILES ***************************
load_channel_widths
%************************END OF IMPORT DATA FROM CSV FILES ***************************

assert(min(px)==0); assert(min(py)==0); assert(minx~=0); assert(miny~=0);
%minx = min(px); miny = min(py);
%px = px - minx; 
%py = py - miny;

%check image:
scatter(px,py,10,z);

%check image:
if exist('rootsemi','var') == 1
    figure
    scatter(px,py,10,rootsemi);
end
    
%********* Pepper Points Uniformly Along Channel Centerlines **************
if sum(dottingstrategy==1)>=1
        ctl(1:max(bid)) = 0; %bank total length
        %Interpolate uniformly in s, dir coordinates
        for l = 1:(size(px,1)-1)
            if bid(l) == bid(l+1)%Is it from the same bank?
                lseg(l)     =   sqrt((px(l+1) - px(l)).^2 + ...
                                     (py(l+1) - py(l)).^2); %length of segment
                ctl(bid(l)) = ctl(bid(l)) + lseg(l); %bank total length
                lac(l+1)    = ctl(bid(l)); %length along bank for this vertex
            else %last point on channel trace
                %lac(l) = ctl(cid(l));
            end
        end
        
        assert(min(ctl(ctl>0))>dbd,'Error. Decrease distance between dots (dbd)')

        for m = 1:max(bid) 
            if sum(bid==m) >= 1
            %position vector of interpolated points
            p(m).x = interp1(lac(bid==m),px(bid==m),[0:dbd:ctl(m)]);
            p(m).y = interp1(lac(bid==m),py(bid==m),[0:dbd:ctl(m)]);
            p(m).z = interp1(lac(bid==m), z(bid==m),[0:dbd:ctl(m)]);
            p(m).bid(1:max(size(p(m).x)))=m;
            end
        end
end

%*********************************FIND CHANNEL WIDTHS ***************************
dcls_temp = NaN;
%*********************************FIND CHANNEL WIDTHS ***************************

for strategy_index = 1:length(dottingstrategy)
    dotuniformly = dottingstrategy(strategy_index)
    
    if dotuniformly == 0
        nstartpoints = size(px,1);
    elseif dotuniformly == 1
        nstartpoints = max(size([p(:).bid]));
        udx = [p(:).x]; udy = [p(:).y]; udz = [p(:).z]; udbid = [p(:).bid]; %uniformly-dotted x, y, z, and bank-ID
    end
    
%if dotuniformly == 0
for l = 1:nstartpoints %for each vertex
    if rem(l,100) == 0;disp(strcat('dotting strategy:',num2str(dotuniformly), ...
                                   ',l:',num2str(l)));end
    %if bid(min(size(px,1),l+1)) == bid(l) && bid(max(1,l-1)) == bid(l) ...
    %   && l > 1 && l < size(px,1)   %omit endpoints of banks
    
    if l > 1 && l < nstartpoints   %include endpoints of banks
    %find distance and coordinates of closest line segment
        if dotuniformly == 0
            x3 = px(l);y3 = py(l);
        elseif dotuniformly ==1
            x3 = udx(l);y3 = udy(l);
        end
        
        pd = sqrt((x3 - px).^2 + (y3 - py).^2);
        %pd(1) = distancecutoff*9e9; 
        pd(end) = distancecutoff*9e9; %prevent evaluation of first and last m-points
        xc = NaN*zeros(size(px));
        yc = xc;
        
        okmlist = find(pd<distancecutoff);
        
    for mindex = 1:length(okmlist) %2:(size(px,1)-1) 
           m = okmlist(mindex);
           %if pd(m) < distancecutoff
           %per http://stackoverflow.com/questions/6176227/for-a-point-in-an-irregular-polygon-what-is-the-most-efficient-way-to-select-th
           %but with sqrt added to denominator of equation for u(m) 
            x1 = px(m)  ; y1 = py(m);
            x2 = px(m+1); y2 = py(m+1);
        
        %segment is finite
            dQ1(m) = sqrt((x1 - x3)^2 + (y1 - y3)^2); %<-- Bugs result if these lines are moved inside the "filter out line segment between different bank ID's" loop
            
            dQ2(m) = sqrt((x2 - x3)^2 + (y2 - y3)^2);
            
        %DO consider the first line segment (for which bid(m-1)~=bid(m))
        if dotuniformly == 0
            acceptable_endpoint = (bid(m+1) == bid(m) && l~=m && bid(m)~=bid(l));
        elseif dotuniformly ==1
            acceptable_endpoint = (bid(m+1) == bid(m) && bid(m)~=udbid(l));
        end
            
        if acceptable_endpoint == 1 
      
%Following http://paulbourke.net/geometry/pointlineplane/
%using lambda instead of u and PC instead of P
lambda(m) = ( (x3-x1).*(x2-x1) + (y3-y1).*(y2-y1) ) ./ ...
             ((x2-x1).^2 + (y2-y1).^2);


if lambda(m)<0
    xc(m) = x1; yc(m) = y1;
    lambda(m) = 0;
elseif lambda(m)>1
    xc(m) = x2; yc(m) = y2;
    lambda(m) = 1;
else
    xc(m) = x1 + lambda(m)*(x2-x1);
    yc(m) = y1 + lambda(m)*(y2-y1);
end
          
        else
            %skip. Zero values of xc will be used to set dcls_temp(m) to NaN
        end
            %skip. Zero values of xc will be used to set dcls_temp(m) to NaN
    end
  
 dcls_temp = sqrt((x3 - xc).^2 + (y3 - yc).^2);        
 dcls_temp(isnan(xc)) = NaN; % Zero values of xc are used to set dcls_temp(m) to NaN
    
    [C,I] = nanmin(dcls_temp);
    if dotuniformly == 0
        lambdaforl(l) = lambda(I);
        clx(l) = xc(I);
        cly(l) = yc(I);
        dcls(l) = C;
        cls(l) = I; %line segment number, not point number
        clsbank(l) = bid(I);
    elseif dotuniformly == 1
        udlambdaforl(l) = lambda(I);
        udclx(l) = xc(I);
        udcly(l) = yc(I);
        uddcls(l) = C;
        udcls(l) = I; %line segment number, not point number
        udclsbank(l) = bid(I);
    end
    
     %We want to suppress the case where the closest *point* (not line segment) is a
            %vertex at the END of a bank (banks are indexed by 'bid')
                dQdebug(l) = 0;
            if dQ1(I) == dcls_temp(I) && bid(I-1)~=bid(I)
                dQdebug(l) =1;
            elseif dQ2(I) == dcls_temp(I) && bid(min(length(bid),I+2))~=bid(I)
                dQdebug(l) =2;
            end
            
            if dQdebug(l) > 0
                if dotuniformly == 0
                    dcls(l) = NaN; cls(l) = NaN; clsbank(l) = NaN;
                elseif dotuniformly == 1
                    udcls(l) = NaN; uddcls(l) = NaN; udclsbank(l) = NaN;
                end
            end
    
    else
        if dotuniformly == 0
            cls(l) = NaN; dcls(l) = NaN; clsbank(l) = NaN;
        else
            udcls(l) = NaN; uddcls(l) = NaN; udclsbank(l) = NaN;
        end
    end
end
end

%***********END OF FIND CHANNEL WIDTHS ***************************

if picturesque == 1
    figure %check plot
    scatter(px,py,'k*')
    hold on
    title('plot for debugging')
    EdwinColor = {'k','b','r','g','y',[.5 .6 .7],[.8 .2 .6]};
            % for l = 2:(size(px,1)-1)
            %     if isnan(cls(l)) == 0
            %         line([px(l) clx(l)],[py(l) cly(l)],'Color', EdwinColor{1+rem(l,6)})
            %         %line([px(l) px(cls(l))],[py(l) py(cls(l))],'Color', EdwinColor{1+rem(l,6)})
            %         scatter(px(l),py(l),20,EdwinColor{1+rem(l,6)},'LineWidth',2)
            %     end
            % end
            %     scatter(px(1:length(dQdebug)),py(1:length(dQdebug)),1500,dQdebug);

%     for l = 2:max(size(udclx))
%         if isnan(udcls(l)) == 0
%             line([udx(l) udclx(l)],[udy(l) udcly(l)],'Color', EdwinColor{1+rem(l,6)})
%             %line([px(l) px(cls(l))],[py(l) py(cls(l))],'Color', EdwinColor{1+rem(l,6)})
%             scatter(udx(l),udy(l),20,EdwinColor{1+rem(l,6)},'LineWidth',2)
%         end
%     end

            %    scatter(px(1:length(dQdebug)),py(1:length(dQdebug)),50,dQdebug);
end 
    
%**********AGGREGATE POINT-BY-POINT MEASUREMENTS TO FIND CHANNEL FIRST MOMENTS ***************************

%now find median width, dz, and rootsemi for each channel
    %chann = floor((bid+2)/2); %channel number
    
    %Note 21FEB2014: This section should be cleaned up so that all
    %measurements go into the 's' (strategy) structure. Check lower down in
    %the script where some of the 'cx', 'cy', e.t.c. are used for plotting.
for nn = 1:(max(bid)/2) %each channel has two banks; they are assumed to be in order of this pairing!
    
    css(nn)  = sum(bid==((nn*2) ) | bid==((nn*2)-1));%channel sample size (for calculating pooled variance, later)
    if css(nn)>0
        cx(nn)   = nanmedian(px(bid==(nn*2) | bid==((nn*2)-1)));
        cy(nn)   = nanmedian(py(bid==(nn*2) | bid==((nn*2)-1)));
        cqual(nn) = nanmedian(pqual((bid==(nn*2) | bid==(nn*2)-1))); %note that FID == a(:,2)
        cz(nn)    = nanmedian(z((bid==(nn*2) | bid==(nn*2)-1) & isnan(dcls')==0));
        czdev(nn) = nanstd(z((bid==(nn*2) | bid==(nn*2)-1) & isnan(dcls')==0));
        
        cwhichdtm(nn) = nanmean(whichdtm(bid==(nn*2) | bid==(nn*2)-1));
        cpstyle(nn) = nanmean(pstyle((bid==(nn*2) | bid==(nn*2)-1)));

        if sum(dottingstrategy==0)>=1
            banktobankdistances = dcls(bid==(nn*2) | bid==((nn*2)-1));
            s(10).cw(nn)   = nanmedian(banktobankdistances);
            s(10).cwrange(nn) = max(banktobankdistances)./min(banktobankdistances);
            s(10).cwmin(nn) = nanmin(banktobankdistances);
            s(10).cwmax(nn) = nanmax(banktobankdistances);
            s(10).cwdev(nn) = nanstd(dcls(bid==(nn*2) | bid==(nn*2)-1));
        end
        
        if sum(dottingstrategy==1)>=1
            banktobankdistances = uddcls(udbid==(nn*2) | udbid==((nn*2)-1));
            s(11).cw(nn)   = nanmedian(banktobankdistances);
            s(11).cwrange(nn) = max(banktobankdistances)./min(banktobankdistances);
            s(11).cwmin(nn) = nanmin(banktobankdistances);
            s(11).cwmax(nn) = nanmax(banktobankdistances);
            s(11).cwdev(nn) = nanstd(uddcls(udbid==(nn*2) | udbid==(nn*2)-1));
        end
        
        
        if exist('dz','var') == 1
            cdz(nn) = nanmedian(dz((bid==(nn*2) | bid==(nn*2)-1) & isnan(dcls')==0));
        end
        if exist('rootsemi','var') == 1;
                crootsemi(nn) = nanmedian(rootsemi((bid==(nn*2) | bid==(nn*2)-1) & isnan(dcls')==0));
        end
    else
        cz(nn) = NaN;
    end
end

%**********END OF AGGREGATE POINT-BY-POINT MEASUREMENTS TO FIND CHANNEL FIRST MOMENTS ***************************

    scatter(cx(s(10).cwrange>5),cy(s(10).cwrange>5),1500,'*')


%Find planar elevation (from f1_f2_contact_second_pass_b20g02_only_plane_fit_one_point_per_contact_pick_24feb14_rms_output ...)    
%of minx, miny

cetl = 6026.681 + minx*(0.00648907) + miny*(-0.0051384);%contact elevation at top left (of transect)

%Kludge fix for mixed preservation styles
cpstyle(cpstyle==9) = 8;
cpstyle(cpstyle==10) = 8;

figure('units','normalized','outerposition',[0 0 1 1])
    errorbar(cz,s(10).cw,s(10).cwdev,'xk','LineWidth',2)
    hold on;grid on; colorbar
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz,s(10).cw,300,cx,'LineWidth',2)
                saveas(gca,'DebugPlanarFit24FEB2014_1.eps','epsc');
    
figure('units','normalized','outerposition',[0 0 1 1])
    errorbar(cz ,s(11).cw,s(11).cwdev,'xk','LineWidth',2)
    hold on;grid on
    set(gca,'FontSize',18)
    title('Raw elevation, color is log2(quality)')
    xlabel('Raw elevation'); ylabel('Channel width')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz ,s(10).cw,300,log(cqual)./log(2),'LineWidth',2)
                    saveas(gca,'DebugPlanarFit24FEB2014_2.eps','epsc');


figure('units','normalized','outerposition',[0 0 1 1])
    errorbar(cz ,s(11).cw,s(11).cwdev,'xk','LineWidth',2)
    hold on;colorbar;grid on
    set(gca,'FontSize',18)
    title('Raw elevation, color is log2(pstyle)')
    xlabel('Raw elevation'); ylabel('Channel width')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz ,s(10).cw,300,log(cpstyle)./log(2),'LineWidth',2)    
                    saveas(gca,'DebugPlanarFit24FEB2014_3.eps','epsc');

    
    
figure('units','normalized','outerposition',[0 0 1 1])
    errorbar(cz - cetl - cx.*0.00648907 + cy.*0.0051384,s(11).cw,s(11).cwdev,'xk','LineWidth',2)
    hold on;grid on;colorbar
    czcorr = cz - cetl - cx.*0.00648907 + cy.*0.0051384;
    set(gca,'FontSize',18)
    title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(quality)'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz - cetl - cx.*0.00648907 + cy.*0.0051384,s(11).cw,300,log(cqual)./log(2),'LineWidth',2)
                    saveas(gca,'DebugPlanarFit24FEB2014_4.eps','epsc');


figure('units','normalized','outerposition',[0 0 1 1])
    errorbar(cz - cetl - cx.*0.00648907 + cy.*0.0051384,s(11).cw,s(11).cwdev,'xk','LineWidth',2)
    hold on;grid on;colorbar
    czcorr = cz - cetl - cx.*0.00648907 + cy.*0.0051384;
    set(gca,'FontSize',18)
    title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(pstyle)'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz - cetl - cx.*0.00648907 + cy.*0.0051384,s(11).cw,300,log(cpstyle)./log(2),'LineWidth',2)    
                    saveas(gca,'DebugPlanarFit24FEB2014_5.eps','epsc');

figure('units','normalized','outerposition',[0 0 1 1])
    errorbar(cz(cqual<=2) - cetl - cx(cqual<=2).*0.00648907 + cy(cqual<=2).*0.0051384, ...
             s(11).cw(cqual<=2),s(11).cwdev(cqual<=2),'xk','LineWidth',2)
    hold on;grid on;colorbar
    czcorr = cz(cqual<=2) - cetl - cx(cqual<=2).*0.00648907 + cy(cqual<=2).*0.0051384;
    set(gca,'FontSize',18)
    title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(quality)'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz(cqual<=2) - cetl - cx(cqual<=2).*0.00648907 + cy(cqual<=2).*0.0051384, ...
            s(11).cw(cqual<=2),300,log(cqual(cqual<=2))./log(2),'LineWidth',2)
                    saveas(gca,'DebugPlanarFit24FEB2014_4_PQual2or1.eps','epsc');


figure('units','normalized','outerposition',[0 0 1 1])
    errorbar(cz(cqual<=2) - cetl - cx(cqual<=2).*0.00648907 + cy(cqual<=2).*0.0051384, ...
             s(11).cw(cqual<=2),s(11).cwdev(cqual<=2),'xk','LineWidth',2)
    hold on;grid on;colorbar
    czcorr = cz(cqual<=2) - cetl - cx(cqual<=2).*0.00648907 + cy(cqual<=2).*0.0051384;
    set(gca,'FontSize',18)
    title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(pstyle)'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz(cqual<=2) - cetl - cx(cqual<=2).*0.00648907 + cy(cqual<=2).*0.0051384, ...
            s(11).cw(cqual<=2),300,log(cpstyle(cqual<=2))./log(2),'LineWidth',2)    
                    saveas(gca,'DebugPlanarFit24FEB2014_5_PQual2or1.eps','epsc');  
                    
                    
figure

    errorbar(cdz,s(10).cw,cwdev,'xk','LineWidth',2)
    hold on
    scatter(cdz,s(10).cw,50,'co','fill')
    scatter(cdz(cdz<0),         s(10).cw(cdz<0),50,'bo','fill')
    scatter(cdz(cdz>0 & cdz<50),s(10).cw(cdz>0 & cdz<50),50,'wo','fill')
    scatter(cdz(cdz>0 & cdz<50),s(10).cw(cdz>0 & cdz<50),50,'ko')

    %errorbar(cdz(cwhichdtm==2),cw(cwhichdtm==2),cwdev(cwhichdtm==2),'xr')

    %herrorbar(cdz,cw,crootsemi,'x')
xlabel('Relative Time (Height Above F1/F2 Contact (m))','FontSize',20)
ylabel('River Channel Width (m)','FontSize',20)
set(gca,'FontSize',20)
set(gca,'xminortick','on')
set(gca,'yminortick','on')
ylim([0 70])
xlim([-100 150])
%grid on
box on
line([-100,0],[29.3,29.3],'Color','b','LineWidth',2)
line([50,150],[15.8616,15.8616],'Color','c','LineWidth',2)
line([0,50],[ 26.8001, 26.8001],'Color','k','LineWidth',2)

%********** SAVE OUTPUT ***************************
switch output_type
    case 'terse'
        save('parseTransect1ChannelWidthMeasurementsSecondPassForLPSC2014_v2_output_terse', ...
             's','cqual','cpstyle','cx','cy','cz','minx','miny')
    case 'complete'
    case 'none'
end

%********** END OF SAVE OUTPUT ***************************

%error('Not doing probability summation while debugging')

%Quick summation of probability distribution (assuming uncorrelated, Guassian error
%bars) From Wikipedia, "Multivariate normal distribution"
p = zeros(601,401);
warning('Assuming fixed sigmax for debug')
for nn = 1:(max(bid)/2) %add the probability kernel of each observation
oldp = p;
nn
            sigmax = 29.0847;%<-- Assuming fixed sigma x for debugging %crootsemi(nn);
            %The value of fixed sigma x comes from
            %f1_f2_contact_second_pass_b20g02_only_plane_fit_one_point_per_contact_pick_24feb14_rms_output
            
            sigmay = s(11).cwdev(nn);
            normfactor = 1/(2*pi*sigmax*sigmay);
            mux = cz(nn) - cetl - cx(nn).*0.00648907 + cy(nn).*0.0051384;
            muy = s(11).cw(nn);
            
            if sum(isnan([sigmax,sigmay,mux,muy]))==0
                for x = -300:1:300 %height above F1/F2 contact (m)
                    for y = -100:1:300 %river channel width (m)

                            p(x+301,y+101) = p(x+301,y+101) + ...
                            normfactor*exp(-1/2*(((x-mux)^2/(sigmax^2))+((y-muy)^2/(sigmay^2))));
                    end
                end
            end
                %For debug only
                %check that the additional probability sums to 1 (note that
                %the discretisation and the cutoffs will introduce errors!)
                %if abs((sum(p(:)) - sum(oldp(:)))-1) > 0.01
                %    error('MVN Dist Checksum Error!')
                %end
 end

 [X,Y] = meshgrid(-300:1:300,-100:1:300);
%contour(X,Y,p');grid on %show the sum of all kernels

%show 1 sigma and 2 sigma error windows for channel width as a function of
%stratigraphic position

pcs = cumsum(p,2);
ps = sum(p');
for l = 1:601 %cumsum normalized
pcsn(l,:) = pcs(l,:)./ps(l);
end

hold on
contour(X,Y,pcsn',[0.05,(0.31731/2),1-((0.31731/2)),0.95]);

%find pooled variance for 50m stratigaphic-position bins: (Note to self:
%should preferably use *fractional* variance?)
stratbin = 1;
for stratz = -200:50:200
    to_use(stratbin).v = find(cdz>stratz & cdz <stratz+50);
    ninbin(stratbin) = length(to_use(stratbin).v);
    pooledvar(stratbin) = sqrt(sum(css(to_use(stratbin).v).*(cwdev(to_use(stratbin).v).^2))/sum(css(to_use(stratbin).v)));
    stratbin = stratbin + 1;
end

%Overall pooled variance:
opooledvar = sqrt(sum(css.*(cwdev.^2))/sum(css))