%% ---------------------------------------------------------------------------
%                    Master file for meandering river workflow
% ----------------------------------------------------------------------------
%
% atom_data_array column order
% 1 - cutoff time
% 2 - preserved size at the end of reworking time
% 3 - original size
% 4 - preservation percent
% 5 - median age of the bar
% 6 - centroid X
% 7 - centroid Y
% 8 continuos size 
% 9 - distance from channel belt axis
% 10 - Major Axis length
% 11 - Minor Axis length
% 12 - St.Dev of ages


%% ---------------------------------------------------------------------
%                    Controls
% ------------------------------------------------------------------
to_run_model = 1;
edit_centerlines = 0;
to_make_grid = 1;
gen_reworking = 1;
gen_ages = 0;
gen_channel_belt_axis = 0;
gen_maps_n_discharge=0;
gen_atoms = 0;
gen_atom_data = 0;
gen_bar_preservation_plot = 1;
to_gen_bar_stats = 0;
gen_data = 0;
gen_age_distribution = 0;
calc_sinuosity = 0;

% Misc variables
years = 10000; % How many simulated years would you like
dt_yrs = 1;
sim_time = ceil(years/dt_yrs);
statstoolbox = 1; % enter 1 if you have the Statistics Toolbox, else enter 0
age_mapping_threshold = 0.005;

grid_width = 700;
grid_height = 50;
%% ------------------------------------------------------------------------
%                     Run the Model
%------------------------------------------------------------------------
channel_belt_enforce = 1;
if to_run_model == 1
    [river, B, E, nodecount] = migration_model(statstoolbox,years,dt_yrs, channel_belt_enforce); 
end

%% --------------------------------------------------------------------------
%                    Centerline editing
% -------------------------------------------------------------------------
if edit_centerlines == 1
    river_truncated(numel(river)).Xcl = [];
    river_truncated(numel(river)).Ycl = [];
    nodecount_truncated(numel(river)).A_cutoff_rem = [];
    nodecount_truncated(numel(river)).B_spacing_ins =[];
    nodecount_truncated(numel(river)).C_duplicate_rem = [];
    for ii = 1:numel(river)-1
        [r,c] = find(river(ii).Ycl < 4200 & river(ii).Ycl > -4200);
        river_truncated(ii).Xcl(r,1) = river(ii).Xcl(r);
        river_truncated(ii).Ycl(r,1) = river(ii).Ycl(r);
        if  (~isempty(river(ii).Ycl(nodecount(ii).A_cutoff_rem))) && max((max(river(ii).Ycl(nodecount(ii).A_cutoff_rem))) < -4200);
                nodecount_truncated(ii).A_cutoff_rem = [];
                nodecount_truncated(ii).B_spacing_ins = nodecount(ii).B_spacing_ins;
                nodecount_truncated(ii).C_duplicate_rem = nodecount(ii).C_duplicate_rem;
        else
                nodecount_truncated(ii).A_cutoff_rem = nodecount.A_cutoff_rem;
                nodecount_truncated(ii).B_spacing_ins = nodecount(ii).B_spacing_ins;
                nodecount_truncated(ii).C_duplicate_rem = nodecount(ii).C_duplicate_rem;
        end
        if  (~isempty(river(ii).Ycl(nodecount(ii).A_cutoff_rem))) && max((max(river(ii).Ycl(nodecount(ii).A_cutoff_rem))) > 6000);
                nodecount_truncated(ii).A_cutoff_rem = [];
                nodecount_truncated(ii).B_spacing_ins = nodecount(ii).B_spacing_ins;
                nodecount_truncated(ii).C_duplicate_rem = nodecount(ii).C_duplicate_rem;
        else
                nodecount_truncated(ii).A_cutoff_rem = nodecount.A_cutoff_rem;
        end
    end
%{
    figure(5)
    subplot(2,1,1)
    imagesc(meandering_grid(100).map);
    subplot(2,1,2)
    plot(river_truncated(100).Xcl,river_truncated(100).Ycl);
    %}
end


%% --------------------------------------------------------------------------
%                    Generate the grid 
% -------------------------------------------------------------------------
% from the river structure centerlines
% Note: I changed river structure to river_truncated because I shorted the
% existing structure
% Depending on if I edit the centerline - I have to change the name of the
% river structure in the make_grid function

if to_make_grid == 1
    [meandering_grid, X_edges, Y_edges] = make_grid(river, grid_width, grid_height);
end


%% -----------------------------------------------------------------------
%                   Fluvial surface reworking
% -------------------------------------------------------------------------
% Note: 
%
starting_iteration = 1000;
if gen_reworking ==1;
    [reworking, Imagesum] = meandering_fluvial_surface_reworking(meandering_grid, grid_height, grid_width,starting_iteration, dt_yrs);

    figure(1)
    imagesc(Imagesum);
    figure(2)
    x = reworking(1:((numel(meandering_grid)-starting_iteration)/2),1);
    y = reworking(1:((numel(meandering_grid)-starting_iteration)/2),2);
    scatter(x,y)
end

%% -----------------------------------------------------------------------
%                   Generate the age map
% ----------------------------------------------------------------------
% Note: 
reworking_time = 29000;
snapshot_iteration = reworking_time + starting_iteration;
snapshot_length = length(starting_iteration:reworking_time);
if gen_ages == 1
    [ages] = age_map_meandering(meandering_grid,grid_height, grid_width,starting_iteration,reworking_time,dt_yrs);
end  

%% ------------------------------------------------------------------------
%                    Generate Channel-belt axis
% ------------------------------------------------------------------------
if gen_channel_belt_axis == 1
    CB_position=(sum(cat(1,river.Ycl)))/(length(cat(1,river.Ycl)));
    CB_grid_height = discretize(CB_position,Y_edges);
    
    ages_keep = ages(reworking_time-1).map;
    CB_grid = zeros(grid_height,grid_width);
    for ii = 1:grid_width
        [r,c] = find(ages_keep(:,ii) >= 0);
        CB(ii) = (sum(r)/length(r));
    end
    CB_smooth = round(smooth(CB,451));
    for ii = 1:grid_width
        CB_grid(CB_smooth(ii),ii) = 1;
    end
end
%% -----------------------------------------------------------------------
%           Color coded age map and critical discharge threshold
% -------------------------------------------------------------------------
%
if gen_maps_n_discharge == 1;
    
    [nr, nc] = size(ages_keep); %ages(reworking_time).map
    cmap = parula(max(ages_keep(:)));
    % Make values 0-5 black:
    yellow = [1 1 0.4];
    cmap(1:10,:) = zeros(10,3);
    colormap(cmap);
    
    figure(3)
    %imagesc(ages(reworking_time).map);
    pcolor([ages_keep nan(nr,1); nan(1,nc+1)]);
    colormap(cmap);
    colorbar
    %caxis([10 14000])
    shading flat;
    set(gca, 'ydir', 'reverse');
    hold on
    axis off;
 
    % Threshold discharge map
    % Note: 
    threshold_discharge = meandering_grid(20000).map;
    %threshold_discharge(threshold_discharge < pctile) = 0;
    %threshold_discharge(threshold_discharge > pctile) = 1;

    figure(4)
    imagesc(threshold_discharge);
    
end    

%% -----------------------------------------------------------------------
%             Atom identification and preservation
% ------------------------------------------------------------------------
if gen_atoms == 1
    ages_keep = ages(reworking_time-1).map;
    %clear ages
    % Run the atom_ages function
    % from this we generate the atom tracks, the gridded atoms, and the
    % atom_data strcuture which has sizes 
    [atom, atom_data,atom_data_array] = atom_ages(river,nodecount, B,grid_height, grid_width,starting_iteration,reworking_time, X_edges, Y_edges)
end

%% -----------------------------------------------------------------------
%            Atom data: atom age, atom preservation, atom size
% -----------------------------------------------------------------------

if gen_atom_data == 1
    % Finds the equivalent grid cell ages for the bar packages and the
    % median age of the bar packages
    %atom_data(numel(atom_data)).preserved_ages = [];   % Commented out
    %because I changed the structure into array
    %atom_data(numel(atom_data)).median_ages = [];
    for ii = 1:numel(atom_data)
        atom_data(ii).preserved_ages = NaN(grid_height, grid_width);
        [r, c] = find(atom_data(ii).preserved_grid == 1);
        row_length = length(r);
        for pp = 1:row_length
            atom_data(ii).preserved_ages(r(pp),c(pp)) = (ages_keep(r(pp),c(pp))./reworking_time-1);
        end
        atom_data_array(ii,5) = median(median(atom_data(ii).preserved_ages, 'omitnan'),'omitnan');
    end

    % STATS for the continuous bars
    % Decide whether I should just move this down to the other data
    % collection seciton... probably will
    for ii = 1:numel(atom_data)
        bar_map_binary = im2bw(atom_data(ii).preserved_grid,0.99);
        CC = bwconncomp(bar_map_binary,4);
        s = regionprops(CC,'centroid','MajorAxisLength','MinorAxisLength','Area');
        centroids = cat(1,s.Centroid);
        MajorAxes = cat(1,s.MajorAxisLength);
        MinorAxes = cat(1,s.MinorAxisLength);
        Area = cat(1,s.Area);
        % Create a single array with all of the output data
        % Columns:
        % 1- bar # from the atom finding
        % 2- centroid X
        % 3- centroid y
        % 4- major axis length
        % 5- minor axis length
        % 6- pixel area
        if isempty(centroids)
            atom_bar_data_full(ii,:) = [0 0 0 0 0 0];
        else
            components_found = length(centroids(:,1));
            bar_number = ones(components_found,1) * ii;
            atom_bar_data = [bar_number centroids Area MajorAxes MinorAxes];
            [M,I] = max(atom_bar_data(:,4));
            [r,c] = ind2sub(size(atom_bar_data),I);
            atom_bar_data_full(ii,:) = atom_bar_data(r,:);
        end
    end
    
    %atom_data(numel(atom_data)).centroidX = [];
    %atom_data(numel(atom_data)).centroidY = [];
    %atom_data(numel(atom_data)).continuous = [];
    
    % Creating cutoff times and st.dev of bar ages
    %atom_data_array = NaN(numel(atom_data),12;
    for ii = 1:numel(atom_data)
        atom_age = cat(1,atom_data(ii).preserved_ages(:));
        age_stdv = std(atom_age,'omitnan');
        atom_data_array(ii,12) = age_stdv;
        atom_data_array(ii,1) = atom(ii).t(1);
    end
   
    atom_data_array(:,6) = atom_bar_data_full(:,2);
    atom_data_array(:,7) = atom_bar_data_full(:,3);
    atom_data_array(:,8) = atom_bar_data_full(:,4);
    atom_data_array(:,10) = atom_bar_data_full(:,5);
    atom_data_array(:,11) = atom_bar_data_full(:,6);
    
    % Calculating Standard deviation of the ages
    for ii = 1:numel(atom_data)
        if atom_data_array(ii,2) == 0
            %atom_data_array(ii,1) = NaN;
            atom_data(ii).preserved_grid = NaN;
             atom_data(ii).preserved_ages = NaN;
            atom_data_array(ii,6) = NaN;
            atom_data_array(ii,7) = NaN;
            atom_data_array(ii,8) = NaN;
            atom_data_array(ii,10) = NaN;
            atom_data_array(ii,11) = NaN;
        end
        if isnan(atom_data_array(ii,7))
            atom_data_array(ii,9) == NaN;
        else
            atom_data_array(ii,9)= abs(atom_data_array(ii,7) - CB_smooth(round(atom_data_array(ii,7))));
        end
    end
    % clear some unused variables 
end  
%clear atom_bar_data atom_bar_data_full CC centroids
%% -----------------------------------------------------------------------
%           Making composite maps of bar preservation and bar age
% --------------------------------------------------------------------------
if gen_bar_preservation_plot == 1   
    % For plotting. I'm generating a map of the specific atoms/bars with
    % their associated ages
    for ii = 1:numel(atom_data)
        if ii == 1
            bar_composite = atom_data(ii).grid;
        else
            % Add the next bar on top of the other, but set the overlapping
            % portion to the ID from the most recent bar
            bar_composite = bar_composite + atom_data(ii).grid;
            [r,c] = find(bar_composite > ii);
            for pp = 1:length(r)
                bar_composite(r(pp),c(pp)) = ii;
            end
        end
    end    
    % Uses the bar ID to reset the ID to the percentage of original atom
    % size preserved
    preservation_composite = bar_composite;
    for ii = 1:numel(atom_data)
        [r,c] = find(preservation_composite==ii);
        if isempty(r)
            ;
        else
            for pp = 1:length(r)
                preservation_composite(r(pp),c(pp)) = atom_data_array(ii,2);
            end
        end
    end
    % Uses the bar ID to reset the ID to the median ages calculated above.
    % Each bar is now coded by median bar age
    age_composite = bar_composite;
    for ii = 1:numel(atom_data)
        [r,c] = find(age_composite==ii);
        if isempty(r)
            ;
        else
            for pp = 1:length(r)
                age_composite(r(pp),c(pp)) = atom_data_array(ii,5);
            end
        end
    end
    
   preservation_composite = preservation_composite(:,100:600);
   age_composite = age_composite(:,100:600);
    
    centroidX = atom_data_array(:,6);
    centroidY = atom_data_array(:,7);    
    
    figure(5)
    subplot(2,1,1)
    imagesc(preservation_composite);
    colorbar
    subplot(2,1,2)
    imagesc(age_composite);
    colorbar
    hold on
    %scatter(centroidX,centroidY,'r*');
    %scatter([1:1:700],CB_smooth,'r*')
    
    colorbar
    
end


%% -----------------------------------------------------------------------
%             Creating the threshold maps and binary file
% ------------------------------------------------------------------------
if gen_data == 1
    threshold_map = ages_keep./(reworking_time); %ages(reworking_time).map
    threshold_map(threshold_map < age_mapping_threshold) = NaN;
    threshold_map(threshold_map >= age_mapping_threshold) = 1;


%% Set up the grids
% Call in a fresh version of the equivalent age_map over an equivalent grid
% area
    ages_map = ages_keep; %ages(reworking_time).map;
    threshold_map_binary = im2bw(threshold_map,0.99); % Convert the threshold map from model output to a binary image
    threshold_map_show = zeros(grid_height,grid_width); % set up fresh map for showing the isolated threshold maps

%% ---------------------------------------------------------------------
%         Labeling the largest bars
% ------------------------------------------------------------------------
% Find connected regions
label_map = bwlabel(threshold_map_binary,4); % creates labeled grid of connected grid components
label_maximum = max(label_map); % finds how many labels the function found
maximum= max(label_maximum);

if maximum == 0
    disp('there are no continuous bars');
else
    label(maximum).list = []; % fresh list to find the largest bar feature

%% Sorting Labeled bars
% Finds the rows and columns associated with each bar feature and stores
% its length to a growing array of index lengths
    for ii = 1:maximum
        [r, c] = find(label_map==ii);
        label(ii).list = [r c];
        label_length(ii) = length(label(ii).list); 
    end
    [sorted_labels,label_idx] = sort(label_length, 'descend'); % indexes the longest label array (therefore the largest bar)

%% Plotting the biggest (most continuous bars) onto grids w/ their age values
    ages_map_show = zeros(grid_height,grid_width); % fresh grid to show the isolated age maps
    
    store = 1;
for pp = 1:length(label_idx)
    row_col = label(label_idx(pp)).list; % cleaning up the path for calling in the indexed label positions
    c = makeColorMap([0.75 0.75 0.74],[0.2 0.2 1], [0 0.6 0.2], 20);
    
% Overprints zero-array ages_map_show with values at the position of the
% biggest bar
    for ii = 1:length(row_col(:,1))
        ages_map_show(row_col(ii,1),row_col(ii,2)) = ages_map(row_col(ii,1),row_col(ii,2));
    end
end
    
end
end

%% ------------------------------------------------------------------------
%           Creating histogram for the age distribution on the bar
% -------------------------------------------------------------------------   
 if gen_age_distribution == 1;
     
     ages_map = ages(reworking_time).map./reworking_time;
     age_distribution(numel(label)).list = [];
    for ii = 1:numel(label)
        run_length = length(label(ii).list(:,1));
        age_distribution(ii).list = zeros(run_length,1);
        for pp = 1:run_length
            age_distribution(ii).list(pp,1) = ages_map(label(ii).list(pp,1),label(ii).list(pp,2));
        end
     
% Deciding which size to plot with individual histogram    
    ages_map_show = zeros(grid_height,grid_width); % fresh grid to show the isolated age maps
   
    for pp = 1:1 % if all length(label_idx), otherwise pick
        row_col = label(label_idx(pp)).list; % cleaning up the path for calling in the indexed label positions
        c = makeColorMap([0.75 0.75 0.74],[0.2 0.2 1], [0 0.6 0.2], 20);
    
% Overprints zero-array ages_map_show with values at the position of the
% biggest bar
        for ii = 1:length(row_col(:,1))
            ages_map_show(row_col(ii,1),row_col(ii,2)) = ages_map(row_col(ii,1),row_col(ii,2));
        end
    end
    % Plotting
    figure(6)
    subplot(2,1,1)
    imagesc(ages_map_show);
    subplot(2,1,2);
    histogram(age_distribution(1).list,30)
    
       
      
    end
 end
 
 %% -----------------------------------------------------------------------
 %             Calculate the sinuosities
 % ------------------------------------------------------------------------
 % Note: I shorted the river structure and renamed some variables
 % 'river_truncated' I will have to change this back fro the the normal
 % output
if calc_sinuosity == 1
    for pp = 1:numel(river)
        river_length = length(river(pp).Xcl);
        X = river(pp).Xcl;
        Y = river(pp).Ycl;
        Xstart = river(pp).Xcl(1,1);
        Ystart = river(pp).Ycl(1,1);
        Xend = river(pp).Xcl(river_length,1);
        Yend = river(pp).Xcl(river_length,1);
        Xdiff = Xend - Xstart;
        Ydiff = Yend - Ystart;
        Lv(pp) = sqrt((Xdiff.^2)+ (Ydiff.^2));
        dS = sqrt(diff(X).^2+diff(Y).^2);
        Lc(pp) = sum(dS);
        Sinuosity(pp) = Lc(pp) / Lv(pp);
    end   
    length_Sin = length(Sinuosity)
    figure(9)
    W = [1:length_Sin];
    scatter(W, Sinuosity, 'filled');
end