参考：<<数字图像处理——使用MATLAB分析与实现>> 清华大学出版社 蔡利梅主编
1.图像分割
2.阈值分割方法
（这些都是阈值的选择方法）
a.基于灰度直方图阈值分割
对于imhist函数，1-256，数组没有0下标，但是颜色范围是0-255，所以差一个1
function [is, peak ] = Bimodal( histgram)
count = 0;
for j= 2:255
if histgram(j-1)<histgram(j) && histgram(j+ 1)< histgram(j) %局部极大值
count = count + 1;
peak(count)= j; %记录峰所在的位置
if count> 2
is = 0;
return; %多个阈值分割，这里不return
end
end
end
if count == 2
is = 1;
else
is = 0;
end
end
Image = rgb2gray(imread( 'face2.JPG'));
figure,imshow( Image) ,title('原图');imhist( Image);
hist1= imhist(Image);
%对于imhist函数，1-256，数组没有0下标，但是颜色范围是0-255，所以差一个1
hist2 = hist1;
iter = 0; %多迭代次数,限制循环次数
while 1
[ is, peak] = Bimodal(hist1); %判断是否为双峰直方图，是则找到峰
if is == 0 %非双峰直方图进行平滑
hist2(1)= (hist1(1)* 2+ hist1(2))/3; %第一个点
for j= 2:255 %所有灰度
hist2(j) = ( hist1(j- 1 ) + hist1(j)+ hist1(j+ 1))/3;%相邻3个点求平均以平滑直方图
end
hist2(256)= (hist1(255) + hist1(256)* 2)/3; %最后一个点
hist1 = hist2;
iter = iter+1;
if iter > 1000
break;
end
else %是双峰直接退出
break;
end
end
%整体对应的灰度应该-1，但对于数组是对的
[ trough, pos ] = min(hist1( peak(1) : peak (2) ) );%min从1开始，找双峰间的波谷，多个阈值分割，for循环寻找多个波谷
thresh = pos + peak(1); %波谷对应的灰度，for循环寻找多个波谷灰度
figure, stem( 1:256, hist1, 'Marker' , 'none ' );
hold on
stem( [ thresh,thresh], [ 0, trough], 'Linewidth',2); %做标记
hold off
result = zeros( size( Image) );
%for循环设置不同颜色
result( Image > thresh) = 1; %阈值化
figure,imshow( result) ,title( '基于双峰直方图的阈值化');
缺点：对于边缘恰好为最值未判断
b.基于模式分类思路的阈值选择
b.1 OTSU（最大类间方差法）
由两类推导值N类
Image = rgb2gray( imread( ' lotus.jpg'));
figure,imshow( Image) ,title('原始图像');
T= graythresh( Image) ; %获取阈值
result = im2bw ( Image, T);%二值化图像
% result = imbinarize( Image) ; %或者用本句代替上两行实现图像二值化，使用OTSU方法
figure,imshow( result) ,title( ' 0TSU方法二值化图像');
b.2 最大熵法
clear,clc
Image = rgb2gray( imread( 'face2.jpg '));
figure,imshow( Image),title('原图');
hist = imhist( Image); %
bottom = min(Image(: )) + 1;%图中最小灰度，+1防止下标为0,同时保证了和imhist下标的一致性
top = max( Image( :))+1;%图中最大灰度
J= zeros( 256,1);
for t = bottom + 1:top - 1 %阈值在bottom和top之间，其他地方没有值
po = sum( hist (bottom:t)); %当前阈值t下，目标区域概率
pb = sum( hist(t + 1:top) ) ;%当前阈值t下，背景区域概率
ho = 0;hb = 0;
for j = bottom:t
ho = ho- log(hist(j)/po +0.01) * hist(j)/po; %当前阈值t下，目标区域熵计算
end
for j=t+1:top
hb = hb - log( hist(j)/pb +0.01) * hist(j)/pb;%当前阈值t下,背景区域熵计算
end
J(t) = ho + hb;%当前阈值t下,熵函数计算
end
[ maxJ , pos ] = max(J( :));%熵函数最大值
result= zeros( size( Image));
result( Image > pos) = 1;%阈值化
figure,imshow( result);
b.3 最小误差法
Image = rgb2gray( imread( 'p3.jpg'));
figure,imshow (Image) , title('原图');hist = imhist( Image);
bottom = min(Image( : ))+1; top = max(Image( : ))+ 1; %灰度值范围
J= zeros(256,1);J=J+10000; %最小误差法误判概率
alpha = 0.25; %目标在图像中所占的比例
scope = find(hist>5);
%估计概率密度时,每一类要保证一定的样本数目,排除直方图两端
%概率很小的的灰度级,避免估计不准确导致计算偏差
minthresh = scope(1); maxthresh= scope(end);
if maxthresh>= top
maxthresh= top-1;
end
for t= minthresh+ 1 : maxthresh %选择最优的t
miuo = 0;sigmaho = 0;
for j= bottom : t
miuo = miuo + hist(j) * double(j);
end
pixelnum = sum(hist(bottom : t)) ;
miuo = miuo/pixelnum; %当前阈值下,求目标区域均值,加权平均
for j = bottom:t
sigmaho = sigmaho + ( double(j) -miuo)^2 * hist(j);
end
sigmaho = sigmaho/pixelnum; %当前阈值下,求目标区域方差
miub= 0;sigmahb = 0;
for j= t+1 :top
miub = miub + hist(j) * double(j);
end
pixelnum = sum( hist(t+1:top));
miub = miub/pixelnum; %当前阈值下,求背景区域均值
for j=t+ 1:top
sigmahb = sigmahb +( double(j) - miub)^2* hist(j);
end
sigmahb = sigmahb/pixelnum; %当前阈值下,求背景区域方差
Epsilonb = 0; Epsilono = 0; %各区域误判概率初始化
for j = bottom : t %当前阈值下,背景区域误判概率
pb = exp( - (double(j) - miub)^ 2/( sigmahb * 2 + eps) )/( sqrt(2 * pi * sigmahb) + eps) ;
Epsilonb = Epsilonb + pb;
end
for j=t+ 1:top %当前阈值下，目标区域误判概率
po = exp( - ( double(j) - miuo)^2/( sigmaho * 2+ eps))/( sqrt(2* pi *sigmaho) + eps);
Epsilono = Epsilono + po;
end
J(t)= alpha * Epsilono +(1- alpha) * Epsilonb; %当前阈值下，整体误判概率
end
[minJ , pos] = min(J(:));
%求最小误判概率
result = zeros(size(Image));
result( Image > pos)= 1;
figure,imshow( result);
3.边界分割
包括三个步骤：边界检测，边界改良，边界追踪
a.基于梯度的边界闭合
b.Hough变换（检测直线）
clear,clc;
Image = rgb2gray( imread( '123.png') );
bw = edge(Image, 'canny'); % canny算子边缘检测得二值边缘图像
figure,imshow(bw);
[h, theta,rho] = hough( bw, 'RhoResolution' ,0.5, 'ThetaResolution' ,0.5); % Hough变换,rho间隔0.5，Theta间隔0.5
%h矩阵，SHT(r,c) 中的 Q 值表示 xy 平面中的 Q 个点位于 theta(c) 和 rho(r) 指定的线条上
figure,imshow(imadjust(mat2gray(h)),'XData' ,theta, 'YData' ,rho,'InitialMagnification' , 'fit');
%显示Hough变换矩阵
xlabel( '\theta') , ylabel( '\rho' );
axis on, axis normal , hold on;
P= houghpeaks(h, 100); %出现最多的theta,rho对于这些点的最可能的图形,该函数返回p，即一个保留峰值的行和列坐标的矩阵。
x=theta(P(:,2));y= rho(P(:,1));
plot( x, y, 's' ,'color', 'r ');%获取并标出参数平面上的峰值点
lines = houghlines(bw,theta,rho,P, 'FillGap',5, 'Minlength' , 7); %检测图像中的直线段
figure,imshow(Image);
hold on;
max_len = 0;
for i= 1 : length(lines)
xy = [lines(i).point1 ; lines(i).point2];
plot(xy( : , 1), xy( :,2) , 'LineWidth' ,2, 'Color' , 'g'); %用绿色线段标注直线段
plot(xy(1,1),xy(1,2), 'x ' , 'Linewidth',2, 'Color' , 'y');
plot(xy(2,1 ) , xy(2,2) , 'x ', 'LineWidth', 2 , 'Color', 'r');%标注直线段端点
end
c.边界追踪
5.区域分割
a.区域生长
包括种子点选取，相似性准则，区域停止生长条件
Image = im2double( imread( 'face2.jpg' ));
[ height,width, channel] = size( Image) ;
if channel == 3
Image = rgb2gray( Image);
end
figure,imshow( Image);
[ seedx , seedy , button] =ginput(1);%交互式获取一个种子点
seedx = round( seedx) ; seedy = round( seedy) ; %取整
region = zeros( height , width); % 生长区域
region(seedy , seedx) = 1;
region_mean = Image( seedy , seedx); %初始化，一开始只有一个值
region_num = 1; %初始区域只有一个种子点
flag = zeros( height, width);
flag( seedy , seedx)= 1; % 处理过的点做标记,避免重复处理
neighbor =[-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];%8邻点
for k= 1:8 %8邻点
y = seedy + neighbor (k, 1); x = seedx + neighbor ( k,2);
waiting(k,: )= [ y, x];
flag(y, x)= 2; %待处理像素点，即当前点的八邻点，待处理的标为2
end
pos = 1; len = length(waiting);
while pos <len %是否存在待处理像素点
len = length( waiting); %不断更新
current = waiting(pos, : );pos = pos +1;
pixel = Image( current(1) , current(2)); %当前要判断的像素点
pdist = abs(pixel- region_mean) ; %当前像素点与区域灰度均值的距离
if pdist <40/255
region( current( 1 ), current(2))=1; %生长出来的像素点
region_mean = region_mean * region_num +pixel; %加入
region_num = region_num + 1;
region_mean = region_mean/region_num; %新区域求灰度均值
for k = 1:8 %对于新加入的点，将其8邻点加入 ,生长的一个具体体现
newpoint = current + neighbor(k, : );
if newpoint(1)> 0 && newpoint(1)<= height && newpoint(2)> 0 &&...
newpoint(2)< width && flag( newpoint(1 ) , newpoint(2)) == 0
waiting( end+1,: ) = newpoint;
flag(newpoint(1 ) , newpoint(2))=2; %新生长出来的点作为种子点,将其邻点备选
end
end
end
end
figure,imshow( region) ,title( '区域生长');
种子点选取：随机选取
相似性准则：新像元与全区域灰度均值之差小于40
停止条件：没有待处理点
思考：多种子点
b.区域合并
clear,clc;
Image = [1 0 4 6 5 1;1 0 4 6 6 2;0 1 5 5 5 1;0 0 5 6 5 0;0 0 1 6 0 1;1 0 1 2 1 1];
% Image = im2gray(imread("face2.jpg"));
[ height, width] = size( Image);
flag = zeros( height , width);%区域标记
thresh = 2; %合并阈值
neighbor =[-1 -1;-1 0;0 -1] ;%左上、上、左邻点，为什么取这个，因为后面的没有进行分类
flag(1,1)= 1; %左上角第一个点设为区域1
number = 1;
%前面的已经进行了归类，看后面的属于前面的哪一类
for j= 1 : height
for i= 1 : width
pdist =[ 300 300 300]; %当前点与左上、上、左邻点距离初始化,像素差不可能大于300
for k=1:3 %当前点与左上、上、左邻点距离
y=j+ neighbor(k,1); x= i+ neighbor(k,2);
if x>=1 && y>=1
pdist(k)= abs(Image(j,i)- Image(y,x));
end
end
[mindist, pos] = min(pdist( : ));
if mindist<= thresh %三个距离中最小的若符合合并规则,将对应邻点的标记赋予当前点
y= j+neighbor( pos,1);
x = i+ neighbor( pos,2);
if flag(y,x)
flag(j,i)= flag(y,x);
end
elseif mindist ~= 300 %若没有相似的点,则赋予当前点新的标记
number = number +1;
flag(j,i) = number;
end
end
end
for j= 2:height
for i= 2:width
if flag(j-1,i) ~= flag(j,i-1)&& ...
( flag(j,i)== flag(j-1,i)||flag(j,i) == flag(j,i-1))
%若上,左邻点标记不一致,但当前点和其中一个邻点标记一致,则判断两个区域是否是同一个。
%若是,则将两个区域标记修改为较小的一个，即区域合并
pdist = abs( Image(j- 1,i) -Image(j,i-1));
if pdist<= thresh
minv = min(flag(j- 1,i) ,flag(j,i-1) );
maxv = max( flag(j-1,i) ,flag(j,i-1) );
flag(flag == maxv) = minv;
end
end
end
end
前面的进行分类，后面的点考虑，左，上，左上点，若该点与其中一个点相似，则合并，否则自己成为一个新的类，对于合并一遍的数据，考虑左和上的数据，若相似，则将大的flag改为小的flag，即合并这两种类型
c.区域分裂
6.基于聚类的图像分割
k-means
1206
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
