function [concave] = freeman_concave(I)
%% Find concave points of object's boundary

B = im2bw(I,graythresh(I));
[M,N] = size(B);
% 8-connectivity, counterclockwise
b = boundaries(B, 8, 'cww');  
d = cellfun('length', b);
[max_d, k] = max(d);
b = b{1};
bim = bound2im(b, M, N, min(b(:,1)), min(b(:,2)));
c = fchcode(b);
% imshow(bim);
% text(c.x0y0(2), c.x0y0(1), 'starting point', 'Color', 'r', 'LineWidth', 16);
 
%%
l = length(c.fcc);

% relative chain code
R(1) = mod((c.fcc(1)-c.fcc(l)+8),8);
for i = 2:l
    R(i) = mod((c.fcc(i)-c.fcc(i-1)+8),8);
end
n = find(R>4);
R(n) = R(n)-8;

% absolute chain code
A0 = 0;
A(1) = R(1);
for i = 2:l
    A(i) = A(i-1)+R(i);
end

% chain code sum
Sum0 = A0+A(l-1)+A(l-2)+16;
Sum(1) = A(1)+A0+A(l-1)+8;
Sum(2) = A(2)+A(1)+A0;
for i = 3:l
    Sum(i) = A(i)+A(i-1)+A(i-2);
end
Sum = Sum-Sum0;
Sum0 = 0;

% chain code difference
D0 = Sum(3)-Sum0;
for i = 1:l-3
    D(i) = Sum(i+3)-Sum(i);
end
D(l-2) = Sum(1)-24-Sum(l-2);
D(l-1) = Sum(2)-24-Sum(l-1);
D(l) = Sum(3)-24-Sum(l);

Diff(1) = D0;
for i = 2:l
    Diff(i) = D(i-1);
end


%%
% concave detection
p = find(Diff>=3 & Diff<=10);
concave = b(p,:);

% figure,imshow(bim);
% hold on
% plot(concave(:,2), concave(:,1), 'r*');

end






    
















 























