Matlab-图像分割之区域生长法

区域增长是一种区域分割方法。假定区域的数目以及在每个区域中种子点的位置已知，则从该已知点开始，加上与已知点相似的邻近点形成一个区域。

(1)相似性准则可以是灰度级、彩色、组织、梯度或其它特性；

(2)相似性的测度可以由所确定的阈值来判定。

(3)从满足检测准则的点开始，在各个方向上生长区域，如果其邻近点满足检测准则就并入区域中；

(4)重复这一过程，直到没有可接受的邻近点后，区域增长过程终止。

区域生长法的其优点是基本思想相对简单，通常能将具有相同特征的联通区域分割出来，并能提供很好的边界信息和分割结果。在没有先验知识可以利用时，可以取得最佳的性能，可以用来分割比较复杂的图象，如自然景物。但是，区域生长法是一种迭代的方法，空间和时间开销都比较大，噪声和灰度不均一可能会导致空洞和过分割，并在对图像中的阴影效果处理上往往不是很好。

代码实现

function [seg, seg_arr, reg_info] = region_grow1(im, seg, criterion, connectivity)
% 输入参数:
% im: 用于分区的图像
% seg: 图像初始标签矩阵，与im大小相同。通常它的每个像素都有不同的标签值
% criterion: 判断区域合并与否的条件值
% connectivity: 邻域连通的大小，取4或8
% 输出参数:
% seg: 分割后的图像
% seg_arr:  每次循环处理的图像分割的结果，是中间图像
% reg_info: 记录各个分割区域的信息的结构
   
rows = size(im, 1); %计算图像行数
cols = size(im, 2); %计算图像列数
reg_info = repmat(struct('size', 0, ...             %记录图像大小
                         'crt', 0, ...              %记录规则
                         'row_coord', [], ...       %记录区域行坐标
                         'col_coord', []), 1, rows*cols);   %记录区域列坐标
 
%从初始划分中收集已有区域的信息
for m = 1:rows
  for n = 1:cols
    c_lab = seg(m, n);
    reg_info(c_lab).size = reg_info(c_lab).size + 1;
 
    %计算每一个区域的均值
    reg_info(c_lab).crt = double(...
            ((reg_info(c_lab).size-1) * reg_info(c_lab).crt + im(m,n)) / ...
            reg_info(c_lab).size);
    reg_info(c_lab).row_coord = [reg_info(c_lab).row_coord m];
    reg_info(c_lab).col_coord = [reg_info(c_lab).col_coord n];
  end
end
 
changed = 1;    %假设第一次循环区域划分会发生变化
iter = 0;       %记录循环次数
 
%图像发生变化时，按照规则继续分割图像
while(changed)    
    disp(['Currently in iteration: ' num2str(iter+1)]);
    changed = 0;
    for m = 1:rows
        for n = 1:cols
            c_lab=seg(m, n);    %获得当前像素的标签值
  
            %获取当前像素的4-连通或8-连通邻域的像素坐标 
            if (connectivity == 4)
                ngb = get4ngb(rows, cols, m, n);
            end
            if (connectivity == 8)
                ngb = get8ngb(rows, cols, m, n);
            end
 
            %检查所有邻域像素的标签值，按照规则依次判断是否与当前像素同区
            for k = 1:size(ngb, 2)
                n_lab = seg(ngb(1, k), ngb(2, k));
                if(n_lab ~= c_lab)
                    if(abs(reg_info(n_lab).crt - reg_info(c_lab).crt) < criterion)
                         
                        %更改标签值等于c_lab的区域的均值
                        reg_info(c_lab).crt = (reg_info(c_lab).size * reg_info(c_lab).crt + ...
                                               reg_info(n_lab).size * reg_info(n_lab).crt) / ...
                                              (reg_info(c_lab).size + reg_info(n_lab).size);
 
                        %更改该区域的大小
                       reg_info(c_lab).size = reg_info(c_lab).size + reg_info(n_lab).size;
 
                        %给该区域中的所有像素赋标签值c_lab 
                        tmp_row_coord = reg_info(n_lab).row_coord;
                        tmp_col_coord = reg_info(n_lab).col_coord;
                        for p = 1:reg_info(n_lab).size
                            seg(tmp_row_coord(p), tmp_col_coord(p)) = c_lab;
                        end
 
                        %更改该区域内所有像素的坐标
                        reg_info(c_lab).row_coord = [reg_info(c_lab).row_coord reg_info(n_lab).row_coord];
                        reg_info(c_lab).col_coord = [reg_info(c_lab).col_coord reg_info(n_lab).col_coord];
                        reg_info(n_lab).size = 0;
                        changed = 1;
                    end
                end
            end
        end
    end
    iter = iter+1;
     
    %将本次循环结果存在放到seg_arr中，用于验证。
%其中第一个图像是初始的分割图像
    seg_arr(:, :, iter) = seg;
end
 
%清空reg_info
index = 1;
for m = 1:max(max(seg))
    if (reg_info(m).size > 0)
        tmp_reg_info(index) = reg_info(m);
        index = index+1;
    end
end
reg_info = tmp_reg_info;

function ngb = get4ngb(rows,cols,x,y)
  % 用来获取矩阵中一个元素的4-连通相邻元素的坐标，结果存放在ngb中。
  % 如果(x,y)是矩阵边界或顶角元素，获取的将是实际的相邻元素的坐标。如边界
  % 像素有3个相邻元素，而顶角元素有2个。 
if ((x ==  1) & (y == 1))       % 处理左上角
    ngb(1:2,1:1) = [2 1]';
    ngb(1:2,2:2) = [1 2]';   
  elseif ((x == rows) & (y == 1))   % 处理左下角    
    ngb(1:2,1:1) = [rows 2]';
    ngb(1:2,2:2) = [(rows - 1) 1]'; 
  elseif (y == 1)                   % 处理左边界
    ngb(1:2,1:1) = [(x+1) 1]';
    ngb(1:2,2:2) = [x 2]';
    ngb(1:2,3:3) = [(x-1) 1]';        
  elseif ((x == 1) & (y == cols))       % 处理右上角
    ngb(1:2,1:1) = [1 (cols-1)]';
    ngb(1:2,2:2) = [2 cols]'; 
  elseif ((x == rows) & (y == cols))    % 处理右下角
    ngb(1:2,1:1) = [rows (cols-1)]';
    ngb(1:2,2:2) = [(rows-1) cols]'; 
  elseif (y == cols)                 % 处理右边界
    ngb(1:2,1:1) = [(x+1) cols]';
    ngb(1:2,2:2) = [x (cols-1)]';
    ngb(1:2,3:3) = [(x-1) cols]'; 
  elseif (x == 1)                   % 处理上边界
    ngb(1:2,1:1) = [1 (y-1)]';
    ngb(1:2,2:2) = [2 y]';
    ngb(1:2,3:3) = [1 (y+1)]'; 
  elseif (x == rows)                % 处理下边界
    ngb(1:2,1:1) = [rows (y-1)]';
    ngb(1:2,2:2) = [(rows-1) y]';
    ngb(1:2,3:3) = [rows (y+1)]'; 
  else                          %其他
    ngb(1:2,1:1) = [(x-1) y]';
    ngb(1:2,2:2) = [x (y-1)]';
    ngb(1:2,3:3) = [(x+1) y]';
    ngb(1:2,4:4) = [x (y+1)]';
  end

function ngb = get8ngb(rows,cols,x,y)
  % 用来获取矩阵中一个元素的8-连通相邻元素的坐标，结果存放在ngb中。
  % 如果(x,y)是矩阵边界或顶角元素，获取的将是实际的相邻元素的坐标。如边界
  % 像素有5个相邻元素，而顶角元素有3个。 
if ((x ==  1) & (y == 1))       % 处理左上角
    ngb(1:2,1:1) = [2 1]';
    ngb(1:2,2:2) = [2 2]';
    ngb(1:2,2:2) = [1 2]';   
  elseif ((x == rows) & (y == 1))   % 处理左下角    
    ngb(1:2,1:1) = [rows 2]';
ngb(1:2,2:2) = [(rows - 1) 2]';
    ngb(1:2,2:2) = [(rows - 1) 1]'; 
  elseif (y == 1)                   % 处理左边界
    ngb(1:2,1:1) = [(x+1) 1]';
    ngb(1:2,2:2) = [(x+1) 2]';
    ngb(1:2,2:2) = [x 2]';
    ngb(1:2,4:4) = [(x-1) 2]';
    ngb(1:2,3:3) = [(x-1) 1]';        
  elseif ((x == 1) & (y == cols))       % 处理右上角
    ngb(1:2,1:1) = [1 (cols-1)]';
    ngb(1:2,2:2) = [2 (cols-1)]';
    ngb(1:2,2:2) = [2 cols]'; 
  elseif ((x == rows) & (y == cols))    % 处理右下角
    ngb(1:2,1:1) = [rows (cols-1)]';
    ngb(1:2,2:2) = [(rows-1) (cols-1)]';
    ngb(1:2,2:2) = [(rows-1) cols]'; 
  elseif (y == cols)                 % 处理右边界
    ngb(1:2,1:1) = [(x+1) cols]';
    ngb(1:2,2:2) = [(x+1) (cols-1)]';
    ngb(1:2,3:3) = [x (cols-1)]';
    ngb(1:2,4:4) = [(x-1) (cols-1)]';
    ngb(1:2,5:5) = [(x-1) cols]';
  elseif (x == 1)                   % 处理上边界
    ngb(1:2,1:1) = [1 (y-1)]';
    ngb(1:2,2:2) = [2 (y-1)]';
    ngb(1:2,3:3) = [2 y]';
    ngb(1:2,4:4) = [2 (y+1)]';
    ngb(1:2,5:5) = [1 (y+1)]';
  elseif (x == rows)                % 处理下边界
    ngb(1:2,1:1) = [rows (y-1)]';
    ngb(1:2,2:2) = [(rows-1) (y-1)]';
    ngb(1:2,3:3) = [(rows-1) y]';
    ngb(1:2,4:4) = [(rows-1) (y+1)]';
    ngb(1:2,5:5) = [rows (y+1)]';
  else                          %其他
    ngb(1:2,1:1) = [(x-1) y]';
    ngb(1:2,2:2) = [(x-1) (y-1)]';
    ngb(1:2,3:3) = [x (y-1)]';
    ngb(1:2,4:4) = [(x+1) (y-1)]';
    ngb(1:2,5:5) = [(x+1) y]';
    ngb(1:2,6:6) = [(x+1) (y+1)]';
    ngb(1:2,7:7) = [x (y+1)]';
    ngb(1:2,8:8) = [(x-1) (y+1)]';
  end

P = imread('cameraman.jpg');
[row, col] = size(P);
i = 1;
 
for con = 4:4:8
    for crt = 2:2:12
        seg = reshape(1:row*col, row, col);  % 生成初始标签矩阵
        [seg,seg_arr,reg_info] = region_grow1(P, seg, crt, con);
        seg = im2uint8(uint16(seg));  % 将double类型数据转换成uint16类型后映射到uint8类型
        subplot(2,6,i),imshow(seg);
        title([num2str(con),'邻域下',' 均值之差<',num2str(crt)])
        i = i+1;
    end
end