rng(0);
cluster = 0;

% code to create a random matrix representing the percolation sites
nx = 512;
ny = 512;
r1 = rand(nx*ny,1);
M = reshape(r1, [nx,ny]);

p = 0.62;

latice = zeros(nx,ny);
latice(M < p) = 1;

% make a copy of our matrix
latice_b = latice*(-1);

% arrays
i_idx = zeros(nx*ny,1);
j_idx = zeros(nx*ny,1);

i = 1;
j = 1;
count = 1;
step = 1;

for j = 1:ny
    if (latice_b(i,j) == -1)
        i_idx(count) = i;
        j_idx(count) = j;
        
        latice_b(i,j) = step;
        count = count +1;
    end
end

n_max = 0;
n_min = -1;

%
% Latice "burning"
%

% while there are unexplored sites to check for connections
while (n_max+1 ~= count)
    step = step + 1;
    n_min = n_max+1;
    n_max = count-1;
    %fprintf('n: %d %d\n', n_min, n_max);
    for (k = n_min:n_max)
        i = i_idx(k);
        j = j_idx(k);
        %fprintf('  %d, %d, %d\n', k, i, j);
        % check to see if the point below has a latice site
        if (i < nx && latice_b(i+1,j) == -1)
            i_idx(count) = i+1;
            j_idx(count) = j;
            latice_b(i+1,j) = step;
            count = count + 1;
        end
        
        % check to see if the point to the right has a latice site
        if (j < ny && latice_b(i,j+1) == -1)
            i_idx(count) = i;
            j_idx(count) = j+1;
            latice_b(i,j+1) = step;
            count = count + 1;
        end
        
        % check to see if the point to the left has a latice site
        if (j > 1 && latice_b(i,j-1) == -1)
            i_idx(count) = i;
            j_idx(count) = j-1;
            latice_b(i,j-1) = step;
            count = count + 1;
        end
    end
end

last_row = latice_b(end,:);
row_idx = last_row > 0;
path = 0;

if (sum(row_idx) > 0)
    path = 1;
    fprintf('path found\n');
else
    fprintf('no path found\n');
end

%
% Path finding
%

if (path)
    latice_p = latice_b;
    last_row = latice_p(end,:);
    
    row_idx = last_row == 0;
    last_row(row_idx) = -(nx*ny+1);
    
    i = nx;
    [max_num, j]=max(last_row(:));
    step = last_row(j);
    length = step;
    
    if (step == 0 || step == -(nx*ny+1))
        fprintf('no path found, error');
    end
    
    latice_p(i,j) = -latice_p(i,j);
    
    while (i > 1)
        if (latice_p(i-1,j) == step -1)
            i = i-1;
        elseif (j > 1 && latice_p(i,j-1) == step -1)
            j = j-1;
        elseif (j < ny && latice_p(i,j+1) == step -1)
            j = j+1;
        elseif (i < nx && latice_p(i+1,j) == step -1)
            i = i+1;
        end
        step = step - 1;
        latice_p(i,j) = -latice_p(i,j);
    end
    
    %
    % Create and image
    %
    
    I = mat2gray(latice_p, [min(min(latice_p)) max(max(latice_p))]);
    figure;
    imshow(I);
end

%
% Create clusters
%

if (cluster)
    
    labels = 1:nx*ny;
    
    count = 0;
    label_latice = zeros(nx,ny);
    for i = 1:nx
        for j = 1:ny
            if (latice(i,j) == 1)
                if (i > 1)
                    above = label_latice(i-1,j);
                else
                    above = 0;
                end
                
                if (j > 1)
                    left = label_latice(i,j-1);
                else
                    left = 0;
                end
                %fprintf('i: %d j: %d\n', i, j);
                
                % new label
                if (left == 0 && above == 0)
                    %fprintf('above: %d left: %d\n', above, left);
                    count = count + 1;
                    label_latice(i,j) = count;
                    %fprintf('i: %d j: %d, a\n', i, j);
                elseif (left ~= 0 && above == 0)
                    %fprintf('above: %d left: %d\n', above, left);
                    [label_latice(i,j), labels] = find_val(left,labels);
                    %fprintf('i: %d j: %d, b\n', i, j);
                elseif (left == 0 && above ~= 0)
                    %fprintf('above: %d left: %d\n', above, left);
                    [label_latice(i,j), labels] = find_val(above,labels);
                    %fprintf('i: %d j: %d, c\n', i, j);
                else
                    labels = union(left,above,labels);
                    [label_latice(i,j),labels] = find_val(left,labels);
                    %fprintf('d');
                end
                
                labels;
            end
        end
    end
    
    for i = 1:nx
        for j = 1:ny
            if (label_latice(i,j) ~= 0)
                while (label_latice(i,j) ~= labels(label_latice(i,j)))
                    label_latice(i,j) = labels(label_latice(i,j));
                end
            end
        end
    end
    
    I = mat2gray(label_latice, [min(min(label_latice)) max(max(label_latice))]);
    figure;
    imshow(I);
end