% correlation length g is proportional to exp(-r/xi),
% where r is distance and xi is critical exponent
% xi is proportional to |p-p_c|^-nu
% so log(xi) = -nu*log(|p-p_c|) + xi_0

clear; close all; more off;

% make list of files %Melt Pond Images

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

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

nf = length(files);

ps = 1301; %pixel dimension of subimages

cent = (ps-1)/2+1;

%setup arrays for results

%area fraction bins
afBinsize = .025;
afBins = 0:afBinsize:1;

for j=1:length(afBins)
    for i=1:(cent-1)
        disRes{j,i} = []; % fraction of image that is water
    end
end

% disRes is a cell array of points binned by the image's area fraction
% and their distance from center point. The value inputted is a 0 or 1
% depending on whether or not it belonged to the center cluster

imgNum = 0;

for k = 1:nf
    
    disp(strcat('imgNum',int2str(k))); %shows which image is being processed
    
    %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
        
    width = (floor(sz1/ps)-1);
    height = (floor(sz2/ps)-1);
    
    for h=1:height

        for w=1:width

            pn = w + (h-1)*width;

            img = BW((w-1)*ps+1:w*ps,(h-1)*ps+1:h*ps);
            
            if img(cent,cent) == 1
                
                %fraction of area that is water
                af = sum(img(:)==1)/length(img(:));
                
                %bin the area fraction
                for j=1:length(afBins)
                    if af < afBins(j+1)
                        af = j;
                        break
                    end
                end
                
                imgNum = imgNum + 1;
                
                disp(imgNum) %tracks progress--how many subimages have been checked         

                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
                
                ls = size(lk);
                
                %calulate distances for each element in array
                
                for i=1:(ls(1)*ls(2))
                    
                    %disp(i) %check program speed
                     
                    %calculates distance from center point
                    dis = sqrt(((i - ls(1)*floor(i/ls(1))-cent)^2)+((floor(i/ls(1))+1-cent)^2));
                    
                    % checks if point is within circle of data and not the
                    % center point
                    if dis <= (cent-1) && dis > 0
                        
                        dis = int16(dis);
                        
                        % bins results--1 if belonging to same cluster
                        % as center point, 0 if not
                        if lk(i - ls(1)*floor(i/ls(1)),floor(i/ls(1))+1) == lk(cent,cent)
                            
                            disRes{af,dis} = cat(1,disRes{af,dis},1);
                            
                        else
                            
                            disRes{af,dis} = cat(1,disRes{af,dis},0);
                            
                        end
                    
                    end
                    
                end

            end

        end

    end

end

%convert results to proportions and graph
for j=1:length(afBins)
    
    %temporary arrays for graphing
    disArray = [];
    propArray = [];
    
    %loop to add results to arrays
    for i=1:(cent-1)
        if ~isempty(disRes{j,i})
            
            summed = sum(disRes{j,i});
            total = length(disRes{j,i});
            
            propArray = cat(1,propArray,summed/total);
            disArray = cat(1,disArray,i);
        end
    end
    
    if ~isempty(disArray)
    
        %graph results by area fraction
        figure( 'Name', strcat('Distance vs. Proportion in Cluster',num2str(afBins(j)))); plot(propArray,disArray,'bo');
        
    end
    
end