isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if isOctave
    pkg load image
    addpath octave_functions
end

close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0
    

% Definisco una maschera
BW = logical([1 1 1 0 0 0 0 0
              1 1 1 0 1 1 0 0
              1 1 1 0 1 1 0 0
              1 1 1 0 0 0 1 0
              1 1 1 0 0 0 1 0
              1 1 1 0 0 0 1 0
              1 1 1 0 0 1 1 0
              1 1 1 0 0 0 0 0]);
figure('color', 'cyan'), imshow(BW) % Comportamento differente Matlab/Octave (dimensioni finestra)
figure, spy(BW)
figure, imagesc(BW), colormap([0 0 0; 1 1 1]) % Comportamento lievemente differente Matlab/Octave

L = bwlabel(BW, 4);   % ,4 oppure ,8 (default)
s = regionprops(L, 'basic');  % area, centroid, bbox (array di strutture)
% Accedo a un campo della i-ma struttura con s(i).campo:
s(1).Centroid
s(3).BoundingBox

% notare che la boubding box parte dal mezzo pixel che precede il primo
% pixel dell'oggetto, e termina mezzo pixel dopo l'ultimo pixel
% dell'oggetto, per cui le dimensioni sono consistenti
%
%    0.5    1   1.5    2    2.5   3   3.5
%      |----O----||----O----||----O----|         

disp ('Premere invio')
pause
close all

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vogliamo contare e caratterizzare gli oggetti presenti nell'immagine
% blobs.gif

[I, map] = imread('blobs.gif');
figure, imshow(I, [])
I = ind2gray(I,map);
figure, imshow(I)
impixelinfo

figure, hist(double(I(:)), 50)

disp ('Premere invio')
pause  % per dare la possibilita' di esplorare l'immagine e individuare la soglia

IM = I < 130; % 130;  % buon valore per evidenziare i blob
figure, imshow(IM)

IM = imclearborder(IM); 
figure, imshow(IM)
IM = imfill(IM, 'holes');  % Utile per soglie basse (e.g. 100)
figure, imshow(IM)

I2 = I .* uint8(IM);
figure, imshow(I2, [])

L = bwlabel(IM);
figure, imshow(L, []) 

impixelinfo

s = regionprops(L, 'basic');  % area, centroid, bbox
length(s)

s(4).Area
areas = [s.Area];   % <------------

figure, hist(areas, 20), title('Areas')

% Cerchiamo oggetti di grandi dimensioni
%idx = find(areas > 700);
idx = find(areas > 500);
L1 = L == idx(1); % idx(1) e' l'indice del primo oggetto "grosso"; L1 = L == idx(1) accende nel risultato i pixel di quest'oggetto
figure, imshow(L1) % conserva l'oggetto "grosso"; potrei cercare l'oggetto di area massima e visualizzarlo

% Marchio con un + verde gli oggetti di area > 500
figure, imshow(I, [])
hold on
for j = 1 : length(idx)
    c = s(idx(j)).Centroid % j esplora idx, idx(j) e' il j-mo indice di oggetto > 500, c ne e' il centroide
    plot(c(1), c(2), '+g')
end
hold off

% senza loop
figure, imshow(I, [])
cs = [s(idx).Centroid];  % confrontare con s(3).Centroid
x = cs(1:2:end);
y = cs(2:2:end);
hold on
plot(x, y, 'sg')
hold off

%{
% marchio con un quadrato tutti gli oggetti, senza loop
cs = [s.Centroid];
x = cs(1:2:end);
y = cs(2:2:end);
hold
plot(x, y, 'sg')
hold
%}

% dall'help:
%centroids = cat(1, s.Centroid);
%hold on
%plot(centroids(:,1), centroids(:,2), 'b*')
%hold off

% Conservare solo oggetti "abbastanza grandi", area > 50

s = s(areas > 50);   % s = s(find(areas > 50))
% gli indici di s NON CORRIPONDONO piu' alle etichette di L!!!!

% Se desidero rimettere in fase L e s, dapprima indviduo gli indici da
% eliminare:
idx = find(areas > 50);
% ismember([1 2 3 4 ; 5 6 7 8], [3 4 8 11])

L1111 = ismember(L, idx); % maschera in cui sono accesi solo i pixel degli oggetti il cui indice e' in idx
figure, imshow(L1111, [])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% falsi positivi e falsi negativi


s = regionprops(L, 'Area', 'Centroid', 'Eccentricity', 'Perimeter');  % area, centroid, etc
% notare: non posso accoppiare basic con altre opzioni!!

% Eccentricity, oltre all'help vedere:
% https://stackoverflow.com/questions/1532168/what-are-the-second-moments-of-a-region
% http://cambridge-archive.blogspot.com/2014/06/how-does-matlab-calculate-eccentricity.html

% Da qui in poi, Matlab e Octave divergono significativamente!

% CIRCULARITY (shape factor)
% https://en.wikipedia.org/wiki/Shape_factor_(image_analysis_and_microscopy)
% circ = 4*pi*A/p^2
% Vale (idealmente) 1 per una circonferenza, e valori inferiori per forme
% diverse via via piu' lontane dalla circonferenza

areas = [s.Area];
perimeters = [s.Perimeter];   % In octave, tutti leggermente piu' grandi!
circularities = 4*pi*areas./perimeters.^2;
eccentricities = [s.Eccentricity]; % In octave, completamente diverse!!

figure, hist(circularities, 50), title('circularities')
figure, hist(eccentricities, 20), title('eccentricities')

% quali sono gli oggetti con circolarita' > 1.5 Ce n'e' uno perche' molto
% piccolo, quindi misura non affidabile!

idx = find(circularities > 1.5);
figure, imshow(I, [])
hold on
for j = 1 : length(idx)
    c = s(idx(j)).Centroid
    plot(c(1), c(2), '+g')
end
hold off

% quali sono gli oggetti con circolarita' > 0.9 ma < 1.5

idx = find(circularities > 0.9 & circularities < 1.5);
figure, imshow(I, [])
hold on
for j = 1 : length(idx)
    c = s(idx(j)).Centroid
    plot(c(1), c(2), '+g')
end
hold off

% che area hanno gli oggetti che stiamo conservando

areaOggettiBuoni = areas(idx);
figure, hist(areaOggettiBuoni, 10)
s = s(idx);
s(10) = [];

% scatter plot
figure('color', 'w'), plot(circularities, eccentricities, '+'), xlim([0,4]), ylim([0,4])
