clear all;
close all;

%I=imread('BW.jpg');
%I=imread('spiral_puzzle_bw-LR.png');
I=imread('real_bw.tif');
bw = im2bw(I, graythresh(I));

%Make an erosion to thicken the thin black lines (if needed)
%se = strel('disk',1);
%erodeBW = imerode(bw,se);
%figure, imshow(erodeBW)
%bw=erodeBW

%Analyze the B&W image and extract the particles
L = bwlabel(bw);
s = regionprops(L,'Area','PixelList','Centroid')

% sort shapes by area size
[ss,idx]=sort([s.Area]);

figure, imshow(L);
hold on
%write number on shapes (1 is the smallest area)
for k = 1:numel(s)
    c=s(idx(k)).Centroid;
    text(c(1), c(2), sprintf('%d', k), ...
        'HorizontalAlignment', 'center', ...
        'VerticalAlignment', 'middle');
end
hold off


BW3=bw.*0;
BW3color=cat(3,bw.*0,bw.*0,bw.*0);

%define colormap as a function of the number of objects
cmp=colormap(jet(numel(s)))
for kk=1:numel(s)
    BW2 = bwpropfilt(bw,'Area',[s(idx(kk)).Area s(idx(kk)).Area]);%s(kk).Area
    BW2color=cat(3,BW2.*cmp(kk,1),BW2.*cmp(kk,2),BW2.*cmp(kk,3));
    BW3=BW2+BW3;
    BW3color=BW2color+BW3color;
end

figure, imshow(BW3);
figure, imshow(BW3color);
hold on

%write number on shapes (1 is the smallest area)
for k = 1:numel(s)
    c=s(idx(k)).Centroid;
    text(c(1), c(2), sprintf('%d', k), ...
        'HorizontalAlignment', 'center', ...
        'VerticalAlignment', 'middle');
end
hold off