demo:图像处理 matlab实现迭代法全局阈值分割和局部自适应阈值分割

全局阈值分割

迭代法阈值获取方法：
1）选取一个的初始估计值T（128)；
2）用T分割图像。这样便会生成两组像素集合：G1由所有灰度值大于T的像素组成，而G2由所有灰度值小于或等于T的像素组成。
3）对G1和G2中所有像素计算平均灰度值u1和u2。
4）计算新的阈值：T=1/2(u1 + u2)。
重复步骤（2）到（4），直到得到的T值之差小于一个事先定义的参数T0。

先导进来图片

I = imread('rice.png');
[width,height] = size(I);
T0 = 1;%设置门限
T1= 128; %初始阈值T1
%设置G1，G2两个列向量，各自统计<T和>T的值
gray_leval_1 = 1;
gray_leval_2 = 1;

matlab里没有do-while语句，所以用
while 1
…
break
代替

while 1
    for i = 1:width
        for j = 1:height
            if I(i,j)>T1
                G1(gray_leval_1) = I(i,j); %得到分组G1
                gray_leval_1 = gray_leval_1 + 1;
            else
                G2(gray_leval_2) = I(i,j); %得到分组G2
                gray_leval_2 = gray_leval_2 + 1;
            end
        end
    end
    %计算G1、G2均值
    avg1 = mean(G1);
    avg2 = mean(G2);
    T2 = (avg1 + avg2)/2;
    if abs(T2 - T1)<T0 
        break;
    end
    T1 = T2;
    gray_leval_1 = 1;
    gray_leval_2 = 1;
end
T1
T1=uint8(T1);

T1=130.7
然后画出来

for i=1:width
    for j=1:height
        if(I(i,j)<T1)
            BW1(i,j)=0;
        else 
            BW1(i,j)=1;
        end
    end
end
subplot(1,3,1);imshow(I)
subplot(1,3,2);imshow(BW1);title('阈值分割');

matlab自带的Otus阈值分割方法

T3=graythresh(I)
BW2=imbinarize(I,T3);%Otus阈值分割
subplot(1,3,3);imshow(BW2);title('Otus阈值分割');

I=130.9,差不多哈。
结果

这个阈值分割方法存在偏差，如果两个峰值数目相差越多，偏差越大
（假设第一次迭代前，T在正确位置，很容易看出来，随着迭代进行，越来越偏离）

计算的T=93，目测值应该在110-120吧。

多说一句，图像处理，一定要double一下，否则可能出现问题。比如，我想把大米图的直方图暗部拉大（偏向0），以体现这种偏差，就用了个一次函数的变换，结果这样了

因为先乘k，结果>255的全归于255，以后全错。

自适应阈值分割

1）将整幅图像分成一系列互相之间有50%重叠的子图像；
2）做出每个子图像的直方图；
3）检测各个子图像的直方图是否为双峰，如果是，则采用最佳阈值法确定一个阈值，否则就不进行处理；
4）根据对直方图为双峰的子图像得到的阈值通过插值得到所有子图像（像素点）的阈值。
仨坑
1、怎么割图像啊（这个还好吧，也不算割，设置好步长来就行了）
2、怎么检测单峰还是双峰啊，单峰的话，空缺imresize是不行的。
3、怎么插值啊（主要是老师上课讲的太玄学，按照距离加权，四点平均，把相距很远的点内部全部插满，听懵了。后来想到，matlab直接imresize啊！然后就码下来啦）
1、分割方式设置好矩形大小和步长就行
2、峰值个数没写，matlab统计工具箱里好像有正态拟合图像，而不是plot（）这种拟合，plot拟合只是画的圆滑了些还是很多局部峰值

3、插值方式采用了imresize函数
A=imresize(I,2,‘nearest’);%最近邻插值
B=imresize(I,2,‘bilinear’);%双线性
C=imresize(I,2,‘bicubic’);%双三次

I=imread('rice.png');
T=zeros(15);
for i=0:14
    for j=0:14
        img = imcrop(I,[1+i*16 1+j*16 31 31])；%内核32*32，步长设置为16
        t=255*graythresh(img);%我这里直接用了
        T(j+1,i+1)=t;%这里注意坐标顺序反着
    end
end
T=uint8(T);
T1=imresize(T,[256 256],'bilinear');%看了下貌似双线性插值最好
%%
for i=1:256
    for j=1:256
        if (I(i,j)<T1(i,j))
            BW(i,j)=0;
        else
            BW(i,j)=1;
        end
    end
end
imshow(BW)

有一个地方要注意，图像[i,j],矩阵得是[j,i]。

阈值图像：

底部更暗，所以阈值也更低才能滤完整
阈值分割结果

嗯，上部噪声少了，下部也更完整了
然后试了下这位兄弟博客的图片
这老哥跑的OPENCV中自带的函数，效果一级棒（怎么看着像边缘检测）
我这个…毕竟没字的地方有黑背景有白背景，一加权，也没法搞，不过还阔以吧

matlab中自带的阈值分割，对大米图还行。全局阈值分割底部大米都割残了，这个方法完整，可是，这跟局部有关系？貌似先调了调灰度，就好分割了

 I=imread('rice.png');
 I=im2double(I);
se=strel('disk',10);%se.Neighborhood:10*10圆盘
ft=imtophat(I,se);% 高帽变换
gt=uint8(255*ft);
Th=graythresh(ft);                       
G=imbinarize(ft,Th);%阈值分割                       
figure,imshow(G),title('局部阈值');

Otsu算法

matlab自带的阈值分割方法
参考动态阈值–大津法和基于Otsu算法的图像自适应阈值分割，重新打了一下
调用：

[counts x] = imhist(img);
[m n] = size(img);
BW1 = otsu(counts, m*n);
第一个元素为直方图，第二个元素为图片大小
或者：
T=graythresh(I)
BW2=imbinarize(I,T);

原理：

让T遍历所有灰度值，找到g的最大值
代码：

function level = otsu(histogramCounts, total)%histogramCounts是图像的直方图, total图像的总像素数
sum0 = 0;
w0 = 0;%像素个数
maximum = 0.0;%保存T的最大值
total_value = sum((0:255).*histogramCounts');
for ii=1:256
    w0 = w0 + histogramCounts(ii);%<T像素个数累加
    if (w0 == 0)
        continue;
    end
    w1 = total - w0;
    if (w1 == 0)%遍历结束>T像素=0
        break;
    end
    sum0 = sum0 +  (ii-1) * histogramCounts(ii);%灰度和
    m0 = sum0 / w0;%平局灰度u0
    m1 = (total_value - sum0) / w1;%平局灰度u1
    icv = w0 * w1 * (m0 - m1) * (m0 - m1);
    if ( icv >= maximum )
        level = ii;
        maximum = icv;
    end
end
end