大津算法的matlab实现

大津算法详解

一、算法功能

​ 图像分割就是把图像分成若干个特定的、具有独特性质的区域并提出感兴趣目标的技术和过程。它是由图像处理到图像分析的关键步骤。

​ 大津算法也称最大类间差法，由大津于1979年提出，被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。

二、图像的简单分类

​ 一般而言，图像分为彩色图，灰度图，二值图。在计算机处理中，一般将图片作为一个包含多个通道的像素矩阵来处理。若将图片长度计为m，图片宽度计为n，有如下：

​ 彩色图：有R,G,B三原色。三个通道，每个通道的矩阵大小为图片大小。矩阵大小为m×n×3

​ 灰度图：只有一个颜色通道，值域为0-255，矩阵大小为m×n

​ 二值图：只有一个颜色通道，值域为0或1，矩阵大小为m×n

三、大津分割算法

​ 大津算法目标为从所有的像素值中寻找类间方差最大的像素值，将其作为分割的界限，即阈值。使得比该阈值小的像素值为0，其余为1。

1.类间方差

      ω0=N0/ (M×N)             (1)
      ω1=N1/ (M×N)             (2)
　　　 N0+N1=M×N                (3)
　　　 ω0+ω1=1　　　             (4)
　　　 μ=ω0*μ0+ω1*μ1            (5)
　　　 g=ω0(μ0-μ)^2+ω1(μ1-μ)^2  (6)
将式(5)代入式(6),得到等价公式:
　　　 g=ω0ω1(μ0-μ1)^2              %类间方差方差 

2.求类间方差最大值对应的像素值

​ 首先，从图片的像素点中寻找出像素最大值与像素最小值。我们要遍历的灰度值便在这最大值与最小值之间。

Max = max(Imag(:));               % 最大值
Min = min(Imag(:));               % 最小值
T = Min:Max;                     % 灰度值范围

​ 其次，遍历这些灰度值，按照公式计算前景比例，背景比例，前景灰度平均值，背景灰度平均值。并按照这些参数计算方差。

for i = 1 : length(T)
TK = T(i);
iFg = 0;          % 前景
iBg = 0;          % 背景
FgSum = 0;    % 前景总数
BgSum = 0;    % 背景总数
W1=Imag>TK;
WK1=W1.*Imag;
iFg1=sum(W1(:));%求前景比例
iBg1=imagSize-iFg1;
FgSum1=sum(WK1(:));
BgSum1=sum(Imag(:))-FgSum1;

w0 = iFg1/imagSize;      % 前景比例
w1 = iBg1/imagSize;     % 背景比例
u0 = FgSum1/iFg1;         % 前景灰度平均值
u1 = BgSum1/iBg1;        % 背景灰度平均值
Tval(i) = w0*w1*(u0 - u1)*(u0 - u1);     % 计算方差
end

​ 得到方差矩阵后，寻找最大下标对应的像素值即为阈值。（部分代码参考(52条消息) MATLAB图像处理——大津阈值分割（附代码）_HNU_刘yuan的博客-CSDN博客_大津阈值，对其进行了优化）

[~, flag] = max(Tval);             % 最大值下标
TValue = T(flag);

3.根据得到的阈值进行图像分割

BW = im2bw(I, T/255);%转为二值图像

4.结果展示

四、完整代码

clear; clc;
F=imread('SF.jpg');%彩图
I=rgb2gray(F);%灰度图
subplot(2, 2, 1)
imshow(I);
xlabel('(a) 灰度图');

T = OTSU(double(I));     %使用大津法计算阈值
disp(['大津法计算灰度阈值:', num2str(T)])
BW = im2bw(I, T/255);%转为二值图像
%阈值分割
subplot(2, 2, 2)
imshow(BW);
xlabel('(b) 大津法');
 
subplot(2, 2, 3)
imshow(F);
xlabel('(c) 彩色图');

function TValue = OTSU(Imag)

Max = max(Imag(:));              % 最大值
Min = min(Imag(:));               % 最小值
T =Min:Max;                        % 灰度值范围
Tval = zeros(size(T));               % 方差


[Row,Col]=size(Imag);


imagSize = Row*Col;            % 像素点数量
% 遍历灰度值
for i = 1 : length(T)
    TK = T(i);
    iFg = 0;          % 前景
    iBg = 0;          % 背景
    FgSum = 0;    % 前景总数
    BgSum = 0;    % 背景总数
    W1=Imag>TK;
    WK1=W1.*Imag;
    iFg1=sum(W1(:));%求前景比例
    iBg1=imagSize-iFg1;
    FgSum1=sum(WK1(:));
    BgSum1=sum(Imag(:))-FgSum1;

    w0 = iFg1/imagSize;      % 前景比例
    w1 = iBg1/imagSize;     % 背景比例
    u0 = FgSum1/iFg1;         % 前景灰度平均值
    u1 = BgSum1/iBg1;        % 背景灰度平均值
    Tval(i) = w0*w1*(u0 - u1)*(u0 - u1);     % 计算方差
end
[~, flag] = max(Tval);             % 最大值下标
TValue = T(flag);
end