close all
clear all

I = imread('sample_4.bmp');	% reading of the picture from the file
figure, imshow (I)
 
I2 = rgb2gray(I);	% conversion of the picture to the greyscale image
figure, imshow (I2)
 
% setting of a threshold and the conversion to the black and white image
threshold = graythresh(I2);
I = im2bw(I2,threshold);

figure, imshow(I)	% depiction of the image
 
% removal of the grid lines in a horizontal position
dim = size(I2);
for i = 1:dim(1)
    j = 1;
    poc = 0;
    pom = I(i,j);
    while (j < dim(2) && pom == 0 && pom == I(i,j+1))
        j = j+1;
        poc = poc+1;
        if poc > (dim(1)/2)
            for k = 1:dim(2)
                I(i,k) = I(i-1,k);
            end
        end
    end
end
 
% removal of the grid lines in a vertical position
for j = 1:dim(2)
    i = 1;
    poc = 0;
    pom = I(i,j);
    while (i < dim(1) && pom == 0 && pom == I(i+1,j))
        i = i+1;
        poc = poc+1;
        if poc > (dim(2)/2)
            for k = 1:dim(1)
                I(k,j) = I(k,j-1);
            end
        end
    end
end

I = imcomplement (I); % inversion of black and white colours
figure, imshow(I)
 
I = bwareaopen(I,100); % removal of objects smaller than 100 pixels
 
I = imfill(I,'holes');	% filling of holes inside objects
 
[B,L] = bwboundaries(I,'noholes'); % finding of boundaries of objects
 
% coloration of objects and their boundaries
map = [0, 0, 0.7];
for k = 1:length(B)
    map = [map; [0, 0, 0.7]];
end
figure
img = label2rgb(L, map, [.5 .5 .5]);
imshow(img)
hold on
for k = 1:length(B)
  boundary = B{k};
  plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 1)
end
 
% gaining of areas and perimeters of objects
stats = regionprops(L,'Area', 'Perimeter');

% conversion from pixels to micrometers (1 pix = 1/1056 mm)
conver = (1/1056)*1000;
for k = 1:length(stats)
    stats(k).Area = stats(k).Area*(conver^2);
    stats(k).Perimeter = stats(k).Perimeter*conver;
end
 
% finding of minimal and maximal areas
min = stats(1).Area;
max = 1;
area = 0;
for k = 1:length(B)
  area(k) = stats(k).Area;
  if min > area(k)
      min = area(k);
      indexmin = k;
  elseif max < area(k)
      max = area(k);
      indexmax = k;
  end
end
 
% highlighting of the smallest and largest objet
boundary = B{indexmin};
plot(boundary(:,2), boundary(:,1), 'y', 'LineWidth', 2)
boundary = B{indexmax};
plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)
 
% calculation of diameters based on areas
for k = 1:length(B)
    da(k) = sqrt(4*area(k)/pi());
end
 
sumda = sum(da);
da10 = sumda/length(B);
 
da2 = da.^2;
sumda2 = sum(da2);
da20 = sqrt(sumda2/length(B));
 
da3 = da.^3;
sumda3 = sum(da3);
da30 = nthroot(sumda3/length(B),3);
daS = sumda3/sumda2;
 
da4 = da.^4;
sumda4 = sum(da4);
daB = sumda4/sumda3;
 
% calculation of diameters based on perimeters
for k=1:length(B)
    perimeter(k) = stats(k).Perimeter;
end
for k = 1:length(B)
    dp(k) = perimeter(k)/(2*pi());
end
 
sumdp = sum(dp);
dp10 = sumdp/length(B);
 
dp2 = dp.^2;
sumdp2 = sum(dp2);
dp20 = sqrt(sumdp2/length(B));
 
dp3 = dp.^3;
sumdp3 = sum(dp3);
dp30 = nthroot(sumdp3/length(B),3);
dpS = sumdp3/sumdp2;
 
dp4 = dp.^4;
sumdp4 = sum(dp4);
dpB = sumdp4/sumdp3;
 
% graphs
dasort = sort(da);
areasort = sort(area);
dpsort = sort(dp);

figure
subplot(2,2,1),
    plot (area,'r.')
    xlim([0 125]);
    title('A) Object''s area - no order');
    xlabel('Object''s index k [-]')
    ylabel('Area [um^2]')
subplot(2,2,3),
    plot (areasort, '.')
    xlim([0 125]);
    title('B) Object''s area - ordered');
    xlabel('Object''s index k [-]')
    ylabel('Area [um^2]')
area2 = areasort(1:90);
subplot(1,2,2),
    plot (area2, '.')
    title('C) Object''s area - ordered, detail');
    xlabel('Object''s index k [-]')
    ylabel('Area [um^2]')
    
kg = [1:length(B)];
figure
    p = plot (dasort,kg/length(B), dasort,(1-kg/length(B)));
    p(1).Marker = '.';
    p(1).MarkerSize = 8;
    p(2).Marker = '.';
    p(2).MarkerSize = 8;
    xlabel('Object''s diameter d_a [um]')
    ylabel('F(d_a), S(d_a) [-]')
    legend('F(d_a)','S(d_a)')
  
% writing of characteristics
pocet = length(B);
formatSpec_pocet = 'Number of objects: %d.';
str_pocet = sprintf(formatSpec_pocet,pocet);
    
formatSpec_min = 'The smallest area: %d um^{2}.';
str_min = sprintf(formatSpec_min,round(areasort(1)));
 
formatSpec_max = 'The largest area: %d um^{2}.';
str_max = sprintf(formatSpec_max,round(areasort(length(B))));
 
formatSpec_mean = 'The avarage area value: %d um^{2}.';
str_mean = sprintf(formatSpec_mean,round(mean(area)));
 
formatSpec_median = 'The median value: %d um^{2}.';
str_median = sprintf(formatSpec_median,round(median(area)));
 
str_da10 = sprintf('d_{a[1,0]}: %d um.', round(da10));
str_da20 = sprintf('d_{a[2,0]}: %d um.', round(da20));
str_da30 = sprintf('d_{a[3,0]}: %d um.', round(da30));
str_daS = sprintf('d_{aS}: %d um.', round(daS));
str_daB = sprintf('d_{aB}: %d um.', round(daB));
str_dp10 = sprintf('d_{p[1,0]}: %d um.', round(dp10));
str_dp20 = sprintf('d_{p[2,0]}: %d um.', round(dp20));
str_dp30 = sprintf('d_{p[3,0]}: %d um.', round(dp30));
str_dpS = sprintf('d_{pS}: %d um.', round(dpS));
str_dpB = sprintf('d_{pB}: %d um.', round(dpB));
 
figure
axis off
text (0.05,0.9, str_pocet)
text (0.05,0.8, str_min)
text (0.05,0.75, str_max)
text (0.05,0.65, str_mean)
text (0.05,0.6, str_median)
text (0.05,0.5, str_da10)
text (0.05,0.45, str_da20)
text (0.05,0.4, str_da30)
text (0.05,0.35, str_daS)
text (0.05,0.3, str_daB)
text (0.4,0.5, str_dp10)
text (0.4,0.45, str_dp20)
text (0.4,0.4, str_dp30)
text (0.4,0.35, str_dpS)
text (0.4,0.3, str_dpB)