% clear all, clc;

se = strel('disk', 5);
r = 15;

[I1,MAP] = imread('C:\Users\qinz\Desktop\OneDrive\OneDrive for Business\CRC\latex\matlab\ch6\Sea_Ice_Floe_Identification\sea_ice_test.jpg');
I = rgb2gray(I1);
bw = im2bw(I', graythresh(I'));
bw1 = bw;
img_Dist = bwdist(~bw1, 'cityblock');
imgDist = -bwdist(~bw1, 'cityblock');
imgDist(~bw1) = -inf;
Dis_img = imregionalmin(imgDist);
dis = Dis_img.*bw1;
dis1 = imdilate(dis, se);

[p, q] = find(dis == 1);

imshow(dis);
figure, imshow(dis1);
figure, imshow(Dis_img);
figure, imshow(img_Dist, []);
figure, imshow(bw), hold on
plot(q, p, 'g+'), hold off;
figure, imshow(bw), hold on
[label, num] = bwlabel(dis1, 8);
t = 0:0.05:6.28;
for n = 1 : num
    cen = regionprops(label == n, 'centroid');   % center
    cen = cat(1, cen.Centroid);   % cen(:,1): horizontal, cen(:,2): vertical
    x = cen(1) + r * cos(t);
    y = cen(2) + r * sin(t);
    plot(cen(1), cen(2), 'r+');
    plot(x, y, 'b');
end

figure, imshow(bw), hold on
for n = 1 : num
    cen = regionprops(label == n, 'centroid');   % center
    cen = cat(1, cen.Centroid);   % cen(:,1): horizontal, cen(:,2): vertical
    r0 = abs(img_Dist(round(cen(2)), round(cen(1))) / sqrt(2));
    if r0 == 0
        r0 = 2;
    end
    x = cen(1) + r0 * cos(t);
    y = cen(2) + r0 * sin(t);
    plot(cen(1), cen(2), 'r+');
    plot(x, y, 'b');
end

[l, nn] = bwlabel(bw, 4);
rgb = label2rgb(l, @jet, [1 1 1], 'shuffle');
figure, imshow(rgb), hold on
for n = 1 : num
    cen = regionprops(label == n, 'centroid');   % center
    cen = cat(1, cen.Centroid);   % cen(:,1): horizontal, cen(:,2): vertical
    r0 = abs(img_Dist(round(cen(2)), round(cen(1))) / sqrt(2));
    if r0 == 0
        r0 = 2;
    end
    x = cen(1) + r0 * cos(t);
    y = cen(2) + r0 * sin(t);
    plot(cen(1), cen(2), 'r+');
    plot(x, y, 'b');
end
     