clear all, clc;

img = imread('q.jpg');
img = im2bw(img, graythresh(img));
% h = fspecial('gaussian', 1, 3);  % gaussian, 3, 2.5
% img = imfilter(img, h, 'replicate');
bw = img;

[m, n] = size(bw);

% watershed ridge lines
gc = ~bw;
D = bwdist(gc, 'cityblock');        % distance transform
L = watershed(-D);
w = L == 0;         % watershed ridge lines of the negtive of the distance
% figure, imshow(w);

f = bitand(bw, w);
% f = bwareaopen(f, 5);           % ice floe junction lines
                
seg = bw & ~w;                  % ice floe segmentation
seg0 = seg;
figure, imshow(seg0); 


[label, num] = bwlabel(f, 4);

% End points of junction lines
wr = [0 -1 0; -1 4 -1; 0 -1 0];

figure, imshow(f);
hold on

for i = 1:num
    p = find(label == i);
    g = zeros(m, n);
    g(p) = 1;
    g2 = abs(imfilter(double(g), wr));
    T = max(g2(:));
    g2 = g2 >= T;
    [x, y] = find(g2 == 1);
    ep = [x, y];             % Ending points of junction lines
    plot(y, x, 'r.');
        
    neighbor = seg;
    neighbor(p) = 1;
        
    connect = bwlabel(neighbor);
    im = imreconstruct(g, connect);

    concave = freeman_concave(im);  % concave points
    c = intersect(ep,concave, 'rows');
    if numel(c) == 0
        seg(p) = 1;
    end
       
end
hold off

figure, imshow(seg);

































