OTSU阈值分割算法

阈值分割算法

二值化

首先以灰度图的 $x,y$ 坐标为二维坐标系的 $x,y$ 坐标，以对应位置的像素灰度值为 $z$ 坐标，建立灰度图的三维坐标系,如下图

二值化处理就是在对应 $0-255$ 范围内找出合适阈值筛选出对我们有用的信息。如果把上图当作一个地形图，阈值分割就是找出合适的水位去淹没对我们无用的信息，保留有用信息。具体效果如下图

上述两图的灰度图像为

现在的难点是如何找到合适的阈值去分割图像信息，准确的说是如何用数学的方法使阈值的选取转化成计算机可以运行的过程。

OTSU阈值分割算法

对于上述灰度等级为 ${0,1...,255}$ 分辨率为 $120*188$ 的灰度图像

设 $n_i$ 为灰度等级为 $i$ 的像素个数，则图像像素总数为 $N=n_1+n_2+...+n_{255}=120*188$ ,

灰度为 $i$ 的像素概率为 $P(i)=\frac{n_i}{N}$

以 $i$ 为 $x$ 轴， $n_i$ 对应的像素个数为 $y$ 轴绘制像素直方图如下图

对于赛道图来说，由于赛道和背景有较大差异(灰度值不同)，且赛道在图像内占比较大，数以直方图内一般会出现两个波峰，分别代表赛道灰度像素与背景灰度像素。

现在我们假设存在一个阈值 $k$ ,对灰度值大于 $k$ 的部分称为前景，小于 $k$ 的部分称为背景

像素分类到前景的概率为 $P_1(k)=\frac{1}{N}\sum_{i=0}^{k}n_i$

分类到背景的概率为 $P_2(k)=1-P_1(k)$

背景的像素灰度均值为：

m e a n_{1} (k) = \frac{1}{n_{1} + . . . . + n_{k}} \sum_{i = 0}^{k} i n_{i} = \frac{N}{n_{1} + . . . . + n_{k}} \frac{1}{N} \sum_{i = 0}^{k} i n_{i} = \frac{N}{n_{1} + . . . . + n_{k}} \sum_{i = 0}^{k} i \frac{n_{i}}{N} = \frac{1}{P_{1} (k)} \sum_{i = 0}^{k} i P (i)

$mean_1(k)=\frac{1}{n_1+....+n_k}\sum_{i=0}^{k}in_i \\ =\frac{N}{n_1+....+n_k}\frac{1}{N}\sum_{i=0}^{k}in_i\\ =\frac{N}{n_1+....+n_k}\sum_{i=0}^{k}i\frac{n_i}{N}\\ =\frac{1}{P_1(k)}\sum_{i=0}^{k}iP(i)$

同理可得前景像素灰度均值为：

m e a n_{2} = \frac{1}{P_{2} (k)} \sum_{i = k + 1}^{255} i P (i)

$mean_2=\frac{1}{P_2(k)}\sum_{i=k+1}^{255}iP(i)$

整张图像的灰度均值为：

m e a n_{G} = \frac{1}{N} \sum_{i = 0}^{255} i n_{i} = \sum_{i = 0}^{255} i P (i) = P_{1} (k) m e a n_{1} (k) + P_{2} (k) m e a n_{2} (k)

$mean_G=\frac{1}{N}\sum_{i=0}^{255}in_i\\ =\sum_{i=0}^{255}iP(i) =P_1(k)mean_1(k)+P_2(k)mean_2(k)$

类间方差定义为阈值两侧数据各自均值距离总均值的加权方差，阈值两侧均值分别代表了前景与背景的主要灰度值，类间方差越大表示前景与背景可以得到更好的分离，二值化效果也更佳。

σ (k) = P_{1} (k) [m e a n_{1} (k) - m e a n_{G}]^{2} + P_{2} (k) [m e a n_{2} (k) - m e a n_{G}]^{2} = P_{1} (k) [m e a n_{1} (k) - P_{1} (k) m e a n_{1} (k) - P_{2} (k) m e a n_{2} (k)]^{2} + P_{2} (k) [m e a n_{2} (k) - P_{1} (k) m e a n_{1} (k) - P_{2} (k) m e a n_{2} (k)]^{2} = P_{1} (k) [P_{2} (k) m e a n_{1} (k) - P_{2} (k) m e a n_{2} (k)]^{2} + P_{2} (k) [P_{1} (k) m e a n_{2} (k) - P_{1} (k) m e a n_{1} (k)]^{2} = P_{1} (k) P_{2} (k)^{2} [m e a n_{1} (k) - m e a n_{2} (k)]^{2} + P_{2} (k) P_{1} (k)^{2} [m e a n_{2} (k) - m e a n_{1} (k)]^{2} = P_{1} (k) P_{2} (k) [[P_{2} (k) + P_{1} (k)] [m e a n_{1} (k) - m e a n_{2} (k)]^{2}] = P_{1} (k) P_{2} (k) [m e a n_{1} (k) - m e a n_{2} (k)]^{2}

$\sigma(k)=P_1(k)[mean_1(k)-mean_G]^2+P_2(k)[mean_2(k)-mean_G]^2\\ =P_1(k)[mean_1(k)-P_1(k)mean_1(k)-P_2(k)mean_2(k)]^2\\ +P_2(k)[mean_2(k)-P_1(k)mean_1(k)-P_2(k)mean_2(k)]^2\\ =P_1(k)[P_2(k)mean_1(k)-P_2(k)mean_2(k)]^2\\ +P_2(k)[P_1(k)mean_2(k)-P_1(k)mean_1(k)]^2\\ =P_1(k)P_2(k)^2[mean_1(k)-mean_2(k)]^2\\ +P_2(k)P_1(k)^2[mean_2(k)-mean_1(k)]^2\\ =P_1(k)P_2(k)[[P_2(k)+P_1(k)][mean_1(k)-mean_2(k)]^2]\\ =P_1(k)P_2(k)[mean_1(k)-mean_2(k)]^2$

最后得出类区间方差公式为：

σ (k) = P_{1} (k) P_{2} (k) [m e a n_{1} (k) - m e a n_{2} (k)]^{2}

$\sigma(k)=P_1(k)P_2(k)[mean_1(k)-mean_2(k)]^2$

遍历灰度的直方图数组，找出使 $\sigma(k)$ 最大的 $k$ 值即为阈值。

代码实现：

function Threshold = OTSU(pic ,W ,H)
HistGram = zeros(256,1);            		%直方图一维矩阵
pixelsum = uint32(0);               		   %像素总个数
PixelIntegral = uint32(0);                             %灰度总个数
PixelFore = uint32(0);             	                   %前景像素个数指灰度值>阈值
PixelBack = uint32(0);           	                 %背景像素个数指灰度值<阈值
P_Fore = double(0);                                      %前景概率
P_Back = double(0);                                     %背景概率
PixelIntegralFore = uint32(0);     	               %前景灰度总值
PixelIntegralBack = uint32(0);                      %背景灰度总值
MIntegralFore = double(0);          		%前景灰度均值
MIntegralBack = double(0);          		%背景灰度均值
tau = double(0);                    			%类间方差
last_tau = double(-1);
for i = 1:W
    for j = 1:H
         HistGram(pic(i,j)) = HistGram(pic(i,j)) + 1;       %遍历图像完善直方图内容
    end
end
for i = 1:256
    pixelsum = pixelsum + HistGram(i);			%获取像素总数可简化为W*H
    PixelIntegral = PixelIntegral + HistGram(i)*i;	   %灰度值总个数
end
for i = 1:256
    PixelFore = PixelFore + HistGram(i);		   %获取前景灰度总值
    PixelBack = pixelsum - PixelFore; 			  %背景灰度值
    P_Fore = double(PixelFore)/double(pixelsum);   %前景概率
    P_Back = 1 - P_Fore;					   %背景概率
    PixelIntegralFore = PixelIntegralFore + HistGram(i)*i;
    PixelIntegralBack = PixelIntegral - PixelIntegralFore;
    MIntegralFore = double(PixelIntegralFore)/double(PixelFore);
    MIntegralBack = double(PixelIntegralBack)/double(PixelBack);
    tau = double(P_Fore*P_Back*(MIntegralFore-MIntegralBack)^2);	%类区间方差公式
    if tau > last_tau											%取阈值为使tau为最大时
        last_tau = tau;
        Threshold = i;
    end     
end
end