%Edwin Kite

%******************* Set Control Parameters **************************
dbip = 10; %distance between interpolated points (m)
smoothn = 10; %number of interpolated points to smooth curvature for (to prevent artifactual inflection picks)
minsinu = 1.1; %minimum half-meander sinuosity to be considered as a half-meander

%******************* Set Graphics & Housekeeping Parameters **************************
EdwinColor = {'k','b','g','y','r',[.5 .6 .7],[.8 .2 .6]};

%******************* Import Data **************************
%1'XCoord' 2'YCoord' 3'FID' 4'ID' 5'PQUALITY' 6'PSTYLE' 7'TRAMBIG'
%8'TRQUALITY' 9'SUREMEANDR' 10'ORIG_FID' 11'POINT_X' 12'POINT_Y' 13'038882HIZ'
%14'548104HIZ' 15'B20G02CTXZ'
a =  dlmread('Transect_1_Channel_Centerlines_for_Third_Pass_Meander_Wavelength_Measurements_csv',',',1,0);
%1 'XCoord' 2'YCoord' 3'FID' 4'ID' 5'ORIG_FID' 6'POINT_X' 7'POINT_Y'
%8 '038882HIZ' 9'548104HIZ' 10'B20G02CTXZ
a2 = dlmread('Transect_1_Pick_2_Third_Pass_Centerlines_for_Meander_Wavelength_Measurements_Points_csv',',',1,0);
a3 = dlmread('Transect_1_Pick_3_Third_Pass_Centerlines_for_Meander_Wavelength_Measurements_csv',',',1,0);

%******************* Preprocess Data **************************

d(1).cid  = a(:,10)+1;           %ArcGIS is zero-indexed for feature IDs
d(1).vx   = a(:,1) - min(a(:,1)); %vertex x
d(1).vy   = a(:,2) - min(a(:,2)); %vertex y
assert(max(max(a(:,13:15)))<=0,'Positive elevation values will cause errors')
d(1).hiz  = min(a(:,13),a(:,14));
d(1).ctxz = a(:,15);
%HiRISE elevations are preferred. CTX gapfilling is needed rarely. CTX frequently is '0' when HiRISE has a
%value.
d(1).ctxz(d(1).ctxz==0) = NaN;
d(1).hiz(d(1).hiz==0)   = NaN;
d(1).ctxz_adjusted = d(1).ctxz - nanmean(d(1).ctxz - d(1).hiz); %Apply offset to bring CTXz mean into line with HiZ
d(1).vz   = nanmin(d(1).ctxz_adjusted,d(1).hiz);

for pick = 2:3
    switch pick
        case 2
            csvdata = a2;
        case 3
            csvdata = a3;
    end
    d(pick).cid  = csvdata(:,5)+1; %ArcGIS is zero-indexed for feature IDs
    d(pick).vx   = csvdata(:,1) - min(csvdata(:,1));
    d(pick).vy   = csvdata(:,2) - min(csvdata(:,2));
    d(pick).hiz  = min(csvdata(:,8),csvdata(:,9));
    d(pick).ctxz = csvdata(:,10);
    d(pick).ctxz(d(pick).ctxz==0) = NaN;
    d(pick).hiz(d(pick).hiz==0)   = NaN;
    d(pick).ctxz_adjusted = d(pick).ctxz - nanmean(d(pick).ctxz - d(pick).hiz); %Apply offset to bring CTXz mean into line with HiZ
    d(pick).vz   = nanmin(d(pick).ctxz_adjusted,d(pick).hiz);
end


%********* Pepper Points Uniformly Along Channel Centerlines **************

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

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

        d(pick).p(m).th1plusth2 = [NaN  d(pick).p(m).smooth(1:end-1) + d(pick).p(m).smooth(2:end)    ];
        d(pick).p(m).th2plusth3 = [     d(pick).p(m).smooth(1:end-1) + d(pick).p(m).smooth(2:end) NaN];

        %Following Hemberger, 1986 Master's thesis, U. Virginia: 
        d(pick).p(m).isinfl = [1 - (sign(d(pick).p(m).th1plusth2) == sign(d(pick).p(m).th2plusth3))]; 
    end
end

hold on

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

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

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

%******************* Define Half-Meanders **************************
for pick = 1:3
    for m = 1:max(d(pick).cid)
        d(pick).p(m).okhmep = d(pick).p(m).isinfl; %O.K. half-meander end points
        d(pick).p(m).okhmep(1) = 0;        %no inflections at end of channel picks
        d(pick).p(m).okhmep(end) = 0;      %no inflections at end of channel picks
        d(pick).hmeplist = find(d(pick).p(m).okhmep==1);
        if length(d(pick).hmeplist)>1
        for epi = 1:(length(d(pick).hmeplist)-1) %end-point index 
            ep = d(pick).hmeplist(epi);
            d(pick).h(m).p(epi,1,1:2) = ...
                [d(pick).p(m).x(d(pick).hmeplist(epi)) d(pick).p(m).x(d(pick).hmeplist(epi+1))];%x,start|end
            d(pick).h(m).p(epi,2,1:2) = ...
                [d(pick).p(m).y(d(pick).hmeplist(epi)) d(pick).p(m).y(d(pick).hmeplist(epi+1))];%y,start|end
            d(pick).h(m).p(epi,3,1:2) = ...
                [d(pick).p(m).z(d(pick).hmeplist(epi)) d(pick).p(m).z(d(pick).hmeplist(epi+1))];%z,start|end
            d(pick).h(m).ls(epi) = ... %half-meander straight-line distance
                sqrt( ...
                (d(pick).h(m).p(epi,2,2) - d(pick).h(m).p(epi,2,1))^2 ...
               +(d(pick).h(m).p(epi,1,2) - d(pick).h(m).p(epi,1,1))^2);
            d(pick).h(m).lc(epi) = (d(pick).hmeplist(epi+1) - d(pick).hmeplist(epi))*dbip;%along-curve distance
            d(pick).h(m).th(epi) = nanmean(d(pick).p(m).th(d(pick).hmeplist(epi):d(pick).hmeplist(epi+1))); %for debugging
            d(pick).h(m).sinu(epi) = d(pick).h(m).lc(epi) ./ d(pick).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(pick).oksinu = d(pick).h(m).sinu>minsinu;

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

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

            for epi = 1:(length(d(pick).hmeplist)-1) %end-point index 
                ep = d(pick).hmeplist(epi);
                d(pick).h(m).p(epi,1,1:2) = ...
                    [d(pick).p(m).x(d(pick).hmeplist(epi)) d(pick).p(m).x(d(pick).hmeplist(epi+1))];%x,start|end
                d(pick).h(m).p(epi,2,1:2) = ...
                    [d(pick).p(m).y(d(pick).hmeplist(epi)) d(pick).p(m).y(d(pick).hmeplist(epi+1))];%y,start|end
                d(pick).h(m).p(epi,3,1:2) = ...
                    [d(pick).p(m).z(d(pick).hmeplist(epi)) d(pick).p(m).z(d(pick).hmeplist(epi+1))];%z,start|end
                d(pick).h(m).ls(epi) = ... %half-meander straight-line distance
                    sqrt( ...
                    (d(pick).h(m).p(epi,2,2) - d(pick).h(m).p(epi,2,1))^2 ...
                   +(d(pick).h(m).p(epi,1,2) - d(pick).h(m).p(epi,1,1))^2);
                d(pick).h(m).lc(epi) = (d(pick).hmeplist(epi+1) - d(pick).hmeplist(epi))*dbip;%along-curve distance
                d(pick).h(m).sinu(epi) = d(pick).h(m).lc(epi) ./ d(pick).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

        end
    end
end

%****** Compute Summary Parameters & Plot Results for Debugging *********
figure
    mi = zeros(size(d)); mreji = zeros(size(d));
for m = 1:max(d(pick).cid) %channel-centerline number
    figure
    for pick = 1:3 %each channel-centerline is picked 3 times

                 scatter(d(pick).p(m).x,d(pick).p(m).y,20,d(pick).p(m).z,'Filled')
         hold on        
                 scatter(d(pick).p(m).x(find(d(pick).p(m).isinfl)),d(pick).p(m).y(find(d(pick).p(m).isinfl)),200,'k*')

    %            scatter(p(m).x,p(m).y,100,p(m).isinfl,'Filled')
                hold on
        for n = 1:size(d(pick).h(m).sinu,2)
            if d(pick).h(m).sinu(n) > minsinu
                line([d(pick).h(m).p(n,1,1) d(pick).h(m).p(n,1,2)], ...
                     [d(pick).h(m).p(n,2,1) d(pick).h(m).p(n,2,2)],'Color',EdwinColor{pick})
                text(d(pick).h(m).p(n,1,1)+50,d(pick).h(m).p(n,2,1),num2str(round(d(pick).h(m).ls(n))))
                mi(pick) = mi(pick) + 1;
                mls(pick,mi(pick)) = d(pick).h(m).ls(n);
                 mx(pick,mi(pick)) = mean(d(pick).h(m).p(n,1,:));
                 my(pick,mi(pick)) = mean(d(pick).h(m).p(n,2,:));
                 mz(pick,mi(pick)) = mean(d(pick).h(m).p(n,3,:));
                 mm(pick,mi(pick)) = m;
              msinu(pick,mi(pick)) = d(pick).h(m).sinu(n);
            else
                line([d(pick).h(m).p(n,1,1) d(pick).h(m).p(n,1,2)], ...
                     [d(pick).h(m).p(n,2,1) d(pick).h(m).p(n,2,2)],'Color','r')    

                mreji(pick) = mreji(pick) + 1;
                mrejsinu(pick,mreji(pick)) = d(pick).h(m).sinu(n); 
            end
        end
    
    end
            title(num2str(m))
            axis equal
end
figure
scatter(mz,mls,30,mx,'Filled');grid on
title('Half-meander straight-line distance versus elevation')

figure
hist([mrejsinu msinu],[1:0.01:4])
title('Histogram of all sinuosities')
