基于阈值分割的图像分割-最小误差-区域生长-标记分水岭（matlab代码）

1
基于阈值分割的图像分割

1.1
阈值分割基本原理

阈值分割因为具有简单、高效、实用的特点在图像分割领域得到了广泛的使用。阈值分割适用于灰度级存在明显差异情况下的前景背景分割。使用一个固定灰度值T作为阈值将前景背景区域进行分割，遍历整幅图像的像素点，灰度值小于或等于阈值T的像素点设置为1，用白色表示；将灰度值大于阈值T的像素点设置为0，用黑色表示。分割后的图像为二值化图像。具体的过程可以用以下公式描述：

阈值分割方法能否准确快速的将图像的前景背景区域分割，其中的核心问题是阈值的选取。阈值的选取有许多的方法：根据在整幅图像中使用同一阈值或者将整幅图像分割成不同区域，在不同区域内使用不同阈值可以将阈值分割方法分为全局阈值方法和局部阈值方法；根据使用阈值的个数可以将阈值分割方法分为单阈值或多阈值方法。

1.2
最小误差法（基于传统阈值分割）

1.2.1最小误差法原理

最小误差阈值分割法是根据图像中背景和目标像素的概率分布密度来实现的，其思想是找到一个阈值，并根据该阈值进行划分，计算出目标点误分为背景的概率和背景点误分为目标点的概率，得出总的误差划分概率。当总的误差划分概率最小时，便得到所需要的最佳阈值。

用表示一副大小为的图像在坐标的像素点灰度值，。图像灰度概率分布用一维直方图h(x)表示。设图像中感兴趣的目标点的灰度也为正态分布，密度为f1(x)），均值和方差分别为miu1和sigma1；图像背景点的灰度也为正态分布，密度为f1(x)，均值和方差分别为miu2 sigma2。因此，整个密度函数可看成是两个单峰密度函数的混合，即双峰密度函数，如图2所示。

图 2 双峰密度函数

1.2.2最小误差法实现步骤

图 3 最小误差阈值算法流程图

根据图20可得最小误差阈值法的算法步骤：

步骤1：假设目标和背景灰度值得密度为，。计算混合概率密度。

步骤2：选定阈值T，计算总的误差概率，对进行求导。

步骤3：根据准则函数，计算使其最小时的T，作为最佳阈值。

1.3
区域生长（基于区域信息分割）

图像分割算法的原理是指根据图像的灰度值不连续性与相似性来实现的。不连续性是基于图像的灰度不连续变化来完成图像的分割，例如图像的边缘。相似性是根据指定的准则把图像分割成相似以及不相似的区域。基于区域生长的图像分割算法就是利用了图像灰度值的相似性特征来实现的。

3.3.1区域生长法原理

区域生长算法的原理就是求图像中相似的像素的最大连通集合。把相似灰度值的像素点集合，然后依次检测该像素点周围的8个邻域像素点，把灰度值相近的点加入到之前的集合中，如果不相似则跳过该点。然后依次浏览该图像中的其他的像素点，直到图像中所有的像素点都被检查完毕为止，这样获得的区域就是按照区域生长算法得到的新区域。

令R表示整幅图像区域，那么图像分割的过程可以看作把区域R划分为n个子区域的过程，并且划分的子区域需要满足以下条件：

区域生长算法是根据待处理目标的特征把像素点或者子区域按照一定的准则聚合成更大区域，并且聚合成的更大区域也满足上述a-e的条件。

这是一个迭代过程，每一步重新计算被扩大区域的物体成员隶属关系并消除弱边界。当没有可以消除的弱边界时，区域合并过程结束。这时，图像分割也就完成。检查这个过程会使人感觉这是一个物体内部的区域不断增长直到其边界对应于物体的真正边界为止的过程。

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

1.3.2区域生长法实现步骤

区域生长算法中需要确定三个条件：生长条件，停止条件和种子点。种子点可以是一个节点，也可以是一段连续区域的集合，解决不同的问题需要根据具体的实际情况来确定种子点的多少。

区域生长的条件是根据图像中像素点的灰度值大小变化率来确定的，区域生长的停止条件固定了区域停止生长时应该满足的条件。

可见区域生长法主要由以下三个步骤组成：

1）选择合适的种子点

2）确定相似性准则（生长准则）

3）确定生长停止条件

令表f(x,y)示一个输入图像阵列；s(x,y)表示一个种子阵列，阵列中种子点位置为1，其他位置为0；Q表示在每个位置处所用的属性。假设阵列f和阵列S的尺寸相同。基于8邻域（8连接）的一个基本区域生长算法可以说明如下。

1） 在s(x,y)中寻找所有连通分量，并把每个连通分量腐蚀为一个像素；把找到的所有这种像素标记为1，把S中的所有其他像素标记为0。

2） 在坐标对(x,y)处形成图像fq：如果输入图像在该坐标处满足给定的属性Q，则令fq(x,y)=1，否则令fq(x,y)=0。

3） 令g是这样形成的图像：即把fq中为8连通种子点的所有1值点，添加到S中的每个种子点。

4） 用不同的区域标记（如1,2,3）标出g中的每个连通分量。这就是由区域生长得到的分割图像。

在区域生长算法中需要设置一个变量max_dist来保存当前图像中已经分割好区域中两灰度值差值最大的距离。算法开始执行时，依次检测图像中的像素点灰度值。设置另外一个变量来保存当前区域已经分割好的区域灰度值的均值，然后设置当前像素点的灰度值，将两者之差的绝对值赋值给D，如果D小于等于max_dist，就把该点加入到已分割的区域中，否则继续检查下一个点。

1.4
标记分水岭（基于混合信息分割）

数学形态学起源于岩相学对岩石结构的定量描述工作，近年来在数字图像处理和机器视觉领域中得到了广泛的应用，已经成为一种独特的数字图像分析方法和理论。目前，对于利用形态学进行图像分割，分水岭算法一直是研究的热点，这是由于分水岭算法能够准确地获得前景物体的边缘，为后期操作提供较好的预分割图像，并且具有许多优点诸如实现简单、具有完整闭合的单像素分割线、易于并行化处理等。

1.4.1标记符分水岭法原理

在使用分水岭算法进行图像分割的时候，如果直接应用往往会发现存在很多的局部极大值，并且这些极大值点的数目多于目标颗粒的数目，而这种现象也是导致分水岭算法出现过分割的主要原因。标记符控制的分水岭算法的主要思想是利用一些附加的方法，在目标颗粒的内部和外部分别寻找标记符，从而来引导分水岭进行分割，进而防止过分割的产生。

标记符是指属于原图像的一个连接分量，包括内部标记符集合和外部标记符集合。其中，内部标记符集合是指标记符处在每一个感兴趣对象的内部，外部标记符集合是指标记符包含在背景中，通过标记符用来修改梯度图像，计算内部和外部标记符集合常用的方法有线性滤波、非线性滤波及形态学处理等等。

基于标记控制的分水岭方法基本思想就是事先分别提取出代表目标区域的内部标记和背景区域的外部标记，然后再利用这些标记去修改原始梯度图，最后再进行分水岭算法。通过该方法极大地限制了局部极小值的个数，很好地控制了过分割现象的发生。

1.4.2标记分水岭法实现步骤

标记符控制分水岭算法的具体方法描述如下：

1） 对原始图像进行预处理操作。通过对比度增强、直方图匹配、形态学去噪等操作改善图像的显示效果，去除干扰噪声。

2） 距离变换。由于经过形态学处理后的图像为二值图像，因此首先应进行距离变换。背景图像可以用一个大小和原图像相同的二维数组来表示，记为，其中aij=0的像素点对应背景，aij=1的点对应目标。设为背景点集合，则距离变换就是对A中所有的像素点求dij：

式中，，从而得到欧式距离变换图像Dmxn。

3） 借助数学形态学中的方法，得到内部标记点及内部标记图像fgm，消除过分割。

4） 利用一次分水岭运算得到外部标记符及外部标记图像fbw；然后基于获取到的内外标记，利用强制最小技术对原梯度图像进行修改，修改后的梯度图像可表示为：

5） 接下来对图像进行分水岭变换，最终得到理想的分割效果，分水岭变换可表示为：

标记分水岭

function Watershed_Fun(fileName)

rgb = imread(fileName);
if ndims(rgb) == 3
    I = rgb2gray(rgb);
else    
    I = rgb;
end
sz = size(I);
if sz(1) ~= 256
    I = imresize(I, 256/sz(1));
    rgb = imresize(rgb, 256/sz(1));
end
hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
se = strel('disk', 3);
Io = imopen(I, se);
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
Ioc = imclose(Io, se);
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
fgm = imregionalmax(Iobrcbr);
se2 = strel(ones(3,3));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
fgm4 = bwareaopen(fgm3, 15);
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;
gradmag2 = imimposemin(gradmag, bgm | fgm4);
L = watershed(gradmag2);
Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');

[pathstr, name, ext] = fileparts(fileName);
filefolder = fullfile(pwd, '实验结果', [name, '_实验截图']);
if ~exist(filefolder, 'dir')
    mkdir(filefolder);
end
h1 = figure(1);
set(h1, 'Name', '图像灰度化', 'NumberTitle', 'off');
subplot(1, 2, 1); imshow(rgb, []); title('原图像');
subplot(1, 2, 2); imshow(I, []); title('灰度图像');
fileurl = fullfile(filefolder, '1');
set(h1,'PaperPositionMode','auto');
print(h1,'-dtiff','-r200',fileurl);
h2 = figure(2);
set(h2, 'Name', '图像形态学操作', 'NumberTitle', 'off');
subplot(1, 2, 1); imshow(Iobrcbr, []); title('图像形态学操作');
subplot(1, 2, 2); imshow(bw, []); title('图像二值化');
fileurl = fullfile(filefolder, '2');
set(h2,'PaperPositionMode','auto');
print(h2,'-dtiff','-r200',fileurl);
h3 = figure(3);
set(h3, 'Name', '图像梯度显示', 'NumberTitle', 'off');
subplot(1, 2, 1); imshow(rgb, []); title('待处理图像');
subplot(1, 2, 2); imshow(gradmag, []); title('梯度图像');
fileurl = fullfile(filefolder, '3');
set(h3,'PaperPositionMode','auto');
print(h3,'-dtiff','-r200',fileurl);
h4 = figure(4); imshow(rgb, []); hold on;
himage = imshow(Lrgb);
set(h4, 'Name', '图像分水岭分割', 'NumberTitle', 'off');
set(himage, 'AlphaData', 0.3);
hold off;
fileurl = fullfile(filefolder, '4');
set(h4,'PaperPositionMode','auto');
print(h4,'-dtiff','-r200',fileurl);

最小误差

clc; clear all; close all;
warning off all;
% 读取图像
filename = fullfile(pwd, 'images/00.jpg');
Img = imread(filename);
% 灰度化
if ndims(Img) == 3
    I = rgb2gray(Img);
else
    I = Img;
end
% 直接二值化
bw_direct = im2bw(I);
figure; imshow(bw_direct); title('直接二值化分割');
% 圈选胃区域空气
c = [1524 1390 1454 1548 1652 1738 1725 1673 1524];
r = [1756 1909 2037 2055 1997 1863 1824 1787 1756];
bw_poly = roipoly(bw_direct, c, r);
figure;
imshow(I, []);
hold on;
plot(c, r, 'r-', 'LineWidth', 2);
hold off;

% 
J = I;
J(bw_poly) = 255;
% 图像增强
J = mat2gray(J);
J = imadjust(J, [0.532 0.72], [0 1]);
J = im2uint8(mat2gray(J));
figure; imshow(J, []); title('图像增强处理');
% 直方图统计
[counts, gray_style] = imhist(J);
% 亮度级别
gray_level = length(gray_style);
% 计算各灰度概率
gray_probability  = counts ./ sum(counts);
% 统计像素均值
gray_mean = gray_style' * gray_probability;
% 初始化
gray_vector = zeros(gray_level, 1);
w = gray_probability(1);
mean_k = 0;
gray_vector(1) = realmax;
ks = gray_level-1;
for k = 1 : ks
    % 迭代计算
    w = w + gray_probability(k+1);
    mean_k = mean_k + k * gray_probability(k+1);
    % 判断是否收敛
    if (w < eps) || (w > 1-eps)
        gray_vector(k+1) = realmax;
    else
        % 计算均值
        mean_k1 = mean_k / w;
        mean_k2 = (gray_mean-mean_k) / (1-w);
        % 计算方差
        var_k1 = (((0 : k)'-mean_k1).^2)' * gray_probability(1 : k+1);
        var_k1 = var_k1 / w;
        var_k2 = (((k+1 : ks)'-mean_k2).^2)' * gray_probability(k+2 : ks+1);
        var_k2 = var_k2 / (1-w);
        % 计算目标函数
        if var_k1 > eps && var_k2 > eps
            gray_vector(k+1) = 1+w * log(var_k1)+(1-w) * log(var_k2)-2*w*log(w)-2*(1-w)*log(1-w);
        else
            gray_vector(k+1) = realmax;
        end
    end
end
% 极值统计
min_gray_index = find(gray_vector == min(gray_vector));
min_gray_index = mean(min_gray_index);
% 计算阈值
threshold_kittler = (min_gray_index-1)/ks;
% 阈值分割
bw__kittler = im2bw(J, threshold_kittler);
% 显示
figure; imshow(bw__kittler, []); title('最小误差法分割');
% 形态学后处理
bw_temp = bw__kittler;
% 反色
bw_temp = ~bw_temp;
% 填充孔洞
bw_temp = imfill(bw_temp, 'holes');
% 去噪
bw_temp = imclose(bw_temp, strel('disk', 5));
bw_temp = imclearborder(bw_temp);
% 区域标记
[L, num] = bwlabel(bw_temp);
% 区域属性
stats = regionprops(L);
Ar = cat(1, stats.Area);
% 提取目标并清理
[Ar, ind] = sort(Ar, 'descend');
bw_temp(L ~= ind(1) & L ~= ind(2)) = 0;
% 去噪
bw_temp = imclose(bw_temp, strel('disk',20));
bw_temp = imfill(bw_temp, 'holes');
figure;
subplot(1, 2, 1); imshow(bw__kittler, []); title('待处理二值图像');
subplot(1, 2, 2); imshow(bw_temp, []); title('形态学后处理图像');
% 提取肺边缘
ed = bwboundaries(bw_temp);
% 显示
figure;
subplot(2, 2, 1); imshow(I, []); title('原图像');
subplot(2, 2, 2); imshow(J, []); title('增强图像');
subplot(2, 2, 3); imshow(bw_temp, []); title('二值化图像');
subplot(2, 2, 4); imshow(I, []); hold on;
for k = 1 : length(ed)
    % 边缘
    boundary = ed{k};
    plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2);
end
title('边缘显示标记');
figure;
subplot(1, 2, 1); imshow(bw_temp, []); title('二值图像');
subplot(1, 2, 2); imshow(I, []); hold on;
for k = 1 : length(ed)
    % 边缘
    boundary = ed{k};
    plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2);
end
title('边缘显示标记');

区域生长

clc; clear all; close all;
I = imread(fullfile(pwd, 'images/00.jpg'));
X = imadjust(I, [0.2 0.8], [0 1]);
% 阈值分割
bw = im2bw(X, graythresh(X));
[r, c] = find(bw);
rect = [min(c) min(r) max(c)-min(c) max(r)-min(r)];
Xt = imcrop(X, rect);
% 自动获取种子点
seed_point = round([size(Xt, 2)*0.15+rect(2) size(Xt, 1)*0.4+rect(1)]);
% 区域生长分割
X = im2double(im2uint8(mat2gray(X)));
X(1:rect(2), :) = 0;
X(:, 1:rect(1)) = 0;
X(rect(2)+rect(4):end, :) = 0;
X(:, rect(1)+rect(3):end) = 0;
[J, seed_point, ts] = Regiongrowing(X, seed_point);
figure(1);
subplot(1, 2, 1); imshow(I, []);
hold on;
plot(seed_point(1), seed_point(2), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
title('种子点自动选择');
hold off;
subplot(1, 2, 2); imshow(J, []); title('区域生长图像');
% 形态学后处理
bw = imfill(J, 'holes');
bw = imopen(bw, strel('disk', 5));
% 提取边缘
ed = bwboundaries(bw);
figure;
subplot(1, 2, 1); imshow(bw, []); title('形态学后处理图像');
subplot(1, 2, 2); imshow(I); 
hold on;
for k = 1 : length(ed)
    % 边缘
    boundary = ed{k};
    plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2);
end
hold off;
title('边缘标记图像');

function [J, seed_point, ts] = Regiongrowing(I, seed_point)
% 统计耗时
t1 = cputime;
% 参数检测
if nargin < 2
    % 显示并选择种子点
    figure; imshow(I,[]);  hold on;
    seed_point = ginput(1);   
    plot(seed_point(1), seed_point(2), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
    title('种子点选择');
    hold off;
end
% 变量初始化
seed_point = round(seed_point);
x = seed_point(2);
y = seed_point(1);
I = double(I);
rc = size(I);
J = zeros(rc(1), rc(2));
% 参数初始化
seed_pixel = I(x,y);
seed_count = 1;
pixel_free = rc(1)*rc(2); 
pixel_index = 0;
pixel_list = zeros(pixel_free, 3);
pixel_similarity_min = 0;
pixel_similarity_limit = 0.1;
% 邻域
neighbor_index = [-1 0;
        1 0;
        0 -1;
        0 1];
% 循环处理
while pixel_similarity_min < pixel_similarity_limit && seed_count < rc(1)*rc(2)   
    % 增加邻域点
    for k = 1 : size(neighbor_index, 1)
        % 计算相邻位置
        xk = x + neighbor_index(k, 1); 
        yk = y + neighbor_index(k, 2);
        % 区域生长
        if xk>=1 && yk>=1 && xk<=rc(1) && yk<=rc(2) && J(xk,yk) == 0
            % 满足条件
            pixel_index = pixel_index+1;
            pixel_list(pixel_index,:) = [xk yk I(xk,yk)]; 
            % 更新状态
            J(xk, yk) = 1;
        end
    end
    % 更新空间
    if pixel_index+10 > pixel_free 
        pixel_free = pixel_free+pixel_free;
        pixel_list(pixel_index+1:pixel_free,:) = 0;
    end
    % 统计迭代
    pixel_similarity = abs(pixel_list(1:pixel_index,3) - seed_pixel);
    [pixel_similarity_min, index] = min(pixel_similarity);
    % 更新状态
    J(x,y) = 1; 
    seed_count = seed_count+1;
    seed_pixel = (seed_pixel*seed_count + pixel_list(index,3))/(seed_count+1);
    % 存储位置
    x = pixel_list(index,1); 
    y = pixel_list(index,2);
    pixel_list(index,:) = pixel_list(pixel_index,:);
    pixel_index = pixel_index-1;
end
% 返回结果
J = mat2gray(J);
J = im2bw(J, graythresh(J));
% 统计耗时
t2 = cputime;
ts = t2 - t1;