如何根据凹性对二值图像进行分割?

How to segment binary image based on concavity?

我有一个二进制图像,例如:

我希望将主要的大椭圆白色部分与顶部的小蘑菇云分开。对于许多不同的图像,这需要是一个自动过程,这些图像可能完全不同,但仍将具有主要斑点和上方或侧面接触的较小斑点的特征。

我正在考虑使用分水岭,但它并不适用于所有情况,具体取决于额外斑点的比例。我现在正在尝试查看是否有一种方法可以找到二值图像的边缘并为该边缘的凹度设置条件,但我找不到如何执行此操作。

最好在 MATLAB 中实现,但如果在 SciPy/Mathematica 中可以做得更好,那也很好。

这是我的尝试。代码很粗糙,但它是关于如何使用边界和有符号曲率找到将 "mushroom top" 与 body 的其余部分分开的两个点,然后使用 turn 谓词的一般思路确定兴趣点。

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);

使用这张图片:

输出:

您可以使用水平和垂直直方图来分割蘑菇。请注意,代码非常简单:图像处理代码只有 10-15 行。

当蘑菇 blob 与其他较小的 blob 链接时,这种方法的优点是稳健(理论上,您需要在其他图像上对其进行测试),只要保持最大的 blob 或较长的轮廓就可以导致错误。

代码中的注释应该说明每一步,但如果有不清楚的地方请评论。

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;