clear; close all; more off;

dr = '~/Documents/Large Melt Pond Images/GFL Images/'; % change this to your directory

files = dir([dr '*.tif']);

nf = length(files);

ps = 100:50:1500; %subimage dimensions to loop over

%setup arrays for results

infRes = NaN(length(ps),1);

for i=1:length(ps)
    
    afRes{i} = []; % fraction of image that is water

    percRes{i} = []; % 1 if percolation

end

for k = 1:nf
    
    disp(k);
    
    %download image
    imageData = imread(files(k).name);

    BW = im2bw(imageData, .7); % water = 0, ice = 1
    %change the second input to adjust how much is converted to 0/1
    %around .5 is recommended to disclude gray areas, while .7 includes these

    BW = imcomplement(BW); % water = 1, ice = 0

    %figure; imshow(BW); %can be used to check bw conversion

    [sz1,sz2] = size(BW);

    %for each image, slice into multiple smaller images
    
    for j=1:length(ps)
        
        %finds maximum number of n by n squares fully fitting into image
        
        width = (floor(sz1/ps(j))-1);
        height = (floor(sz2/ps(j))-1);
        
        perc = NaN(width*height,1);
        af = perc;
        
        for h=1:height
            
            for w=1:width
                
                %loop for subimage calculation
                
                pn = w + (h-1)*width; %number assignment for each subimage (left to right and down rows)
                
                img = BW((w-1)*ps(j)+1:w*ps(j),(h-1)*ps(j)+1:h*ps(j)); %extract subimage

                af(pn) = sum(img(:)==1)/length(img(:)); % fraction of image that is water

                lk = bwlabel(img,4); %labels groups of nearest neighbors (water)

                lk(lk==0) = NaN; %converts 0s to NaN so that the ismember function won't read it
                
                %checks for percolation

                ptb = find(ismember(lk(1,:),lk(end,:)));

                pss = find(ismember(lk(:,1),lk(:,end)));

                if ~isempty(ptb) || ~isempty(pss)

                    perc(pn) = 1;
                    %returns 1 if percolation has occurred

                else

                    perc(pn) = 0;

                end
                
            end

        end
        
        % store image info into arrays categorized by subimage pixel width 
        %this decreases runtime because images are only loaded once
        
        afRes{j} = cat(1,afRes{j},af);
        percRes{j} = cat(1,percRes{j},perc);

    end

end

for k=1:length(ps)
    
    %figure( 'Name', strcat('Area Fraction vs. p',int2str(ps(k)))); plot(afRes{k},percRes{k},'bo');
    
    %binning results

    binsize = .025;
    bins = 0:binsize:1;

    p = discretize(afRes{k},bins); %organizes results into bins

    %counts the number of instances in each bin
    y = accumarray(p(:),1);
    z = zeros(length(y),1);

    %convert bin res back to area fraction
    %counts number of perc=1 in each bin
    for i = 1:length(p)
        z(p(i)) = z(p(i)) + percRes{k}(i);
        p(i) = p(i)*binsize - binsize*.5; %+ .25); %%%remember--min value must be added
    end

    %averages instances of percolation per bin
    for j = 1:length(z)
        z(j) = z(j)/y(j);
    end

    %remove NaNs from z (empty bins)
    z = z(~isnan(z));

    %remove duplicate values to get one point for each filled bin
    p = unique(p);

    %figure('Name', strcat('Area Fraction vs. p',int2str(ps(k)))); plot(p,z,'bo');
    
    %% Fit: 'Area Fraction vs. p'.
    [xData, yData] = prepareCurveData(p,z);

    % Set up fittype and options.
    ft = fittype( 'a * atan((x*(pi/2) - b)*c)+d', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    
    %determine bounds for fit
    if ps(k) >= 1000
        opts.Lower = [0.33 -Inf 30 0.5]; %provides a better fit when subimages get larger and
        opts.Upper = [0.35 Inf Inf 0.5]; %there are fewer data points to fit the curve to
    end
    opts.StartPoint = [0.33 0.5 100 0.5];

    % Fit model to data.
    [fitresult, gof] = fit( xData, yData, ft, opts );

    % Plot fit with data.
    figure( 'Name', strcat('Percolation Threshold of Real Ponds','_',int2str(ps(k))));
    h = plot( fitresult, xData, yData );
    legend(h, 'Binned Data', 'Arctan Fit','location','southeast');
    % Label axes
    xlabel 'Area Fraction'
    ylabel 'Probability of Percolation'
    grid on
    
    %find inflection point/calculate percolation threshold
    cf = coeffvalues(fitresult);
    syms x
    f = cf(1) * atan((x*(pi/2) - cf(2))*cf(3)) + cf(4);
    f1 = diff(f);
    f2 = diff(f1);
    infRes(k) = fzero(matlabFunction(f2),0.34); %store result of inflection point in array
    
    %export figures
%     folder = '~/Percolation Threshold Graphs';
%     pngFileName = sprintf(strcat('PercolationThresholdRealPonds_',num2str(ps(k)),'.png'));
% 	fullFileName = fullfile(folder, pngFileName);
% 	export_fig(fullFileName);
    
end

figure( 'Name', strcat('Sub Image Size vs. Percolation Threshold'));
h = plot(ps,infRes,'bo');
xlabel('Image Width in Pixels')
ylabel('Percolation Threshold')
grid on

%export_fig Image_Width_vs_Percolation_Threshold_Real_Ponds.png