数字图像处理孔洞填充算法实现

两种解决思路

关于这个孔洞填充算法，我先是按照中规中矩的正向思路进行孔洞填充，后来机缘巧合之下得知网上有一种更为简单的思路，是反向求解的，这激起了我的好奇心，终于看懂理解了，也加深了对空洞填充算法的理解，在此记录一下。

普通正向思路

为求方便，我在这里直接上ppt原图吧。

算法实现思路

在这里小小的解释一下原理：这个算法的应用场景是“0-1-0”问题，即背景色为暗色，像素为0。孔洞前景色为亮色，像素为1。孔洞内部（即需要填充起来的部分）也是暗色，像素值为0。而我们要做的就是把010转化成011。

然后就是具体的实现思路，首先找到每个孔洞内部的一个点做起始点，然后再挨个点开始填补对应的孔洞。那么如何根据一个起始点填补一个孔洞呢，这里主要采用形态学的手段，首先利用特定的核将起始点向四周膨胀（“四周”是因为核的形状和内容），这其实就是要填补洞洞的节奏，但是不做约束的话，很快就会“多做功”或做“无用功”，比如你扩张填补的位置像素本来就是1，或者你膨胀到空洞外面去了，多填充了背景。当然前者的影响不是很大，但后者的影响是不可忽略的。所以我们扩张完后需要与A的补集做交集，目的就是只允许在孔洞内部扩张。然后在已经填补一部分的基础上，按上述思路一直填补，直到当前孔洞填满为止。

那么为什么每个孔洞都需要一个起始点呢？为什么不能一个起始点一直填，然后走遍天下填完全部的洞呢？

首先提几个注意事项：

1. 所有孔洞都属于连通分量，而整张图片像素值为0的连通分量除了空洞以外，其实还有整个背景，这一点很重要，在之后的代码段里会提到。
2. 由于公式里是与A的补集取交集，所以这个公式主要针对的是0-1-0模型，而不是1-0-1模型。所以使用公式的时候要注意自己的问题场景。

然后谈谈我对这个算法的本质的理解。这个算法其实就是实现了输入一个连通区域内点的坐标，它会将这个连通区域填补上，这样一个功能。由于每个连通区域之间是独立的，所以一个连通区域内部的点受算法限制，基本上是不会扩张到外面（即背景）的（除非孔洞边缘连一个像素宽都没有），所以必须要一个连通区域对应一个起始点。

具体代码实现

%% myfill_hole.m
% 对一个空洞进行迭代填充
function img = myfill_hole(img, hole)
    img(img==255)=1;
    [height, width] = size(img);
    X0 = zeros(height,width);
    X0(hole(1),hole(2)) = 1;    %获取X0
    mask=[0,1,0;1,1,1;0,1,0]; %膨胀操作的结构元
    SE=strel('arbitrary',mask);
    %求补集
    Ac = 1-img;
    %迭代填充
    Xk=X0;
    Xk_1=ones(size(Xk));
    while (~isequal(Xk_1,Xk))
        Xk_1=Xk;
        %具体操作
        Xk = imdilate(Xk_1,SE);
        Xk = Xk&Ac;
    end
    img = img|Xk;
end

%% find_holes.m
% 找到所有空洞位置的一个点
function myholes = find_holes(img)
    % 将图片转换为二值图像
    bw = im2bw(img);

    % 找到所有的孔洞
    holes = bwconncomp(~bw);
    
    % 计算每个连通区域的属性
    props = regionprops(holes, 'Area', 'Centroid');
    
    % 去掉面积最大的连通区域，即为空洞
    for k=1:length(props)
        data(k)=props(k).Area;
    end
    hole_idx = find([props.Area] < max(data));

    %建立元胞数组存放孔洞坐标
    myholes={};
    % 遍历每个孔洞
    
    for i = 1:length(hole_idx)
        % 找到当前孔洞的像素索引
        holePixels = holes.PixelIdxList{hole_idx(i)};

        % 随机选择一个像素作为孔洞内部的点
        randIndex = randi(length(holePixels));
        randPixel = holePixels(randIndex);

        % 将像素索引转换为坐标
        [x, y] = ind2sub(size(bw), randPixel);
        myholes{i}=[x,y];
        % 输出点的坐标
        fprintf('Hole %d: (%d, %d)\n', i, x, y);
    end
end

%% main.m
clear;clc;
%读取图像
img = imread('region-filling.tif');
% img = rgb2gray(img);   %转化为灰度图
% img = imbinarize(img);  %将灰度图转化为二值图

%显示原图
figure('name','原图'); 
imshow(img);
% 调用函数找到所有空洞位置
holes = find_holes(img);

%实现空洞填充操作
% 对每个空洞进行迭代填充
for i = 1:length(holes) 
    hole = holes{i};
    img = myfill_hole(img, hole);
end

%显示填充后的图像
figure('name','填充后的图像'); 
imshow(img); 

接下来主要讲讲我在代码编程中遇到的一些困难以及我是怎么解决的。

首先是我对于二值图像的疑惑，我们接收到的图像通常数据类型是uint型，像素值也不仅仅是0和1，而是0和255分别代表黑色和白色。所以这个时候为了方便之后处理通常需要用imbinarize(img)函数将图像转成数据类型为逻辑值0和1的图像，此时0表示黑色，1表示白色。所有矩阵和别的矩阵做集合操作（逻辑运算），如交集&和并集 | 后也都会转化成逻辑值。

然后最大的问题就是如何找到每个孔洞内部的一个点且不包括背景图这个最大的连通区域了（前面注意事项已经提过）。一开始我想用holes = bwconncomp(~bw);来找到每个连通区域，并从中任取一点做起始点，但是很快我就发现，这个连通区域也包含着背景图，所以一填充发现一片白色。后来我发现通常情况下，一张图片中的孔洞区域是比较小的，而背景对应的连通区域是最大的，所以我只需要将最大的连通区域去掉就好了，具体见上述代码。

上效果图-：

简单反向思路

后来我得知网上有一种更为简单的方法，而且填充效果也非常nice。具体代码如下：

clear;clc;
A = imread('region-filling.tif');
% A = rgb2gray(A);   %转化为灰度图
A = imbinarize(A); %二值化

figure(1);
subplot(1,2,1);
imshow(A);
title('二值化后的图');
A_ct = ~A; %A的补集
[m,n] = size(A);
X = zeros(m,n);
[k,z] = find(A==1); %找到A中灰度值为1的点的位置 
X(k(1),z(1)) = 1; %选取第一个灰度值为1的点为初始点
se = [0,1,0;1,1,1;0,1,0]; %结构元
over = 0;  %判断是否迭代结束
i = 0;    %迭代次数
while(~over)
    i = i + 1;
    Xp = X;
    X = imdilate(X,se) & A_ct;  %膨胀后与补集求交集
    if X == Xp  %二者相等，迭代结束
        over = 1;         
    end
end
A = ~X | A; %求并集
subplot(1,2,2);
imshow(A);
title('填充结果');
%figure('name','并集前');
%imshow(X);

没错，竟然这么短！而且为啥他只用求一个起始点？！而且为啥他这个起始点选的还是边缘点？！而且为啥他最后要对X取反后再求并集呢？！仔细分析之后才发现，原来他采用的是逆向思维，化复杂为简单，真的是太巧妙了！

那咱来一个个解决这个理解上的问题：

首先，为啥此处的公式在并集前需要取反呢？为了清晰起见，我把并集前的X画了出来，结果如下图。也就是说这个算法迭代得到的结果是背景图，结合前面注意事项提到的，背景图也是一个连通区域，然后再结合这个算法的本质（即根据连通区域内部的一个起始点将该连通区域填满），就不难理解了，这里比起正向的一个个填充每个孔洞的连通区域，他选择反其道而行，将整个背景看成一个大连通区域进行填充，最后在再取反与原图取并集就可以得到填充好的图像了。从填补多个到填补一个，只需要一个起始点，真是巧妙！让我想到了求逆Z变换的留数辅助定理，也是这样的逆向思路。

不过这里还有一个问题，我们都知道这个算法的一个难点就是找到连通域内部的起始点，但是这里为什么仅仅取了图像中的第一个亮点，也就是边缘点做起始点呢？

仔细体会膨胀和迭代的过程，不难发现，通常孔洞的亮值区域是有一定厚度的，所以当边缘点向四周扩张，并与A的补集取交集时，它也会受到约束不会进入孔洞内部，而是不断填满背景这个连通区域。不过，由此也可以看出它的隐患，如果孔洞的亮值区域（即0-1-0模型中的圆环1）太薄，膨胀的步长较大的话，也会存在填补到圆孔内部的隐患。

效果图如下：

两种方法的区别和对比

所以对比两种思路，本质上原理是一样的，只是解决问题的方向不同，方案一是将所有孔洞的起始点都找到再一一填补；方案二是只用填补背景，再将背景取反即可，起始点也仅用边缘点代替。很明显，后者的思路更简洁，实现起来也更简单，但是由于边缘点做起始点所以存在填补失误的隐患，比较适合形状比较规则，厚度均匀的孔洞。

PS：上述很多都是基于本人的理解，并没有进行严格证明，如有错误和疏忽的地方，欢迎指正！
第二种思路的代码来源于网络，如有侵权，请联系删除。