How to segment concavity based binary image?

Question

I have a binary image, for example:

I want to separate the main large elliptical white section from the small cloud of mushrooms on top. This should be an automatic process for many different images, which may be completely different, but will still have this characteristic of the main blob and a touching smaller drop from above or from the side.

I thought about using the watershed, but it does not work in all cases, depending on the proportions of the additional blob. Now I'm trying to figure out if there is a way to find the edge of the binary image and put a condition on the concavity of this edge, but I can not find how to do it.

Ideal implementation in MATLAB, but if it can be done better in SciPy / Mathematica, that’s good too.

Answer 1

. , , , , " " , , .

clear; clc;

binary_img = imread('bin.jpg') > 100;

% Get boundaries
b = bwboundaries(binary_img);

% Get largest boundary
b = b{cellfun(@length,b) == max(cellfun(@length,b))};

% Filter boundary - use circular convolution
b(:,1) = cconv(b(:,1),fspecial('gaussian',[1 81],40)',size(b,1));
b(:,2) = cconv(b(:,2),fspecial('gaussian',[1 81],40)',size(b,1));

% Find curvature
curv_vec = zeros(length(b),1);
for i = 0:size(b,1)-1
    p_b = b(mod(i-25,length(b))+1,:); % p_b = point before
    p_m = b(mod(i,length(b))+1,:);    % p_m = point middle
    p_a = b(mod(i+25,length(b))+1,:); % p_a = point after

    dx_ds = p_a(1)-p_m(1);              % First derivative
    dy_ds = p_a(2)-p_m(2);              % First derivative
    ddx_ds = p_a(1)-2*p_m(1)+p_b(1);    % Second derivative
    ddy_ds = p_a(2)-2*p_m(2)+p_b(2);    % Second derivative
    curv_vec(i+1) = dx_ds*ddy_ds-dy_ds*ddx_ds;
end

% Find local maxima for curvature
[pks,locs] = findpeaks(curv_vec);
[pks,pks_idx] = sort(pks);

% Select two largest curvatures
p1_max = b(curv_vec == pks(end),:);
p2_max = b(curv_vec == pks(end-1),:);

% Paint biggest contiguous region
rp = regionprops(binary_img,'Area','PixelIdxList','PixelList');
rp = rp(max(vertcat(rp.Area)) == vertcat(rp.Area));

% Paint all points to the left of the line
img = zeros(size(binary_img));
img(rp.PixelIdxList) = 0.5;
for i = 1:length(rp.PixelList)     
    turn = sign(det([1 p1_max(1) p1_max(2);
                     1 p2_max(1) p2_max(2);
                     1 rp.PixelList(i,2) rp.PixelList(i,1);]));

    if (turn > 0) 
        img(rp.PixelList(i,2),rp.PixelList(i,1)) = 1;
    end
end

figure(1);
subplot(1,3,1);
plot(b(:,1), b(:,2),'o');
hold on;
plot(p1_max(1), p1_max(2),'ro','Markersize',5,'LineWidth', 5);
plot(p2_max(1), p2_max(2),'ro','Markersize',5,'LineWidth', 5);

subplot(1,3,2);
plot(curv_vec);

subplot(1,3,3);
imshow(img);

:

:

Answer 2

. , : 10-15 .

, (, ), , .

, , - .

close all; clear all;
img = imread('mushroom .jpg');

% Threshold
binary = img > 100;

% Apply open morphology operator to "enlarge" holes and remove small blobs
se = strel('disk',3);        
opened= imopen(binary, se);

% Compute vertical projection
projectionVer = sum(opened,1);

% Find max value on vertical projection
% Is basically the vertical simmetry axis of the mushroom
[~, centerX] = max(projectionVer);

% Find limits of the mushroom
% This procedure can be helpful if the central blob (the mushroom)
% is linked with the small left and right blobs.
% In this case simply taking the larget boundary, or the largest
% blob will fail. 
% But since you said "touching smaller blobs either above of to the
% sides"...

leftMin = min(projectionVer(1 : centerX));
rightMin = min(projectionVer(centerX+1 : end));
leftX = find(projectionVer(1 : centerX) == leftMin, 1, 'last');
rightX = find(projectionVer(centerX+1 : end) == rightMin, 1, 'first');
rightX = centerX + rightX;

% Crop the image to keep only the mushroom
mushroom = img(:, leftX : rightX);

% Compute horizontal projection on mushroom
projectionHor = sum(mushroom, 2);

% Find first minimum peak
[pks, loc] = findpeaks(- projectionHor);
minY = loc(1); % << You are looking for this!
minYVal = -pks(1);

% Segmentation
topY = find(projectionHor>0, 1);

result = uint8(opened);

% probably you can do better than for loops, but ok for now...
for y=leftX:rightX
    for x=topY:size(result ,1)
        if(opened(x,y))            
            if(x<=minY)
                %top
                result(x,y) = 127;
            else
                %bottom
                result(x,y) = 200;
            end
        end
    end
end

%Plotting

imshow(result);
figure();

subplot(221);
imshow(img);
title('Image');

subplot(222);
hold on;
plot(flip(projectionHor), (1 : length(projectionHor)));
plot(minYVal, size(img,1) - minY, 'or');
title('Horizontal Projection');
axis([0, +Inf, 0, size(img,1)]);
hold off;

subplot(223);
hold on;
plot(projectionVer);
plot(leftX, leftMin, 'or');
plot(rightX, rightMin, 'or');
title('Vertical Projection');
axis([0, size(img,2), 0, +Inf]);
hold off;

subplot(224);
imshow(img);
hold on;
plot((1:size(img,2)), ones(1,size(img,2))*minY, 'r'); 
plot(ones(1,size(img,1))*leftX, (1:size(img,1)), 'g');
plot(ones(1,size(img,1))*rightX, (1:size(img,1)), 'b');
title('Result');
hold off;