%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 - but only for meanders that are 
%'reproducible' between picks -
%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 **********************************************
MeanderColor = {'k','b','g','y','r',[.5 .6 .7],[.8 .2 .6]};
savedebugfigures = 0;
picturesque = 2; %2 to show all graphics. 1 to show some graphics. 0 to not show graphics.
output_type='terse';
%******************* Import Data *************************************************************************
WhichTransect = 2; %A transect consists of meander-wavelength data gathered from 1 or more nearby DTMs.
load_meander_wavelengths %Currently set up to read ArcGIS 'shapefile' format.
%******************* Temporary Fix ***********************************************************************
if WhichTransect == 1
    disp('Applying Temporary Fix to Transect 1 x vertex coordinates!')
    d(1).vx = d(1).vx + 150;
end

%******************* Find "Channel Pairs" so that they can be displayed together **************************
%The pick-by-pick data in the ArcGIS shapefiles is not sorted. For each
%channel-centerline trace, we need to find the indices for each pick
%corresponding to that channel centerline. We use the first pick as the
%reference pick.

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
            end
            cpd(pick,m,n) = mean(sqrt(cl2));
        end
    end
    [null_dump cp(:,pick)] = min(squeeze(cpd(pick,:,:))');
end
    clear lov lovpp

%********** Using "Channel Pair" Information, Scrub 'Candidate' (Low) Quality Data Points ************************    
for m = 1:max(d(1).cid);
    IsBad(1,m) = (nanmean(d(1).TrQuality(d(1).cid==m)))>=8;%Quality = 8 is 'Candidate' (Low) Quality; Lower Values are Higher Quality
end

for pk = 1:3
    BadPoint(pk).b = 0.*[1:length(d(pk).cid)];
end

for bad = find(IsBad(1,:)==1)
    BadPoint(1).b(d(1).cid==bad) = 1;
    BadPoint(2).b(d(2).cid==cp(bad,2)) = 1;
    BadPoint(3).b(d(3).cid==cp(bad,3)) = 1;
end
    
%Clear 'Candidate' quality data
for pk = 1:3
    d(pk).cid(BadPoint(pk).b==1) = [];
     d(pk).vx(BadPoint(pk).b==1) = [];    
     d(pk).vy(BadPoint(pk).b==1) = [];    
     d(pk).vz(BadPoint(pk).b==1) = [];    
     d(pk).hiz(BadPoint(pk).b==1) = []; 
     d(pk).ctxz(BadPoint(pk).b==1) = [];    
     d(pk).cz(BadPoint(pk).b==1,:) = [];    
end

%********* Pepper Points Uniformly Along Channel Centerlines & Find Inflection Points ********************************************************

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)
        if sum(d(pk).cid == m) > 0;
        %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)]);
        if (WhichTransect == 1 || WhichTransect == 2) %Interpolated stratigraphic surfaces exist
            for iss = 1:5 %iss = interpolated stratigraphic surface
                d(pk).p(m).cz(:,iss) = ...
                interp1(d(pk).lac(d(pk).cid==m),d(pk).cz(d(pk).cid==m,iss),[0:dbip:d(pk).ctl(m)]);
            end
        end
        
        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
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 Using Inflection Points, Setting Aside Low-Sinuosity Wiggles **********************************************

for pk = 1:3 %for each pick
    for m = 1:max(d(pk).cid) %for each centerline trace
        if sum(d(pk).cid == m) > 0;
        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
           if (WhichTransect == 1 || WhichTransect == 2) %Stratigraphic surfaces are defined
           for iss = 1:5
                d(pk).h(m).cz(epi,iss,1:2) = ...
                    [d(pk).p(m).cz(d(pk).p(m).hmeplist(epi),iss) d(pk).p(m).cz(d(pk).p(m).hmeplist(epi+1),iss)];%z,start|end
           end
           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. 
%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;

        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
                if (WhichTransect == 1 || WhichTransect == 2) %Stratigraphic surfaces are defined
                for iss = 1:5
                d(pk).h(m).cz(epi,iss,1:2) = ...
                    [d(pk).p(m).cz(d(pk).p(m).hmeplist(epi),iss) d(pk).p(m).cz(d(pk).p(m).hmeplist(epi+1),iss)];%z,start|end
                end
                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
end

%****************** Compute Summary Parameters & Plot Results ***************************************************************

    mi = zeros(size(d)); mreji = zeros(size(d));
for m = 1:max(d(pk).cid) %channel-centerline number
 
        if picturesque>=2;figure;end
    for pk = 1:3 %each channel-centerline is picked 3 times
    if sum(d(pk).cid == cp(m,pk)) > 0;

        %PLOTTING
        if picturesque>=1
            if isfield(d(pk).p(m),'z') == 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')
            else
                 scatter(d(pk).p(cp(m,pk)).x,d(pk).p(cp(m,pk)).y,10,'Filled')
            end
         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
         %END OF PLOTTING       
                
        
        %FIND SUMMARY PARAMETERS FOR PICK-BY-PICK HALF-MEANDER
        if  d(pk).h(cp(m,pk)).halfmeanderspresent ~= -1;
        for n = 1:size(d(pk).h(cp(m,pk)).sinu,2) %For each half-meander
            if d(pk).h(cp(m,pk)).sinu(n) > minsinu %provided that sinuosity is OK
                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',MeanderColor{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,:);
             if (WhichTransect == 1 || WhichTransect == 2); %A stratigraphic elevation is defined;
                 for iss = 1:5 %interpolated stratigraphic surface
                    mcz(pk,mi(pk),iss,:) = d(pk).h(cp(m,pk)).cz(n,iss,:);
                 end
             end
                 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 OF FINDING SUMMARY PARAMETERS FOR PICK-BY-PICK HALF-MEANDER

    end
    end
    
    
    
    if picturesque>=2
    for pk = 1:3 %plot text labels on top
    if sum(d(pk).cid == cp(m,pk)) > 0;

        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
    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 (WhichTransect == 1 || WhichTransect == 2); %A stratigraphic elevation is defined;
    meancz = mean(mcz,4);
end


%************************* Identify Replicable Half-Meanders By Pairing-Up Single Picks *************************
% A real half-meander should be identifiable in repeated picks of the same
% channel. Therefore, we compare the three picks to seek half-meanders that
% are reproduced on multiple 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 (possible) 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',MeanderColor{rem(ghm,7)+1})
    %                                  hold on
    %                                  line([mean([x1 x2]) mean(x3)],[mean([y1 y2]) mean(y3)],'LineWidth',2,'Color',MeanderColor{rem(ghm,7)+1})
    %                                  text(mean([x1 x2]),mean([y1 y2]), ...
    %                                  num2str(dcls_fr(ghm,wi)),'FontSize',25,'Color',MeanderColor{rem(ghm,7)+1})
    %                                  
    %                                  
    %                                  title(strcat('ghm:',num2str(ghm),',wi:',num2str(wi)))
    %                                  scatter(hmslx,hmsly,30,MeanderColor{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);

    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));
        meanz(length(meanz):length(meanx)) = NaN;
    if (WhichTransect == 1 || WhichTransect == 2);
        for iss = 1:5
            meanczthisiss = squeeze(meancz(:,:,iss));
            rhm(l).cz(iss) = nanmean(  meanczthisiss(mrhmi==l));
        end
    end
    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 Centerline Trace ******************************

%Combine the mean and standard deviation of replicable half-meanders on the same centerline traces to
%obtain a characteristic half-meander wavelength and error for
%each picked channel.

for m = 1:length(d(1).p)%For each channel ID
    rhmtc = find([rhm(:).m] == m);%replicable half-meanders, this channel
        
    if isempty(rhmtc)==0

    chm(m).lambda = sum([rhm(rhmtc).quality].*[rhm(rhmtc).lambda]) ...
            ./sum([rhm(rhmtc).quality]);
    chm(m).x = sum([rhm(rhmtc).quality].*[rhm(rhmtc).x]) ...
            ./sum([rhm(rhmtc).quality]);
    chm(m).y = sum([rhm(rhmtc).quality].*[rhm(rhmtc).y]) ...
            ./sum([rhm(rhmtc).quality]);   
    chm(m).z = sum([rhm(rhmtc).quality].*[rhm(rhmtc).z]) ...
            ./sum([rhm(rhmtc).quality]);    
        
    if (WhichTransect == 1 || WhichTransect == 2);
    for iss = 1:5
       chm(m).cz(iss) = sum([rhm(rhmtc).quality].*[rhm(rhmtc).cz(iss)]) ...
                        ./sum([rhm(rhmtc).quality]);  
    end
    end    
    %Based on Wikipedia article
    %'Standard_deviation#Combining_standard_deviations' as of 24FEB2014
    %Non-overlapping populations
    chm(m).err = sqrt( ( sum([rhm(rhmtc).quality].*( [rhm(rhmtc).err].^2 + [rhm(rhmtc).lambda].^2 )) ...
                    / sum([rhm(rhmtc).quality]))...
                    - chm(m).lambda.^2);
    else
        chm(m).lambda   = NaN;
        chm(m).err = NaN;
        chm(m).x = NaN;
        chm(m).y = NaN;
        chm(m).z = NaN;
    end
end

%Debug plot
figure
    errorbar([1:length(d(1).p)],[chm.lambda],[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  ******************************
%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).

%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.

%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:max(size(cp))
    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
        if picturesque == 1
        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
        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],2.*[rhm(:).lambda],[rhm(:).err],'xk')
    hold on
    scatter([rhm(:).z],2.*[rhm(:).lambda],30,[rhm(:).x],'Filled')
    box on;grid on; colorbar
    set(gca,'FontSize',18);xlabel('Raw elevation');ylabel('Meander wavelength')
    title('Meander wavelength versus raw elevation (individual half-meanders)')

    figure
    errorbar([chm(:).z],2.*[chm(:).lambda],[chm(:).err],'xk')
    hold on
    scatter([chm(:).z],2.*[chm(:).lambda],30,[chm(:).x],'Filled')
    box on;grid on; colorbar
    set(gca,'FontSize',18);xlabel('Raw elevation');ylabel('Meander wavelength')
    title('Meander wavelength versus raw elevation (collating along picks)')
    
    
    if (WhichTransect == 1 || WhichTransect == 2) %Elevation corrections are available
        figure
        for irhm = 1:length(rhm)
            czadjusted(irhm) = [rhm(irhm).z] - [rhm(irhm).cz(2)];
        end
        errorbar(czadjusted,2.*[rhm(:).lambda],[rhm(:).err],'xk')
        hold on
        scatter(czadjusted,2.*[rhm(:).lambda],30,[rhm(:).x],'Filled')
        box on;grid on; colorbar
        title('Meander wavelength versus corrected elevation')
        set(gca,'FontSize',18);xlabel('Corrected elevation (individual half-meanders)');ylabel('Meander wavelength')
        
        for irhm = 1:length(rhm)
            czadjusted(irhm) = [chm(irhm).z] - [chm(irhm).cz(2)];
        end
        errorbar(czadjusted,2.*[chm(:).lambda],[chm(:).err],'xk')
        hold on
        scatter(czadjusted,2.*[chm(:).lambda],30,[chm(:).x],'Filled')
        box on;grid on; colorbar
        title('Meander wavelength versus corrected elevation (collating along picks)')
        set(gca,'FontSize',18);xlabel('Corrected elevation');ylabel('Meander wavelength')

        
        
        figure
        czadjusted = [rhm(:).z] - (cetl + [rhm(:).x].*pfxc(WhichTransect) + [rhm(:).y].*pfyc(WhichTransect));
        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)')
    end
    
    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'
        meanderoutputname = strcat('parseMeanderWavelengthMeasurementsMarsRivers_output_terse',num2str(WhichTransect))
        save(meanderoutputname,'rhm','d','chm')
    case 'complete'
        disp('Not saving complete output.')
    case 'none'
end