%Edwin Kite
%Purpose: Take as input hand-picked
%sinuous paleo-channel centerlines (three picks per centerline)
%and output aggregate half-meander
%wavelength and sinuosity information.
% - using the multiple picks as an internal quality check
%and imposing some smoothing and a minimum sinuosity.

%******************* Set Control Parameters *******************************
dbip = 10; %distance between interpolated points (m) for inflection-point determination
dbip_hmp = 5; %distance between interpolated points (m) for half-meander pairing
smoothn = 10; %number of interpolated points to smooth curvature for (to prevent artifactual inflection pks)

%See rivers notebook entries for 18FEB2014 for rationale on the following two parameter
%choices.
minsinu = 1.1; %minimum half-meander sinuosity to be considered as a half-meander
max_dcls_fr = 0.25; %maximum fractional median distance between half-meander straight lines to be considered replicable
importpicks = 'all';
assert(strcmp(importpicks,'all')==1)
%******************* Set Graphics & Housekeeping Parameters **************************
EdwinColor = {'k','b','g','y','r',[.5 .6 .7],[.8 .2 .6]};
savedebugfigures = 0;
picturesque = 1; %1 to show graphics. 0 to not show graphics.
output_type='none';
%******************* Import Data **********************************
load_meander_wavelengths
%********* Pepper Points Uniformly Along Channel Centerlines **************

for pk = 1:3
    d(pk).ctl(1:max(d(pk).cid)) = 0; %channel total length
    %Interpolate uniformly in s, dir coordinates
    for l = 1:(size(d(pk).vx,1)-1)
        if d(pk).cid(l) == d(pk).cid(l+1)%Is it from the same channel?
            d(pk).lseg(l)     =   sqrt((d(pk).vx(l+1) - d(pk).vx(l)).^2 + ...
                                       (d(pk).vy(l+1) - d(pk).vy(l)).^2); %length of segment
            d(pk).ctl(d(pk).cid(l)) = d(pk).ctl(d(pk).cid(l)) + d(pk).lseg(l); %channel total length
            d(pk).lac(l+1)    = d(pk).ctl(d(pk).cid(l)); %length along channel for this vertex
        else %last point on channel trace
            %lac(l) = ctl(cid(l));
        end
    end
end

for pk = 1:3
    for m = 1:max(d(pk).cid) 
        %position vector of interpolated points
        d(pk).p(m).x = interp1(d(pk).lac(d(pk).cid==m),d(pk).vx(d(pk).cid==m),[0:dbip:d(pk).ctl(m)]);
        d(pk).p(m).y = interp1(d(pk).lac(d(pk).cid==m),d(pk).vy(d(pk).cid==m),[0:dbip:d(pk).ctl(m)]);
        d(pk).p(m).z = interp1(d(pk).lac(d(pk).cid==m),d(pk).vz(d(pk).cid==m),[0:dbip:d(pk).ctl(m)]);
        d(pk).p(m).dir = [    atan2d(d(pk).p(m).y(2:end) - d(pk).p(m).y(1:end-1), ...
                                       d(pk).p(m).x(2:end) - d(pk).p(m).x(1:end-1)) NaN]; %returns direction in degrees anticlockwise from E. Not defined for last point.
        d(pk).p(m).th =  [NaN d(pk).p(m).dir(2:end) - d(pk).p(m).dir(1:end-1) ]; %change in direction (degrees). Not defined for first or last points.
        d(pk).p(m).th(d(pk).p(m).th<-180) =  (360 + d(pk).p(m).th(d(pk).p(m).th<-180));
        d(pk).p(m).th(d(pk).p(m).th> 180) = -(360 - d(pk).p(m).th(d(pk).p(m).th> 180));
        d(pk).p(m).smooth = smooth(d(pk).p(m).th,smoothn)';
        d(pk).p(m).smooth([1,end]) = NaN; %smooth incorrectly smooths over NaNs
        
        d(pk).p(m).th1plusth2 = [NaN  d(pk).p(m).smooth(1:end-1) + d(pk).p(m).smooth(2:end)    ];
        d(pk).p(m).th2plusth3 = [     d(pk).p(m).smooth(1:end-1) + d(pk).p(m).smooth(2:end) NaN];

        %Following Hemberger, 1986 Master's thesis, U. Virginia.
        d(pk).p(m).isinfl = [1 - (sign(d(pk).p(m).th1plusth2) == sign(d(pk).p(m).th2plusth3))];
        %Using this definition of inflection point, an inflection point cannot occur before the 3rd interpolated point in the channel.
        d(pk).p(m).isinfl((isnan(d(pk).p(m).th1plusth2 + d(pk).p(m).th2plusth3))>0) = 0;
    end
end

if picturesque==1
hold on

for pk = 1:3
    scatter(d(pk).vx,d(pk).vy,50,'kx') %rem(a(:,10)*1e9,11)

    for l = 1:size(d(pk).p,2)
        scatter(d(pk).p(l).x,d(pk).p(l).y,30,d(pk).p(l).isinfl*3+pk,'Filled')
    end
    %scatter([1:length(d(pk).p(1).th)],d(pk).p(1).th,100,d(pk).p(1).isinfl,'Filled')
end

hold on; title('debugging plot'); grid on
end
%plot(p(1).th)

%******************* Define Half-Meanders **************************
for pk = 1:3
    for m = 1:max(d(pk).cid)
        d(pk).p(m).okhmep = d(pk).p(m).isinfl; %O.K. half-meander end points
        assert(d(pk).p(m).okhmep(1)   == 0, 'inflections must not occur at ends of channel traces')
        assert(d(pk).p(m).okhmep(end) == 0, 'inflections must not occur at ends of channel traces')
        d(pk).p(m).hmeplist = find(d(pk).p(m).okhmep==1);
        if length(d(pk).p(m).hmeplist)>1
            d(pk).h(m).halfmeanderspresent = 1;
        for epi = 1:(length(d(pk).p(m).hmeplist)-1) %end-point index 
            ep = d(pk).p(m).hmeplist(epi);
%            d(pk).h(m).iip(epi) = [d(pk).p(m).hmeplist(epi):d(pk).p(m).hmeplist(epi+1)];
            d(pk).h(m).p(epi,1,1:2) = ...
                [d(pk).p(m).x(d(pk).p(m).hmeplist(epi)) d(pk).p(m).x(d(pk).p(m).hmeplist(epi+1))];%x,start|end
            d(pk).h(m).p(epi,2,1:2) = ...
                [d(pk).p(m).y(d(pk).p(m).hmeplist(epi)) d(pk).p(m).y(d(pk).p(m).hmeplist(epi+1))];%y,start|end
            d(pk).h(m).p(epi,3,1:2) = ...
                [d(pk).p(m).z(d(pk).p(m).hmeplist(epi)) d(pk).p(m).z(d(pk).p(m).hmeplist(epi+1))];%z,start|end
            d(pk).h(m).ls(epi) = ... %half-meander straight-line distance
                sqrt( ...
                (d(pk).h(m).p(epi,2,2) - d(pk).h(m).p(epi,2,1))^2 ...
               +(d(pk).h(m).p(epi,1,2) - d(pk).h(m).p(epi,1,1))^2);
            d(pk).h(m).lc(epi) = (d(pk).p(m).hmeplist(epi+1) - d(pk).p(m).hmeplist(epi))*dbip;%along-curve distance
            d(pk).h(m).th(epi) = nanmean(d(pk).p(m).th(d(pk).p(m).hmeplist(epi):d(pk).p(m).hmeplist(epi+1))); %for debugging
            d(pk).h(m).sinu(epi) = d(pk).h(m).lc(epi) ./ d(pk).h(m).ls(epi);
        end

        %Half-meanders below a sinuosity threshold 'minsinu' will be rejected. 
        %Remove end-points
        %that are the end-points for TWO half meanders below the sinuosity
        %threshhold from the end-point index 
        %end-point index. 
        %Then recalculate half-meander sinuosity between the remaining
        %end-points.

        %This is done removing one end-point at a time in order of minimum
        %"average" sinuosity

        half_meander_clipping_is_complete = 0;

        if m == 25 %debugging
            m = m;
        end

        while half_meander_clipping_is_complete == 0
            d(pk).oksinu = d(pk).h(m).sinu>minsinu;

        d(pk).h(m).toberemoved = ...
            1+find((d(pk).oksinu(1:end-1) + d(pk).oksinu(2:end))<1); %returns indices. "1+" because the first half-meander links end-points 1 & 2.
        if size(d(pk).h(m).sinu,2) >= 2
        [null d(pk).min_mean_sinu(m)] = min(0.5.*(d(pk).h(m).sinu(1:end-1) + d(pk).h(m).sinu(2:end)));
        %"1+" because 1st half-meander links end-points 1 & 2:
        d(pk).min_mean_sinu(m) = d(pk).min_mean_sinu(m) + 1; 
        if sum(d(pk).min_mean_sinu(m)==d(pk).h(m).toberemoved) == 1 %minimum mean sinuosity is for an endpoint bracketed by two half-meanders with less than critical sinuosity?
                d(pk).p(m).hmeplist(d(pk).min_mean_sinu(m)) = []; %half-meander end-point list

            d(pk).h(m).p = []; d(pk).h(m).ls = []; d(pk).h(m).lc = []; d(pk).h(m).sinu = [];

            for epi = 1:(length(d(pk).p(m).hmeplist)-1) %end-point index 
                ep = d(pk).p(m).hmeplist(epi);
%                d(pk).h(m).iip(epi) = [d(pk).p(m).hmeplist(epi):d(pk).p(m).hmeplist(epi+1)]; %indices of interpolated points between end-points
                d(pk).h(m).p(epi,1,1:2) = ...
                    [d(pk).p(m).x(d(pk).p(m).hmeplist(epi)) d(pk).p(m).x(d(pk).p(m).hmeplist(epi+1))];%x,start|end
                d(pk).h(m).p(epi,2,1:2) = ...
                    [d(pk).p(m).y(d(pk).p(m).hmeplist(epi)) d(pk).p(m).y(d(pk).p(m).hmeplist(epi+1))];%y,start|end
                d(pk).h(m).p(epi,3,1:2) = ...
                    [d(pk).p(m).z(d(pk).p(m).hmeplist(epi)) d(pk).p(m).z(d(pk).p(m).hmeplist(epi+1))];%z,start|end
                d(pk).h(m).ls(epi) = ... %half-meander straight-line distance
                    sqrt( ...
                    (d(pk).h(m).p(epi,2,2) - d(pk).h(m).p(epi,2,1))^2 ...
                   +(d(pk).h(m).p(epi,1,2) - d(pk).h(m).p(epi,1,1))^2);
                d(pk).h(m).lc(epi) = (d(pk).p(m).hmeplist(epi+1) - d(pk).p(m).hmeplist(epi))*dbip;%along-curve distance
                d(pk).h(m).sinu(epi) = d(pk).h(m).lc(epi) ./ d(pk).h(m).ls(epi);
            end 
        else %Remaining half-meanders exceed minimum sinuosity
            half_meander_clipping_is_complete = 1;
        end
        else %Only one half-meander (of any sinuosity) remains.
            half_meander_clipping_is_complete = 1;
        end

        end

        else %One or zero half-meander end-points for this pick
            d(pk).h(m).halfmeanderspresent = -1;
        end
    end
end

%****** Find "Channel Pairs" so that they can be displayed together ***************
cp(:,1) = [1:max(d(pk).cid)]; 
for pick = 2:3
    for m =1:max(d(1).cid) %First Pick channel ID
        m
        %Inexact, but should work except for VERY close channels
        lov = find(d(1).cid==m); %list of vertices
        for n = 1:max(d(1).cid)%Find the "typical" distance to the vertices of another channel
                lovpp = find(d(pick).cid==n); %list of vertices of potential pair
                cl2 = zeros(1,length(lov));
            for q = 1:length(lov); %First Pick channel vertices
                
             cl2(q) = min((d(1).vx(lov(q)) - d(pick).vx(lovpp)).^2 + ...
                          (d(1).vy(lov(q)) - d(pick).vy(lovpp)).^2); %distance-squared to potential-pair vertices
                %closest distance-squared to vertex with this CID

%                 d2tppv = ((d(1).vx(lov(q)) - d(pick).vx(lovpp)).^2 + ...
%                           (d(1).vy(lov(q)) - d(pick).vy(lovpp)).^2); %distance-squared to potential-pair vertices
%                 cl2(q) = min(d2tppv);%closest distance-squared to vertex with this CID
            end
            cpd(pick,m,n) = mean(sqrt(cl2));
        end
    end
    [null_dump cp(:,pick)] = min(squeeze(cpd(pick,:,:))');
end
    clear lov lovpp

%****************** Compute Summary Parameters & Plot Results for Debugging *********************
%figure
    mi = zeros(size(d)); mreji = zeros(size(d));
for m = 1:max(d(pk).cid) %channel-centerline number
    
    if picturesque==1;figure;end
    for pk = 1:3 %each channel-centerline is pk'ed 3 times
        if picturesque==1
                 scatter(d(pk).p(cp(m,pk)).x,d(pk).p(cp(m,pk)).y,20,d(pk).p(cp(m,pk)).z,'Filled')
         hold on        
                 scatter(d(pk).p(cp(m,pk)).x(find(d(pk).p(cp(m,pk)).isinfl)), ...
                         d(pk).p(cp(m,pk)).y(find(d(pk).p(cp(m,pk)).isinfl)),200,'k*')
        end
    %            scatter(p(m).x,p(m).y,100,p(m).isinfl,'Filled')
                hold on
        if  d(pk).h(cp(m,pk)).halfmeanderspresent ~= -1;
        for n = 1:size(d(pk).h(cp(m,pk)).sinu,2)
            if d(pk).h(cp(m,pk)).sinu(n) > minsinu
                if picturesque==1
                line([d(pk).h(cp(m,pk)).p(n,1,1) d(pk).h(cp(m,pk)).p(n,1,2)], ...
                     [d(pk).h(cp(m,pk)).p(n,2,1) d(pk).h(cp(m,pk)).p(n,2,2)],'Color',EdwinColor{pk})
                end
                       mi(pk)    = mi(pk) + 1;
                mls(pk,mi(pk))   = d(pk).h(cp(m,pk)).ls(n);         
                 mx(pk,mi(pk),:) = d(pk).h(cp(m,pk)).p(n,1,:);
                 my(pk,mi(pk),:) = d(pk).h(cp(m,pk)).p(n,2,:);
                 mz(pk,mi(pk),:) = d(pk).h(cp(m,pk)).p(n,3,:);
                 mm(pk,mi(pk))   = m;  %centerline ID for half-meander
                mpk(pk,mi(pk))   = pk; %indexing information useful for pairing up good half-meanders
                 mn(pk,mi(pk))   = n;  %indexing information useful for pairing up good half-meanders
              msinu(pk,mi(pk))   = d(pk).h(cp(m,pk)).sinu(n);
              mrhmi(pk,mi(pk))   = 0;  %Replicable meander index (to be assigned later)
            else
                if picturesque==1
                line([d(pk).h(cp(m,pk)).p(n,1,1) d(pk).h(cp(m,pk)).p(n,1,2)], ...
                     [d(pk).h(cp(m,pk)).p(n,2,1) d(pk).h(cp(m,pk)).p(n,2,2)],'Color','r')    
                end
                mreji(pk) = mreji(pk) + 1;
                mrejsinu(pk,mreji(pk)) = d(pk).h(cp(m,pk)).sinu(n); 
            end
        end
        end
    
    end
    if picturesque==1
    for pk = 1:3 %plot text labels on top
        for n = 1:size(d(pk).h(cp(m,pk)).sinu,2)
                        if d(pk).h(cp(m,pk)).sinu(n) > minsinu
            text(0.5*(d(pk).h(cp(m,pk)).p(n,1,1)+d(pk).h(cp(m,pk)).p(n,1,2)), ...
                         0.5*(d(pk).h(cp(m,pk)).p(n,2,1)+d(pk).h(cp(m,pk)).p(n,2,2)), ...
                          num2str(round(d(pk).h(cp(m,pk)).ls(n))))
                        end
        end
    end

    
            title({'Centerlines Picked 3 Times, Then Uniformly Interpolated;'; ...
                   'Results Smoothed; then Half-Meanders IDd using Inflection Points (*)'; ...
                   'Half-Meanders with Sinuosity <1.1 not measured (red).';...
                   strcat('Channel Centerline ID:',num2str(cp(m,1)))})
            %title(strcat('Pass 1:',num2str(cp(m,1)),',Pass 2:',num2str(cp(m,2)),',Pass 3:',num2str(cp(m,3))))
            axis equal;box on;xlabel('distance x (m)');ylabel('distance y (m)')
            colorbar
    end
    if savedebugfigures == 1 
            %save to eps color file
            saveas(gca, ...
            strcat('DebugParseTransect1Centerline19FEB2014v2',num2str(cp(m,1)),'.eps'), ...
            'epsc');
    end

end

meanx = mean(mx,3); meany = mean(my,3); meanz = nanmean(mz,3); 


% if picturesque==1
%     figure
%     scatter(mz(:),mls(:),30,meanx(:),'Filled');grid on
%     title('Half-meander straight-line distance versus elevation')
% 
%     figure
%     hist([mrejsinu(:)' msinu(:)'],[1:0.01:4])
%     title('Histogram of all sinuosities')
% end

%************************* Identify Replicable Half-Meanders By Pairing-Up Single Picks *************************
rhmi = 0;%replicable meander index
dcls_fr = zeros(prod(size(mls))); %square matrix for fractional distances between half-meanders
for m = 1:max(d(pk).cid)%For each channel ID
    m
        %figure; axis equal
    
        ghmlist = find(mm(:)==m); %from all picks
      if isempty(ghmlist) == 0
        assert(sum(mn(ghmlist)==0)==0,'Good half-meander assignment error')
        assert(sum(mpk(ghmlist)==0)==0,'Good half-meander assignment error')
                
      for ghmi = 1:length(ghmlist) %For each good half-meander assoc. with this channel ID
            picklist = [1 2 3];

            ghm = ghmlist(ghmi); %good half-meander, from all picks
            picklist(picklist == mpk(ghm)) = []; %No self-pairing
            assert(length(picklist)==2,'Code not tested for no. of picks not equal to 3')
            
            %Pepper points uniformly along half-meander straight-line
            %length
            hmx1 = mx(:,:,1); hmx2 = mx(:,:,2); hmy1 = my(:,:,1); hmy2 = my(:,:,2);
            hmslx = interp1([0 mls(ghm)],...
                            [hmx1(ghm) hmx2(ghm)], ...
                            [0:dbip_hmp:mls(ghm)]);
            hmsly = interp1([0 mls(ghm)],...
                            [hmy1(ghm) hmy2(ghm)], ...
                            [0:dbip_hmp:mls(ghm)]);          
            
            for pfpi = 1:length(picklist) %pick for pairing    
                pfp = picklist(pfpi);
            %Are there any good half-meanders that one could pair with?   
            if sum(mpk(ghmlist)==pfp)>=1
                cplist = find(mm(:)==m & mpk(:)==pfp); %candidate pairs list
            %Look for the closest (ok-sinuosity) half-meander in 
            %units of fractional length
            
%///////%start of copy from width code
            for w = 1:length(cplist) %2:(size(px,1)-1) 
                       wi = cplist(w); %index in list of all half-meanders
                        x1 = hmx1(wi); y1 = hmy1(wi);
                        x2 = hmx2(wi); y2 = hmy2(wi);
                    clear x3 y3 xc yc
            for n = 1:length(hmslx)     
                    x3(n) = hmslx(n); y3(n) = hmsly(n);

                lambda(n) = ( (x3(n)-x1).*(x2-x1) + (y3(n)-y1).*(y2-y1) ) ./ ...
                             ((x2-x1).^2 + (y2-y1).^2);

                if lambda(n)<0
                    xc(n) = x1; yc(n) = y1;
                    lambda(n) = 0;
                elseif lambda(n)>1
                    xc(n) = x2; yc(n) = y2;
                    lambda(n) = 1;
                else
                    xc(n) = x1 + lambda(n)*(x2-x1);
                    yc(n) = y1 + lambda(n)*(y2-y1);
                end


            end
                clear dcls_temp
                dcls_temp = sqrt((x3 - xc).^2 + (y3 - yc).^2);
                %assert(ghm~=wi,'Indexing error')
                dcls(ghm,wi) = nanmedian(dcls_temp); % Zero values of xc are used to set dcls_temp(m) to NaN
             dcls_fr(ghm,wi) = dcls(ghm,wi)./mls(ghm);

             %debug plot             
%                                  line([x3(1) x3(end)],[y3(1) y3(end)],'Color',EdwinColor{rem(ghm,7)+1})
%                                  hold on
%                                  line([mean([x1 x2]) mean(x3)],[mean([y1 y2]) mean(y3)],'LineWidth',2,'Color',EdwinColor{rem(ghm,7)+1})
%                                  text(mean([x1 x2]),mean([y1 y2]), ...
%                                  num2str(dcls_fr(ghm,wi)),'FontSize',25,'Color',EdwinColor{rem(ghm,7)+1})
%                                  
%                                  
%                                  title(strcat('ghm:',num2str(ghm),',wi:',num2str(wi)))
%                                  scatter(hmslx,hmsly,30,EdwinColor{rem(ghm,7)+1})
             %for n = 1:length(hmslx)
             %    line([hmslx(n),xc(n)],[hmsly(n),yc(n)],'Color','r')
             %end
             %text(hmslx(round(length(hmslx)/2)),hmsly(round(length(hmsly)/2)), ...
             %     num2str(dcls_fr(ghm,wi)),'FontSize',25)
             
            end    
        %/////// % end of copy from width code    
            end

            end
      end
         
% Assign pairs if min distance is acceptable for both direct and reciprocal fractional distances
for ghmi = 1:length(ghmlist)
    ghm = ghmlist(ghmi);
    if m == 2
        m
    end
    for w = 1:length(ghmlist)
        wi = ghmlist(w);    
        if (dcls_fr(ghm,wi) < max_dcls_fr) & (dcls_fr(wi,ghm) < max_dcls_fr) & ...
           (dcls_fr(ghm,wi) > 0) &           (dcls_fr(wi,ghm) > 0)
           %OK pairing
              if mrhmi(ghm) > 0
                 mrhmi(wi) = mrhmi(ghm);
              elseif mrhmi(wi) > 0
                 mrhmi(ghm) = mrhmi(wi);
              else %New replicable half-meander
                 rhmi = rhmi + 1;%replicable half-meander index
                 mrhmi(ghm) =  rhmi;
                 mrhmi(wi)  =  rhmi;
              end 
        end
    end
end
% Display pairings showing the pick number
% and the distance
      end
end

%******************** Retain Only Replicable Half-Meanders ******************************
for l = 1:rhmi
    rhm(l).quality =   sum(mrhmi(:)==l); %Number of picks contributing to meander wavelength 
    rhm(l).m       = mean( mm(mrhmi==l));
    assert(range(mm(mrhmi==l))==0,'Error in channel assignments!')
    rhm(l).lambda  =    mean(mls(mrhmi==l));
    rhm(l).err     =     std(mls(mrhmi==l));
    rhm(l).x       =    mean( meanx(mrhmi==l));
    rhm(l).y       =    mean( meany(mrhmi==l));
    rhm(l).z       = nanmean( meanz(mrhmi==l));
    rhm(l).hmx1    =    mean(  hmx1(mrhmi==l));
    rhm(l).hmx2    =    mean(  hmx2(mrhmi==l));
    rhm(l).hmy1    =    mean(  hmy1(mrhmi==l));
    rhm(l).hmy2    =    mean(  hmy2(mrhmi==l));
    
    for pki = 1:3
        rhm(l).pk(pki).n = find(mrhmi(:)==l & mpk(:) == pki); %Can be empty
        if length(rhm(l).pk(pki).n) > 1
            disp(strcat('More than 1 match for the following:', ...
                        'reproducible half-meander:',num2str(l), ...
                        ', pick:',       num2str(pki), ...
                        ', channel ID:', num2str(rhm(l).m)))
            rhm(l).pk(pki).n = rhm(l).pk(pki).n(1);
        end
    end
end

%******************** Aggregate Replicable Half-Meanders on A Single Channel ******************************

%Combine the MEAN and STANDARD DEVIATION of replicable half-meanders to
%obtain a characteristic half-meander wavelength for
%each picked channel.

for m = 1:length(d(1).p)%For each channel ID
    %Based on Wikipedia article
    rhmtc = find([rhm(:).m] == m);%replicable half-meanders, this channel
    
    if isempty(rhmtc)==0
    %'Standard_deviation#Combining_standard_deviations' as of 24FEB2014
    %Non-overlapping populations

    chm(m).w = sum([rhm(rhmtc).quality].*[rhm(rhmtc).lambda]) ...
            ./sum([rhm(rhmtc).quality]);
    chm(m).err = sqrt( ( sum([rhm(rhmtc).quality].*( [rhm(rhmtc).err].^2 + [rhm(rhmtc).lambda].^2 )) ...
                    / sum([rhm(rhmtc).quality]))...
                    - chm(m).w.^2);
    else
        chm(m).w   = NaN;
        chm(m).err = NaN;
    end
end

%Debug plot
figure
    errorbar([1:length(d(1).p)],[chm.w],[chm.err],'r.')
    hold on
for m = 1:length(d(1).p)
    rhmtc = find([rhm(:).m] == m);
    errorbar(0.25+m.*ones(size(rhmtc)),[rhm(rhmtc).lambda],[rhm(rhmtc).err],'*');
end

%******************** Assign All Interpolated Points to the Closest Replicable Half-Meander  ******************************
%This is needed because some channels are long (with doubt about
%stratigraphic continuity), so assigning interpolated points to the
%*channel-averaged* wavelength may not be appropriate.


%To enable width-wavelength comparison for "bracketed" half-meanders,
%Assign all interpolated points from all picks either to the CLOSEST 
%replicable half-meander, or NOTHING RELEVANT (-1).
%figure;scatter(cqual(lambda>0)',2.*lambda(lambda>0)'./w(lambda>0)',100,w(lambda>0)','Filled')
for pki = 1:3
    for m = 1:max(size(d(pki).p))
            d(pki).p(cp(m,pki)).arhm = []; %Reset
        if sum([rhm(:).m] == m) == 0;
            %associated reproducible half-meander
            d(pki).p(cp(m,pki)).arhm = -1 + zeros(size(d(pki).p(cp(m,pki)).x)); %-1: none associated.
        elseif sum([rhm(:).m] == m) == 1; %only one associated reproducible half-meander
            d(pki).p(cp(m,pki)).arhm = find([rhm(:).m] == m) + zeros(size(d(pki).p(cp(m,pki)).x));
        elseif sum([rhm(:).m] == m) >  1; 
            cahmlist = find([rhm(:).m] == m); %candidate associated half-meander list
            %There is more than one reproducible half-meander on this
            %channel segment, so 
            %associate each interpolated point with the closest
            %half-meander (measured along-channel-thread).
                    d(pki).p(cp(m,pki)).arhm = zeros(size(d(pki).p(cp(m,pki)).x));
            for q = 1:length(cahmlist)
                    clear dq dsp dep
            for n = 1:length(d(pki).p(cp(m,pki)).x)
                    dsp(n,q) = (rhm(cahmlist(q)).hmx1 - d(pki).p(cp(m,pki)).x(n)).^2 ...
                             + (rhm(cahmlist(q)).hmy1 - d(pki).p(cp(m,pki)).y(n)).^2; %Distance to start point
                    dep(n,q) = (rhm(cahmlist(q)).hmx2 - d(pki).p(cp(m,pki)).x(n)).^2 ...
                             + (rhm(cahmlist(q)).hmy2 - d(pki).p(cp(m,pki)).y(n)).^2; %Distance to end point
                    dcp(n,q) = (rhm(cahmlist(q)).x - d(pki).p(cp(m,pki)).x(n)).^2 ...
                             + (rhm(cahmlist(q)).y - d(pki).p(cp(m,pki)).y(n)).^2;
            end           
                %[null I(q)] = min(sqrt(dsp(:,q))+sqrt(dep(:,q))); %closest 'n' to this half-meander
                 [null I(q)] = min(sqrt(dcp(:,q)));
                 [null Is(q)]= min(sqrt(dsp(:,q)));
                 [null Ie(q)]= min(sqrt(dep(:,q)));
                 d(pki).p(cp(m,pki)).arhm(min([Is(q),Ie(q)]):max([Ie(q),Is(q)])) = cahmlist(q);

            end
            %Assign to nearest half-meander:
                already_assigned = find(d(pki).p(cp(m,pki)).arhm > 0);
            for n = 1:length(d(pki).p(cp(m,pki)).x)
                [C2 I2] = min(abs(already_assigned - n));
                if C2 > 0
                d(pki).p(cp(m,pki)).arhm(n) = d(pki).p(cp(m,pki)).arhm(already_assigned(I2));
                end
            end
            
        end
    end
end

figure
for pki = 1:3
for m = 1:68
    if max(d(pki).p(cp(m,pki)).arhm) > -1 

scatter(d(pki).p(cp(m,pki)).x,d(pki).p(cp(m,pki)).y,10, ...
        d(pki).p(cp(m,pki)).arhm - min(d(pki).p(cp(m,pki)).arhm),'Filled')
%        scatter(d(pki).p(cp(m,pki)).x,d(pki).p(cp(m,pki)).y,10, ...
%        [rhm(d(pki).p(cp(m,pki)).arhm).lambda])
end
        hold on
        if pki==1
           text(mean(d(pki).p(cp(m,pki)).x),mean(d(pki).p(cp(m,pki)).y),num2str(cp(m,pki)))
        end
end
end
%End of assigning interpolated points to the closest replicable
%half-meander.



    for l = 1:rhmi
        mlist = find(mrhmi==l);
         for mindex = 1:length(mlist)
            m = mlist(mindex);
            line([hmx1(m),hmx2(m)],[hmy1(m),hmy2(m)]);
            hold on
         end
    end
    
    for m = 1:max(d(pk).cid)
        rhmtc = find([rhm(:).m] == m); %Replicable half-meanders, this channel
        %figure
        
        for l = [rhmtc]
        mlist = find(mrhmi==l);
         for mindex = 1:length(mlist)
            n = mlist(mindex);
            line([hmx1(n),hmx2(n)],[hmy1(n),hmy2(n)]);
            hold on
         end
        end
        
        axis equal; box on
        if picturesque==1
            for pk = 1:3
                line(d(pk).vx(d(pk).cid==cp(m,pk)),d(pk).vy(d(pk).cid==cp(m,pk)),'Color','k');
            end
            text(mean(d(1).vx(d(1).cid==m)),mean(d(1).vy(d(1).cid==m)),num2str(m));
            for l = [rhmtc]
                scatter([rhm(l).x],[rhm(l).y],250,[rhm(l).lambda],'Filled')
                   text([rhm(l).x],[rhm(l).y],strcat(num2str(rhm(l).lambda),'\pm',num2str(rhm(l).err)))
            end
            title(strcat('Channel Centerline ID:',num2str(cp(m,1))))
            saveas(gca, ...
            strcat('DebugParseTransect1CenterlinePairings',num2str(cp(m,1)),'.eps'), ...
            'epsc');
        end 
    end
    

%******************** Summarize and Plot Replicable Half-Meander Data ****************
if picturesque==1
    figure
    errorbar([rhm(:).z],[rhm(:).lambda],[rhm(:).err],'xk')
    hold on
    scatter([rhm(:).z],[rhm(:).lambda],30,[rhm(:).x],'Filled')
    box on;grid on; colorbar
    set(gca,'FontSize',18);xlabel('Corrected elevation');ylabel('Half-meander wavelength')
    title('Half-meander wavelength versus raw elevation')
    
    figure
    czadjusted = [rhm(:).z] - [rhm(:).x].*0.00648907 + [rhm(:).y].*0.0051384;
    errorbar(czadjusted,[rhm(:).lambda],[rhm(:).err],'xk')
    hold on
    scatter(czadjusted,[rhm(:).lambda],30,[rhm(:).x],'Filled')
    box on;grid on; colorbar
    title('Half-meander wavelength versus corrected elevation')
    set(gca,'FontSize',18);xlabel('Corrected elevation');ylabel('Half-meander wavelength')
    
    figure
    czadjusted = [rhm(:).z] - [rhm(:).x].*0.00648907 + [rhm(:).y].*0.0051384;
    scatter(czadjusted,[rhm(:).lambda],30,log10([rhm(:).quality])./log(2),'Filled')
    box on;grid on; colorbar
    set(gca,'FontSize',18);xlabel('Corrected elevation');ylabel('Half-meander wavelength')
    title('Half-meander wavelength versus corrected elevation: color is log2(quality)')

    
    figure
    scatter([rhm(:).x],[rhm(:).y],30,[rhm(:).lambda],'Filled')
    for l = 1:max(size(rhm))
        text([rhm(l).x],[rhm(l).y],num2str([rhm(l).m]));
    end
    box on;grid on; colorbar
end

%********** SAVE OUTPUT ***************************
switch output_type
    case 'terse'
        save('parseTransect1MeanderWavelengthMeasurements_output_terse','rhm','d','chm','minx','miny')
    case 'complete'
    case 'none'
end