图像形态学处理——形态学区域填充（Matlab）

clc;
clear all;
close all;
%区域填充测试图像
I=im2double(imread('D:\Gray Files\9-16.tif'));
%获得图像大小
[M,N]=size(I);

%----------------------------形态学区域填充--------------------------------
%===============================边界提取===================================
%结构元素
n=3;
B=ones(n,n);
n_B=length(find(B==1));
%这里需要B对其原点进行翻转，因为B是对称的，所以翻转后的结果与其本身相同
n_l=floor(n/2);
%腐蚀操作
%存放腐蚀后的图像
J=zeros(M,N);
I_pad=padarray(I,[n_l,n_l],'symmetric');
for x=1:M
    for y=1:N
        %从扩展图像中取出子图像
         Block=I_pad(x:x+2*n_l,y:y+2*n_l);
         %将结构元素与子图像点乘,即逻辑“与”操作
         c=B.*Block;
         %比较结构元素与c中的1的数量，如果一样多，则该点的值为1
         ind=find(c==1);
         if length(ind)==n_B
             J(x,y)=1;
         end
    end
end
Beta=I-J;
% imshow(Beta,[]);
% imwrite(Beta,'D:\Gray Files\9-14-BoundaryExtration.tif','tif');
%%
%================用8连通找出边界对象，并对符合条件的对象进行填充===========
%保存所有连通对象集
Objs={};
%对边界图进行扩充，四周各加0，目的是为了处理边界点
Beta_pad=padarray(Beta,[n_l,n_l]);
%找到所有1（白色）点
[rows,cols]=find(Beta_pad==1);
%将行列（x,y）组合成一个下标矩阵
ind=cat(2,rows,cols);
set(0,'RecursionLimit',700);
%递归寻找连通对象
while 1
    %取出下标矩阵中的第一个值
    if isempty(ind)
        break;
    end
    %取出下标矩阵中的第一个对象
    p=ind(1,:);
    A=FindConnectedPoint(p,Beta_pad,ind,n_l);
    %取出边界线对象
    Obj=A{1,1};
    %找出x,y的最大最小值
    x_max=max(Obj(:,1));
    x_min=min(Obj(:,1));
    y_max=max(Obj(:,2));
    y_min=min(Obj(:,2));
    %从原图中取出该对象区域,这里调整了坐标，因为上面对边界图进行了扩充
    Block=I(x_min-n_l:x_max-n_l,y_min-n_l:y_max-n_l);
    [m,n]=size(Block);
    %获得对象区域中1的数量
    num=length(find(Block)==1);
    if num/(m*n)<0.5   
        %将对象区域内填充
        I(x_min-n_l:x_max-n_l,y_min-n_l:y_max-n_l)=1;
        Objs=cat(1,Objs,A{1,1});
    end
    ind=A{1,2}; 
    Beta_pad=A{1,3};  
end
imshow(I)

FindConnectedPoint函数如下：

function R=FindConnectedPoint(p,Beta_pad,ind,n_l)
    %保存边界坐标对象
    Obj=[];
    %以p为中心点，从边界图中取出3*3大小的块区域
    Block=Beta_pad(p(1,1)-n_l:p(1,1)+n_l,p(1,2)-n_l:p(1,2)+n_l);
    Obj=p;
    %将图像中的p点置零
    Beta_pad(p(1,1),p(1,2))=0;
    %将块区域中的p点置零
    Block(n_l+1,n_l+1)=0;
    A=cell(1,3);
    %将p点从下标矩阵中删除
    for j=1:size(ind,1)
        if ind(j,:)==p
            ind(j,:)=[];
            break;
        end
    end
    %寻找块区域中为1的点的下标
    [rows,cols]=find(Block==1);
    ind_sub=cat(2,rows,cols);
    if ~isempty(ind_sub)
        %确定块中数值为1的点的坐标
        for i=1:size(ind_sub,1) 
            p_next=[];
            if ind_sub(i,1)<=n_l+1
                p_next(1,1)=p(1,1)-abs(n_l+1-ind_sub(i,1));
                if ind_sub(i,2)<=n_l+1                    
                    p_next(1,2)=p(1,2)-abs(n_l+1-ind_sub(i,2));
                else
                    p_next(1,2)=p(1,2)+abs(n_l+1-ind_sub(i,2));                    
                end
            else
                p_next(1,1)=p(1,1)+abs(n_l+1-ind_sub(i,1));
                if ind_sub(i,2)<=n_l+1                    
                    p_next(1,2)=p(1,2)-abs(n_l+1-ind_sub(i,2));
                else
                    p_next(1,2)=p(1,2)+abs(n_l+1-ind_sub(i,2));                    
                end                
            end
            if Beta_pad(p_next(1,1),p_next(1,2))~=0                        
                A=FindConnectedPoint(p_next,Beta_pad,ind,n_l);
                Obj=cat(1,A{1,1},Obj);   
            end                       
            %更新下标矩阵
            ind=A{1,2};
            Beta_pad=A{1,3};
        end
    end
    %返回结果
    R={Obj,ind,Beta_pad};
end