Ostu算法其实就是遍历每个灰度级,判断哪个灰度级的阈值分割的效果最佳,判断效果好坏的指标就是类间方差的大小,类将方差越大效果越好。
- 计算图像的归一化直方图
```javascript
//
grayimg = rgb2gray(img);%将图像转为灰度图
grayimg = double(grayimg);
%计算直方图,每个灰度级的概率
[r c] = size(grayimg);
histogram = zeros(1,256);
for i=1:r
for j=1:c
histogram(grayimg(i,j)+1) = histogram(grayimg(i,j)+1)+1;
end
end
%归一化
histogram = histogram./(r*c);
```
- 遍历每一个灰度级,根据灰度级将灰度图分为两个部分,这里记为类1和类2,根据直方图计算类1和类2的概率记为p1、p2
```javascript
%计算两个类的累计概率和,k代表当前灰度级,k+1是因为灰度级从0开始,而matlab数组下标从1开始
p1 = sum(histogram(1:k+1));%类1概率的累加概率
p2 = sum(histogram(k+2:256));%类2概率的累加
```
- 计算类1、类2和全局的平均灰度值(类似期望)--这里会用到贝叶斯公式P(A/B)=P(B/A)P(A)/P(B) 其实这个公式的意思就是在B的情况下A发生的概率。
```javascript
%计算两个类的均值和全局的均值
m1 = 0;
m2 = 0;
mG = 0;
for i=0:k
m1 = m1 + i*histogram(i+1)/p1; % histogram(i+1)/p1为i灰度级在类1的概率,histogram(i+1)是在全局的概率。
end
for i=k+1:255
m2 = m2 + i*histogram(i+1)/p2;%同上
end
for i=0:255
mG = mG + i*histogram(i+1);%同上
end
```
- 根据上面算出的平均灰度值计算类间方差和全局方差
```javascript
%全局方差
for i=0:255
Gg = Gg + (i-mG)^2*histogram(i+1);
end
%类间方差
Gb = p1*(m1-mG)^2+p2*(m2-mG);
```
- 利用全局方差和类间方差相除可以得到一个可分性度量,用来表示阈值分割效果的好坏,这里就不计算了,下面给出matlab的全部代码
```javascript
%otsu算法和进行阈值分割
function oimg = ostuFnc(img)
grayimg = rgb2gray(img);%将图像转为灰度图
grayimg = double(grayimg);
%计算直方图,每个灰度级的概率
[r c] = size(grayimg);
histogram = zeros(1,256);
for i=1:r
for j=1:c
histogram(grayimg(i,j)+1) = histogram(grayimg(i,j)+1)+1;
end
end
%归一化
histogram = histogram./(r*c);
%初始化一个阈值,随后遍历每一个灰度级K
k = 0;
result = zeros(1,256);
while (k<256)
%计算两个类的累计概率和
p1 = sum(histogram(1:k+1));%类1概率的累加概率
p2 = sum(histogram(k+2:256));%类2概率的累加
%计算两个类的均值和全局的均值
m1 = 0;
m2 = 0;
mG = 0;
for i=0:k
m1 = m1 + i*histogram(i+1)/p1;
end
for i=k+1:255
m2 = m2 + i*histogram(i+1)/p2;
end
for i=0:255
mG = mG + i*histogram(i+1);
end
%计算全局方差和类间方差
Gg = 0;
Gb = 0;
%全局方差
for i=0:255
Gg = Gg + (i-mG)^2*histogram(i+1);
end
%类间方差
Gb = p1*(m1-mG)^2+p2*(m2-mG);
%获得使类间方差最大的K值,可能有多个,多个则求平均
result(k+1) = Gb;
k = k+1;
end
[maxGb KK] = max(result);
oimg = grayimg;
%阈值分割
for i=1:r
for j=1:c
if(grayimg(i,j)>KK)
oimg(i,j) = 0;
else
oimg(i,j) = 1;
end
end
end
end
```
- 处理结果