
function [af,l1] = dispersion()
clear
clc
close all
 
% V5 %
% October 20, 2010 %


fprintf('DISPERSION CHARACTERIZATION CODE \nMaterials Tribology Laboratory\nUniversity of Delaware\n\n')

fprintf('Please refer to accompanying document for further instructions \n\tand illustrations on using this code.\n\n')

syms logi1 y n;

% Version check for SORT function
logi1 = input('Is the version of MATLAB 7.0 or later? (Y/N) : ', 's');

if strcmpi(logi1, 'y') ==1;
    vzn = 1;
else
    vzn = 0;
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Read TEM b/w bitmap and generate pixel matrix 'A' %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[filename,PathName] = uigetfile('*.bmp','Select the .BMP image');
 
A = imread([PathName,filename],'bmp');
 
offx=ceil(size(A,1)/2); %Offset for periodic boundary
offy=ceil(size(A,2)/2); %Offset for periodic boundary
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Area fraction calculation %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filledcells=0;
totalcells=0;
 

for i=1:size(A,1);
    for j=1:size(A,2);
        totalcells=totalcells+1;
        if A(i,j)<150;
            filledcells = filledcells+1;
        end
    end
end
af=filledcells/totalcells;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Create data matrix which contains periodic boundary condition %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%% Top layer: Top quarter of the data matrix
for i = 1:offx;
    for j = 1:offy;
        Matrix(i,j)=A(offx+i,(offy+j));
    end
end
 
for i = 1:size(A,1);
    for j = 1:offy;
        Matrix(i+offx,j)=A(i,j+offy);
    end
end
 
for i = 1:offx;
    for j = 1:offy;
        Matrix(i+offx+size(A,1),j)=A(i,j+offy);
    end
end
 
%%%%%%%%% Middle layer: Middle two-rows of data matrix
for i = 1:offx;
    for j = 1:size(A,2);
        Matrix(i,j+offy)=A(i+offx,j);
    end
end
 
for i = 1:size(A,1);
    for j = 1:size(A,2);
        Matrix(i+offx,j+offy)=A(i,j);
    end
end
 
for i = 1:offx;
    for j = 1:size(A,2);
        Matrix(i+offx+size(A,1),j+offy)=A(i,j);
    end
end
 
%%%%%%%%%%%%% Bottom layer: Bottom quarter of data matrix
for i = 1:offx;
    for j = 1:offy;
        Matrix(i,j+offy+size(A,2))=A(i+offx,j);
    end
end
 
for i = 1:size(A,1);
    for j = 1:offy;
        Matrix(i+offx,j+offy+size(A,2))=A(i,j);
    end
end
 
for i = 1:offx;
    for j = 1:offy;
        Matrix(i+offx+size(A,1),j+offy+size(A,2))=A(i,j);
    end
end
 
syms logi y n;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Define initial square parameters %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


scale1= input('Enter approximate width of TEM image in nanometers:  ');
scale= scale1/size(A,2); 
xf1 = input('Enter initial guess of characteristic square width in nanometers:  ');
xf = xf1/scale1;
Nsquares = input('Enter number of random squares to use for analysis:  ');

 
fprintf('\n\n')
fprintf('Computation may take several seconds to a few minutes, depending on input parameters\n Press Ctrl+C any time to terminate computation\n\n')


s2=size(A,1)*size(A,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% First count to Nsquares (corresponding to manual input of 'xf') %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


length = ceil(xf*size(A,2));
for i = 1:Nsquares;
    x=ceil(rand*size(A,1))+offx;
    y=ceil(rand*size(A,2))+offy;
    counter=0;
    for s = 1:length;
        for t = 1:length;
            mx=x-floor(length/2)+s;
            my=y-floor(length/2)+t;
            if Matrix(mx,my)<150;
                counter=counter+1;
            end
        end
    end
    P(i)=counter;
end
data=P';
bimreal=mode(data);

h1=hist(P,2500);
[m,n1]=max(h1);
l=xf;


boxar=length^2;

if boxar==bimreal
    isbimo=1;
    h1(1,2500)=0;
    [m,n1]=max(h1);
else
    isbimo=0;
end



    
 
 if n1~=1
     
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% Automatic interation to obtain Free Space Length %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
margin=1;
while margin > .05
    
    while n1~=1
        xfnew=xf/2;
        l=xf-xfnew; %==xf/2
        length = ceil(l*size(A,1));
        for i = 1:Nsquares;
          x=ceil(rand*size(A,1))+offx;
          y=ceil(rand*size(A,2))+offy;
          counter=0;
          for s = 1:length;
              for t = 1:length;
                mx=x-floor(length/2)+s;
                my=y-floor(length/2)+t;
                if Matrix(mx,my)<150;
                  counter=counter+1;
                end
              end
          end
        P(i)=counter;
        end
    data=P';
    bimreal=mode(data);
    boxar=length^2;
    
    h2=hist(P,2500);
    
    if boxar==bimreal
        isbimo=1;
        h2(1,2500)=0;
        [m,n1]=max(h2);
    else
        isbimo=0;
        [m,n1]=max(h2);
    end

             
    xf=xfnew;
           
    end
    
    if n1==1
        if vzn == 1;
            SRT=sort(h2,'descend');
        else
            SRT=fliplr(sort((h2)));
        end
        
        margin=1-(SRT(2)/SRT(1));
        
        if margin > 0.05
            xfjumpup=2.5*xf;
            xf=xfjumpup;
           
            xfnew=xf/2;
             l=xf-xfnew; %==xf/2
             length = ceil(l*size(A,1));
             for i = 1:Nsquares;
                x=ceil(rand*size(A,1))+offx;
                y=ceil(rand*size(A,2))+offy;
               counter=0;
                for s = 1:length;
                     for t = 1:length;
                     mx=x-floor(length/2)+s;
                      my=y-floor(length/2)+t;
                     if Matrix(mx,my)<150;
                         counter=counter+1;
                     end
                     end
                end
                    P(i)=counter;
               end
         data=P';
         bimreal=mode(data);
         boxar=length^2;
    
    
         h2=hist(P,2500);
         
         if boxar==bimreal
            isbimo=1;
            h2(1,2500)=0;
            [m,n1]=max(h2);
         else
            isbimo=0;
            [m,n1]=max(h2);
         end
    

         
         xf=xfnew;    
    
        else
        end
    else
        continue
        
 
    end

    
end
end
 
    
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Display results %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



    
    l1=l*size(A,2)*scale;
    disp('Value of required box width is (in nanometers): ')
    disp(ceil(l1))
    

 
disp('Particle area fraction (where particles are represented by black pixels): ')
disp(af)


    
hist(P,2500), xlabel('# of particle pixels in a box'), ylabel('# of occurrences')





 
 

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Manual Mode %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
 logi = input('Automatic iteration completed. Continue to manual mode? (Y/N)  :', 's');
 
while strcmpi(logi, 'y') == 1;
    xf1 = input('Enter desired square width in nanometers:  ');
    xf = xf1/scale1;
    Nsquares = input('Enter number of random squares to use for analysis:  ');

     s2=size(A,1)*size(A,1);
     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%% Calculate box size for manual input %%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    length = ceil(xf*size(A,2));
    for i = 1:Nsquares;
     x=ceil(rand*size(A,1))+offx;
     y=ceil(rand*size(A,2))+offy;
      counter=0;
      for s = 1:length;
          for t = 1:length;
              mx=x-floor(length/2)+s;
              my=y-floor(length/2)+t;
              if Matrix(mx,my)<150;
                  counter=counter+1;
              end
          end
      end
       P(i)=counter;
    end
    data=P';
    bimreal=mode(data);
    boxar=length^2;
         
         
    h1=hist(P,2500);
    
         if boxar==bimreal
            isbimo=1;
            h1(1,2500)=0;
            [m,n1]=max(h1);
         else
            isbimo=0;
            [m,n1]=max(h1);
         end
         
         
    
    zerocount=0;
    for j=1:size(P,2);
      if P(j)==0;
           zerocount=zerocount+1;
      end
    end

    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%% Display results %%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 
 
    hist(P,2500), xlabel('# of particle pixels in a box'), ylabel('# of occurrences')
    disp('Number of boxes with a zero black-pixel count: ')
    disp(zerocount)
    
    
    
    logi = input('Enter parameters again? (Y/N)  :', 's');
 
end


end

function [result,percents]=mode(x)
%MODE
% Finds the mode of a 2d matrix.
% [result perecents]=mode(matrix)
% where result is the mode of the matrix
% and percents is the amount of difference within the mode
% ORIGNALLY TABULATE.m by B.A. Jones
% Changes by David Li, UCSB updated: 4-8-2004

[Mo,No]=size(x);
x=reshape(x,Mo*No,1);

y = x(find(~isnan(x)))+1;

maxlevels = max(y(:));
minlevels = min(y(:));
[counts values] = hist(y,(minlevels:maxlevels)); 
total = sum(counts);

result=-1;
index=1; 
while(counts(index) ~= max(counts))
    index=index+1;    
end
result=values(index)-1; %disp(result);

percents =counts(index)/total;

end


 
 
 
 
 
 
 


