
clear;
clc;
close all;
C1=0;
rez=[0;0];
global cmap;
global Plot;


Plot=true;
for res=25%[25,50,100, 200, 400]%, 800, 1000];
  probx=0.26%[0.35:0.01:0.4]%0.38%0.42%0.38%[0.35:0.01:0.41]%[0.21:0.1:0.51] %0.55;

  for i=1:length(probx)
perc_prob=1-probx(i)
prob=probx(i)
for j=1%1:10%2%1:10

    %     if(Plot) figure %('position',[-1200,200,800,800]);
%     end
        
% randSeed=6;
% rand('state', randSeed);
%s tem so narejene originalne slike, ki dajo napake
% randSeed=9;
% rand('state', randSeed);
randSeed=j;
rand('state', randSeed);

T1=rand(res,res);
% % 
II=find(T1>prob);  
% % 
T1(:,:)=zeros(res,res);
T1(II)=1;
% 
%0 is white, 1 is black
%set nord boundary and south boundary to 1
T1(1,:)=1;
T1(size(T1,1),:)=1;
%set left boundary and right boundary to 1
T1(:,1)=0;
T1(:,size(T1,2))=0;

% %%if the system is obtained from a figure
% T1=[];
% % T11=imread('Kicma4Bojenje5 view.png');
% %T11=imread('Two_ cluters.png');
% T11=imread('Objasnjenje.png');
% 
% T11=T11(:,:,1);
% T11=T11/max(max(T11));
% for i=1 : size(T11,1) 
%     for j=1:size(T11,2)
%     if T11(i,j)==0 T11(i,j)=1; else T11(i,j)=0; end
%     end
% end
% 
% k=1;
% % j=6 %6 pixels per sqare in 'Kicma4Bojenje5 view.png'
% % j=20 %20 pixels per sqare in 'Two_ cluters.png'
% j=20 %20 pixels per sqare in 'Objasnjenje.png'
% 
% for i=1 :j: size(T11,1)  %20 pixels per sqare
%     l=1;
%     for m=1:j:size(T11,2)
%         T1(k,l)=T11(i,m);
%         l=l+1;
%     end
%     k=k+1;
% end
% 
% %ADES paper patch for black cells
% T1=double(T1);
% T1(3,5)=1; 
% T1(4:5,2)=1;
% T1(5,9)=1; 
% T1(6,14:16)=1;T1(6,19)=1;
% T1(7,18)=1; T1(7,20)=1;
% T1(12,14)=1;
% T1(14,15:20)=1;T1(14,8)=1;
% T1(15,9)=1;
% T1(16,2:3)=1;T1(16,10)=1;T1(16,19)=1;
% T1(18,20)=1;
% %%end from figure

%     %Figure 1: My case with error - modified by Biljana (one column inserted on right)
%     T1=[0 1 1 1 1 1 0;   0 0 0 1 0 1 0;  0 1 1 0 1 1 0;   0 1 1 1 0 1 0; ...
%         0 0 1 1 1 1 0;   0 0 1 0 1 0 0;  0 1 1 1 1 1 0 ];
% %     ADES Figure 1
%      T1=[0 1 1 1 1 1 1 1 1 1 0;  0 1 1 0 0 0 0 1 1 0 0;  0 0 1 0 0 0 0 0 1 1 0; ...
%          0 0 1 0 1 0 1 1 1 0 0;  0 0 0 1 1 1 1 0 1 1 0;  0 0 0 0 0 0 1 1 1 1 0; ...
%          0 1 1 1 0 1 1 0 1 1 0;  0 1 0 1 0 1 0 0 0 0 0;  0 1 1 1 0 1 0 1 1 1 0; ...
%          0 0 1 1 0 1 0 1 0 1 0;  0 1 1 1 1 1 1 1 1 1 0 ];
%     %A bar
%     T1=[0 1 1 1 1 1 0;   0 0 0 1 0 0 0;  0 0 0 1 0 0 0;   0 0 0 1 1 0 0; ...
%         0 0 0 0 1 0 0;   0 0 0 0 1 0 0;  0 1 1 1 1 1 0 ];
% %     Loop
%      T1=[0 1 1 1 1 1 1 1 1 1 0;  0 0 1 0 0 0 0 0 0 0 0;  0 0 1 1 1 1 1 0 1 1 0; ...
%          0 0 0 0 1 0 1 0 1 0 0;  0 0 0 1 1 1 1 0 1 1 0;  0 0 0 0 0 0 1 1 1 1 0; ...
%          0 1 1 1 0 1 1 0 1 1 0;  0 1 0 1 0 1 0 0 0 0 0;  0 1 1 1 0 1 0 1 1 1 0; ...
%          0 0 1 1 0 1 0 1 0 1 0;  0 1 1 1 1 1 1 1 1 1 0 ];
%      T1=flipud(T1);
% %     Biljana case: simple case
%      T1=[0 1 1 1 1 1 1 1 1 1 0;  0 0 1 0 0 0 0 0 0 0 0;  0 0 1 1 1 1 1 0 1 1 0; ...
%          0 0 0 0 1 1 0 0 1 0 0;  0 0 0 1 1 1 1 0 1 1 0;  0 0 0 1 0 0 1 1 1 1 0; ...
%          0 1 1 1 1 1 1 0 1 1 0;  0 1 0 1 0 1 1 0 0 0 0;  0 1 1 1 0 1 1 0 0 0 0; ...
%          0 0 0 0 0 1 0 0 0 0 0;  0 1 1 1 1 1 1 1 1 1 0 ];
% %     Biljana case: simple case 1
%      T1=[0 1 1 1 1 1  1 1 1 0;  0 0 1 0 0 0  0 0 0 0;  0 0 1 1 1 1  1 1 0 0; ...
%          0 0 0 0 1 1  0 1 0 0;  0 0 0 1 1 1  0 1 1 0;  0 0 0 1 0 1  1 0 1 0; ...
%          0 1 1 1 1 1  0 1 1 0;  0 1 0 1 0 1  0 0 0 0;  0 1 1 1 0 1  0 0 0 0; ...
%          0 1 1 1 1 1  1 1 1 0 ];
% %     %Biljana case: donji èošak
%      T1=[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0; ...
%          0 0 0 0 0 1 1 0 1 0 0 0 1 0 0 0 1 0 1 0 0; ...
%          0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0; ...
%          0 0 0 0 0 1 0 0 1 1 1 1 1 0 1 1 1 1 1 0 0; ...
%          0 0 0 0 0 1 0 0 1 1 1 0 1 0 1 0 1 0 0 0 0; ...
%          0 0 0 0 0 1 1 1 1 1 1 0 1 0 1 1 1 0 0 0 0; ...
%          0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0 0 0 0 0; ...
%          0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0; ...
%          0 0 0 0 0 0 0 0 1 1 0 0 1 1 1 1 0 0 0 0 0; ...
%          0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0; ...
%          0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0; ...
%          0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0; ...
%          0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0; ...
%          0 1 1 1 1 1 1 1 0 0 1 1 0 1 1 1 1 0 0 0 0; ...
%          0 0 0 1 0 0 0 1 1 1 1 0 1 1 1 0 0 0 0 0 0; ...
%          0 0 1 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0; ...
%          0 0 1 1 0 0 1 0 1 1 1 0 1 1 0 0 0 0 0 0 0; ...
%          0 0 0 1 1 1 1 0 1 1 1 1 1 1 1 0 0 0 0 0 0; ...
%          0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0; ...
%          0 0 1 1 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0; ...
%          0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0];

% % %     %Figure 2 left - Biljana case: in paper a)
% T1=[0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	0; ...
% 0	1	1	1	0	1	0	0	1	1	0	1	1	0	0	1	0	0	1	0	0; ...
% 0	1	0	1	0	0	0	1	1	0	0	0	1	0	1	0	1	0	0	1	0; ...
% 0	0	1	1	1	1	1	0	1	1	1	1	1	0	1	1	1	0	1	0	0; ...
% 0	0	1	1	1	1	1	1	1	1	0	0	1	1	1	1	0	0	0	0	0; ...
% 0	1	1	0	1	0	0	0	0	0	0	1	0	1	0	0	0	0	1	1	0; ...
% 0	1	0	1	0	1	1	1	1	0	1	1	1	1	1	0	1	1	1	0	0; ...
% 0	0	1	1	1	1	1	0	1	0	1	1	1	0	0	0	1	1	0	1	0; ...
% 0	0	1	0	0	0	0	1	0	1	0	1	1	1	1	1	0	1	0	1	0; ...
% 0	1	0	1	1	1	1	0	0	0	1	1	1	0	0	1	1	1	1	0	0; ...
% 0	0	1	1	0	1	1	0	1	0	1	0	0	0	1	1	0	0	0	1	0; ...
% 0	0	0	0	1	1	0	0	1	1	0	0	1	0	0	0	0	0	1	1	0; ...
% 0	1	0	1	0	1	0	1	1	0	1	1	0	0	1	1	1	1	1	1	0; ...
% 0	1	1	1	1	1	0	0	0	1	1	0	1	0	1	0	1	0	0	1	0; ...
% 0	1	1	1	1	0	0	1	0	1	1	1	1	1	1	0	0	1	0	0	0; ...
% 0	1	1	1	0	1	0	1	1	1	1	1	1	0	1	0	1	1	1	1	0; ...
% 0	0	1	1	1	0	1	1	1	1	1	1	1	0	1	0	1	1	1	1	0; ...
% 0	1	1	1	1	1	1	0	0	1	0	1	1	1	1	0	1	1	0	0	0; ...
% 0	1	1	0	0	0	0	1	0	0	1	1	1	1	0	1	0	0	1	1	0; ...
% 0	1	1	1	0	0	1	0	1	1	0	1	0	1	0	1	1	1	0	1	0; ...
% 0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	0]; 


% % %     %Figure 2 right -Biljana case: in paper b)
%      T1=[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0; ...
%          0 0 0 1 1 1 1 0 1 1 0 0 0 0 1 0 1 0 1 1 0; ...
%          0 0 1 0 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 0; ...
%          0 0 0 1 0 0 1 0 0 0 0 0 0 1 1 1 0 0 1 1 0; ...
%          0 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 1 1 0; ...
%          0 1 0 1 1 1 1 0 0 1 1 1 0 0 1 1 0 0 1 0 0; ...
%          0 1 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 0 1 0; ...
%          0 0 1 1 0 1 1 1 0 0 1 1 1 0 1 1 1 1 0 1 0; ...
%          0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 0 0 0 0; ...
%          0 0 1 0 0 0 1 0 1 1 0 0 0 1 1 1 0 0 1 0 0; ...
%          0 0 0 0 1 1 1 1 1 0 1 1 1 0 0 1 0 1 1 1 0; ...
%          0 1 0 0 0 1 1 1 1 0 1 1 0 0 1 1 1 0 0 1 0; ...
%          0 1 1 1 0 0 0 1 0 0 0 1 0 0 0 1 1 1 0 1 0; ...
%          0 0 1 1 1 1 1 1 1 0 0 0 1 0 1 0 1 0 1 1 0; ...
%          0 1 0 0 0 0 0 0 1 0 1 1 1 0 0 0 1 0 1 1 0; ...
%          0 1 1 1 0 1 0 1 1 1 0 1 0 0 1 1 1 1 1 1 0; ...
%          0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 1 0 1 1 0 0; ...
%          0 0 0 1 1 0 0 1 0 1 1 1 0 1 1 0 1 0 0 0 0; ...
%          0 1 0 1 0 0 1 0 0 1 1 0 1 1 1 1 1 1 1 1 0; ...
%          0 1 1 1 1 1 1 0 1 0 0 0 0 0 1 0 1 0 0 1 0; ...
%          0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0];
% 


%%%%  Envelope of zeros is added because of %nnnn$ neighborhood
T1= [zeros(size(T1,1),1),T1,zeros(size(T1,1),1)];
T1= [zeros(1,size(T1,2));T1;zeros(1,size(T1,2))];


   Num_allc = ((size(T1,1)-2)*(size(T1,2)-2))
   Num_nzc = nnz(T1)  %number of black cells
    
        %Colors used in automata
        %     colormap('gray');
    %     colormap('default');
    %    cmap=colormap; 
    cmap_dim = (((size(T1,1)-1)*(size(T1,2))-1) + 12 + 4); % X times Y unique colors 
    %5 special colors white, black, green, alpha, beta and 12 tagging  colors
    cmap= colormap(hot(cmap_dim)); 
%     cmap= colormap(cool(cmap_dim)); 
%     cmap= colormap(jet(cmap_dim)); 
%     cmap= colormap(gray(cmap_dim)); 

%     cmap = zeros(cmap_dim,3);
%     cmap(:,1) = 0:1/(cmap_dim-1):1;
%     cmap(:,2) = 1:-1/(cmap_dim-1):0;
%     cmap(:,3) = 0:1/(cmap_dim-1):1;
    
    cmap(1,:) = [0.8 0.8 1];  %white (the smalest value equal to 0) becomes purple 
    cmap(length(cmap),:) = [0 1 0]; %green is the largest valuer reserved for the backbone color equal to length(cmap)
    cmap(length(cmap)-1,:) = [0 0 0]; %black is next smaller value length(cmap)-1
    cmap(length(cmap)-2,:) = [0 0 1]; %alpha (blue) is next smaller value length(cmap)-2
    cmap(length(cmap)-3,:) = [1 0 0]; %beta (red) is next smaller value length(cmap)-2
    
    cmap(length(cmap)-4,:) = [0 0 0.35];
    cmap(length(cmap)-5,:) = [0 0 0.4];
    cmap(length(cmap)-6,:) = [0 0 0.45];
    cmap(length(cmap)-7,:) = [0 0 0.5];
    cmap(length(cmap)-8,:) = [0 0 0.55];
    cmap(length(cmap)-9,:) = [0 0 0.6];
    cmap(length(cmap)-10,:) = [0 0 0.65];
    cmap(length(cmap)-11,:) = [0 0 0.7];
    cmap(length(cmap)-12,:) = [0 0 0.75];
    cmap(length(cmap)-13,:) = [0 0 0.8];
    cmap(length(cmap)-14,:) = [0 0 0.85];
    cmap(length(cmap)-15,:) = [0 0 0.9];
    
    white = 1;%length(cmap)-2;
    green = length(cmap);
    black = length(cmap)-1;
    alpha = length(cmap)-2;
    beta = length(cmap)-3;
   
    BW=T1;
    BWW=BW+1;
    if(Plot) figure(1); colormap(cmap); image(BWW((2:size(BWW,1)-1),(2:size(BWW,2)-1)));  set(gcf,'renderer','OpenGL');end %show 0 as black and 1 as white
%     if(Plot) figure(1); colormap(cmap); image(BW+1); title(sprintf('Percolation site p= %f', prob));  set(gcf,'renderer','OpenGL'); end %show 0 as black and 1 as white
    CC = bwconncomp(BW,4);
    numPixels = cellfun(@numel,CC.PixelIdxList);
    [biggest,idx] = max(numPixels);
    
    BW(CC.PixelIdxList{idx}) = 3; %maximal cluster is now black
    if(Plot) figure; colormap(cmap); image(BW((2:size(BW,1)-1),(2:size(BW,2)-1)));% image(BW),  
        %title('Infinite cluster'); 
    end
        
    T_C=T1*black; %convert the values 1 to black color (length(cmap)-1) and 0 to white color (length(cmap)-2)
    for lines=1:size(T_C,1)
        for columns=1:size(T_C,2)
            if (T_C(lines, columns)== 0) T_C(lines, columns)= white; end
        end
    end
    if(Plot)  figure; colormap(cmap); image(T_C((2:size(T_C,1)-1),(2:size(T_C,2)-1))); title('original'); colorbar; set(gcf,'renderer','OpenGL'); end
    
%CA_A1 - Color uniqly white cells (that are surrounded on e,se and s with black cells)
    tic;[T2,N] = CA_A1(T_C, white, black); C1=[C1,toc];
    N_CA1=N/Num_allc
    N_total =N_CA1;   
%     if(Plot) figure; colormap(cmap); image(T2((2:size(T2,1)-1),(2:size(T2,2)-1))); title('CA\_A1'); colorbar; set(gcf,'renderer','OpenGL'); end
    if(Plot) figure; colormap(cmap); image(T2((2:size(T2,1)-1),(2:size(T2,2)-1)));  set(gcf,'renderer','OpenGL'); end
    %title('CA\_A1');
    
%CA A2 boji bijele nnn komponente razlicitim bojama tako sto za svaku ne crnu celiju bira najvecu ne bijelu i ne crnu boju iz nnn susjedstva
    tic;[T3,N] = CA_A2(T2,white, black); C1=[C1,toc];
    N_CA2=N/Num_allc
    N_total = N_total + N_CA2;
    if(Plot)  figure; colormap(cmap); image(T3((2:size(T3,1)-1),(2:size(T3,2)-1)));  set(gcf,'renderer','OpenGL'); end
    if (T3(size(T3,1)-1,1)==T3(size(T3,1)-1,size(T3,2)-1)) display('No infinite cluster'); end
    %title('CA\_A2'); 
    
%CA A3 boji artikulacione cvorove
    tic;[T4,N] = CA_A3(T3,white, black); C1=[C1,toc];
    N_CA3=N/Num_allc
    N_total = N_total +N_CA3;
    if(Plot)  figure; colormap(cmap); image(T4((2:size(T4,1)-1),(2:size(T4,2)-1)));   set(gcf,'renderer','OpenGL'); end
    %title('CA\_A3');
    
%CA A4 boji dio backbone-a u rezerviranu boju (zelenu)
    tic;[T5,N] = CA_A4(T4,white, black); C1=[C1,toc];
    N_CA4=N/((size(T5,1)-3)*(size(T5,1)-2))
    N_total = N_total + N_CA4;
    if(Plot)  figure; colormap(cmap); image(T5((2:size(T5,1)-1),(2:size(T5,2)-1)));  set(gcf,'renderer','OpenGL'); end
    %title('CA\_A4'); 
    
%CA_A5 - Color uniqly white cells (that are surrounded on e,se and s with black cells)
    tic;[T6,N] = CA_A5(T5,white, black); C1=[C1,toc];
    N_CA5=N/Num_allc
    N_total = N_total + N_CA5;
    if(Plot)  figure; colormap(cmap); image(T6((2:size(T6,1)-1),(2:size(T6,2)-1)));  set(gcf,'renderer','OpenGL'); end
    %title('CA\_A5');
    
%%ADES
%CA_A6 - Color green ceells by \alpha
    tic;[T7,N] = CA_A6(T6,white, black); C1=[C1,toc];
    N_CA6=N/((size(T7,1)-4)*(size(T7,1)-2))
    N_total = N_total + N_CA6;
    if(Plot)  figure; colormap(cmap); image(T7((2:size(T7,1)-1),(2:size(T7,2)-1)));  set(gcf,'renderer','OpenGL'); end
    %title('CA\_A6');
    
%CA_A7 - Color real BB cells with \beta
    tic;[T8,N] = CA_A7(T7,white, black); C1=[C1,toc];
    N_CA7=N/((size(T8,1)-4)*(size(T8,1)-2))
    N_total = N_total + N_CA7
    if(Plot)  figure; colormap(cmap); image(T8((2:size(T8,1)-1),(2:size(T8,2)-1)));  set(gcf,'renderer','OpenGL'); end
    %title('CA\_A7');
%%end ADES

 %Clear the remaining cells and make an overlay on the infinite cluster
    T9=T8;
    for lines=1:size(T9,1)
        for columns=1:size(T9,2)
            if (T8(lines, columns)==beta) 
                BW(lines, columns)= beta;
            else
                T9(lines, columns)= white; 
            end
        end
    end
    
    if(Plot)  figure; colormap(cmap);  image(T9((2:size(T9,1)-1),(2:size(T9,2)-1)));%image(T9); 
        %title('BB');  
        set(gcf,'renderer','OpenGL'); end
    if(Plot) figure; colormap(cmap); image((BW((2:size(BW,1)-1),(2:size(BW,2)-1))));   set(gcf,'renderer','OpenGL'); end
    %title('Infinite cluster (black) overlayed by Backbone (green)');

    Elapsed_time_step(j) = sum(C1); C1=[];
    N_ts(j)=N_total;
    Num_black(j)=Num_nzc;
    IC(j)=biggest;
    BB(j)=length(find(T8==beta));
    Ratio_IC(j)= IC(j)/BB(j);  %Ratio of IC cells : BB cells
    Ratio_total(j)= Num_black(j)/BB(j);   %Ratio of total nuber of cells : BB
%     if (T3(end,1)==T3(end,end)) N_ts(j)=NaN; end
end
    Elapsed_time(i) = sum(Elapsed_time_step);
    N_ts_m(i)=round(mean(N_ts));
    Num_black_m(i)=round(mean(Num_black));
    IC_m(i)=round(mean(IC));
    BB_m(i)=round(mean(BB));
    Ratio_IC_m(i)= IC_m(i)/BB_m(i);  %Ratio of IC cells : BB cells
    Ratio_total_m(i)= Num_black_m(i)/BB_m(i);   %Ratio of total nuber of cells : BB
end
res
Elapsed_time
perc_prob_m=1-probx    
N_ts_m
Num_black_m
IC_m
BB_m
Ratio_IC_m
Ratio_total_m
end

