迭代法阈值获取方法：
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('局部阈值');
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