一、实验目的
1、理解图像分割的基本概念。
2、掌握阈值法、K-means聚类方法、边缘提取及区域生长和分裂方法进行图像分割。
二、实验环境
Matlab 2020B
三、实验内容
题目
1、分别采用两种阈值选取方法实现图像分割（如全局阈值、OTSU等），要求根据阈值选取的思想自己写代码。（分割图像可自由选择）
2、采用K-means聚类算法实现图像分割，要求根据K-means的思想自己写代码。（分割图像可自由选择）
3、分别用Roberts,Sobel和拉普拉斯高斯算子对图像进行边缘检测(可使用系统函数)，比较三种算子处理的不同之处。
4、选择适当方法实现肺的分割，结果包括两部分：肺（白色显示）和背景（黑色显示）。
相关知识
部分核心代码
I=imread("couplenew.png"); I=im2gray(I); I=im2double(I); I=im2bw(I,127/255); imshow(I)
随意地选取了一个灰度级作为阈值，对图像进行分割。
I=imread("lab41.jpg"); I=im2gray(I); I=im2double(I); [h,w]=size(I); eps=1e-6; T=graythresh(I); I=im2bw(I,T); imshow(I)
调用函数graythresh求解阈值，对图像进行分割。graythresh用于计算图像的Otsu阈值。
I=imread("couplenew.png"); I=im2gray(I); I=im2double(I); [h,w]=size(I); eps=1e-6; T1=0.5; T2=1; while(abs(T2-T1)>=eps) T2=T1; G1=0; G1cnt=0; G2=0; G2cnt=0; for i=1:h for j=1:w if I(i,j)>T1 G2=G2+I(i,j); G2cnt=G2cnt+1; else G1=G1+I(i,j); G1cnt=G1cnt+1; end end end m1=G1/G1cnt; m2=G2/G2cnt; T1=1/2*(m1+m2); end I=im2bw(I,T1); imshow(I);
利用基本的全局阈值处理方法，计算图像的阈值。
I=imread("couplenew.png"); I=im2gray(I); [h,w]=size(I); %hist=histogram(I); nq=zeros(1,256); for i=1:h for j=1:w nq(I(i,j)+1)=nq(I(i,j)+1)+1; end end p=zeros(1,256); for i=1:256 p(i)=nq(i)/(h*w); end maxpos=1; maxsigma=0; P=zeros(1,256); P(1)=p(1); for i=2:256 P(i)=P(i-1)+p(i); end mg=sum(sum(I))/(h*w); for k=1:256 mk=0; for i=1:k mk=mk+(i-1)*p(i); end if (P(k)==0||P(k)==1) sigma=0; else sigma=(mg*P(k)-mk).^2/(P(k)*(1-P(k))); end if sigma>maxsigma maxsigma=sigma; maxpos=k; end end I=im2double(I); I=im2bw(I,maxpos/255); imshow(I)
利用Otsu算法计算图像的阈值。算法的思路、过程与前文是完全一致的。
%K-Means算法实现图像分割 k=2; I = imread('lab41.jpg'); figure; subplot(2,2,1); imshow(I),title('原始图像'); gray=rgb2gray(I); [m,n]=size(gray); I = reshape(I(:, :, 1), m*n, 1); res = kmeans(double(I), k); result = reshape(res, m, n); subplot(2,2,2); imshow(label2rgb(result)),title(strcat('K=',num2str(k),'时RGB通道分割结果')); res = kmeans(double(I), k+1); result = reshape(res, m, n); subplot(2,2,3); imshow(label2rgb(result)),title(strcat('K=',num2str(k+1),'时RGB通道分割结果')); res = kmeans(double(I), k+2); result = reshape(res, m, n); subplot(2,2,4); imshow(label2rgb(result)),title(strcat('K=',num2str(k+2),'时RGB通道分割结果'));
调用系统自带的K-means算法，实现图像分割。
I=imread("lab41.jpg"); I=im2gray(I); k=3; Mu=[1,1,1;127,127,127;255,255,255]; eps=1e-6; normMuDiff=[]; m=size(I,1); n=size(I,2); I=double(I); M=reshape(I,m*n,size(I,3)); flag=1; Cov=repmat(eye(size(I,3)),1,1,k); %k个聚类的协方差 Label=zeros(m*n,1); iter=1; figure; for iter=1:5 old_Mu=Mu; Num=zeros(1,k); ClassedPixs=zeros(m*n,size(I,3),k); %分到k个类中的像素值 w=zeros(m*n,k); %各像素点对于聚类的权重矩阵 C=zeros(size(I,3),size(I,3),k); for j=1:k invCov=inv(Cov(:,:,j)); C(:,:,j)=invCov/norm(invCov);%使用Mahalanobis距离(即马氏距离)度量 end for i=1:m*n dis=zeros(k,1); for j=1:k dis(j)=(M(i,:)-Mu(j,:))*C(:,:,j)*(M(i,:)-Mu(j,:))'; end maxIdx=find(dis==min(dis)); label=maxIdx(1); Label(i)=label; Num(label)=Num(label)+1; ClassedPixs(Num(label),:,label)=M(i,:); w(Num(label),label)=1/(sqrt(min(dis))+1); %避免被0除的情况 end w=w./repmat(sum(w),m*n,1); %权重归一化 % 更新协方差Cov for j=1:k PixsThisClass=ClassedPixs(1:Num(j),:,j); if Num(j)~=0 Cov(:,:,j)=cov(PixsThisClass)+0.01*rand(size(I,3)); %加随机量以防止矩阵奇异 Mu(j,:)=sum(repmat(w(1:Num(j),j),1,size(I,3)).*PixsThisClass); end end normDiff=norm(old_Mu-Mu); if normDiff<eps flag=0; end normMuDiff=[normMuDiff;normDiff]; sg=reshape(Label,m,n); end imshow(mat2gray(sg))
实现了K-means算法，将其展开实现。先随机选取k kk个聚类中心，度量每个像素点到每个聚类中心的“距离”并按“距离”分类。对所有像素点的分类完成后，根据每个像素点的类型计算更新其“聚类中心”的灰度级值。
I=imread("couplenew.png"); I=rgb2gray(I); I=im2double(I); I=imbinarize(I); [h,w]=size(I); newimg=zeros(h,w); for i=1:h-1 for j=1:w-1 temp=(I(i,j)-I(i+1,j+1)).^2+(I(i+1,j)-I(i,j+1)).^2; newimg(i,j)=sqrt(temp); end end subplot(1,2,1);imshow(I); subplot(1,2,2);imshow(newimg);
实现了利用Roberts算子对图像进行边缘检测。
I=imread('couplenew.png'); I=im2gray(I); I=im2double(I); G=I; [h,w]=size(I); res=zeros(h,w); for i=2:h-1 for j=2:w-1 gx=-G(i-1,j-1)-2*G(i-1,j)-G(i-1,j+1)+G(i+1,j-1)+2*G(i+1,j)+G(i+1,j+1); gy=-G(i-1,j-1)-2*G(i,j-1)-G(i+1,j-1)+G(i-1,j+1)+2*G(i,j+1)+G(i+1,j+1); res(i,j)=sqrt(gx.^2+gy.^2); end end res=res-G; subplot(1,2,1);imshow(G); subplot(1,2,2);imshow(res);
实现了利用Sobel算子对图像进行边缘检测。
I = imread('couplenew.png'); I=im2gray(I); I=im2double(I); I1=edge(I,"Roberts"); subplot(2,2,1);imshow(I);title("原图") subplot(2,2,2);imshow(I1);title("Roberts") I2=edge(I,"Sobel"); subplot(2,2,3);imshow(I2);title("Sobel") I3=edge(I,"log"); subplot(2,2,4);imshow(I3);title("log")
调用系统函数实现了Roberts算子、Sobel算子、log算子的图像边缘检测。
I=imread("lung.png"); I=im2gray(I); I=im2double(I); subplot(2,2,1);imshow(I);title("原图") [h,w]=size(I); eps=1e-6; T1=0.5; T2=1; while(abs(T2-T1)>=eps) T2=T1; G1=0; G1cnt=0; G2=0; G2cnt=0; for i=1:h for j=1:w if I(i,j)>T1 G2=G2+I(i,j); G2cnt=G2cnt+1; else G1=G1+I(i,j); G1cnt=G1cnt+1; end end end m1=G1/G1cnt; m2=G2/G2cnt; T1=1/2*(m1+m2); end I=im2bw(I,T1); I=imcomplement(I); subplot(2,2,2);imshow(I);title("二值化") L=bwlabel(I,8); %连通标记 s = regionprops(I,'Area');%计算各个连通区域面积 I=ismember(L,find([s.Area]<=18000&[s.Area]>=400)); subplot(2,2,3);imshow(I);title("去背景") se = strel('disk',7); I=imclose(I,se); subplot(2,2,4);imshow(I);title("闭操作")
实现了对肺部的分割，首先利用迭代得到分割的阈值，之后将黑白部分翻转。找出图像中所有8连通域，计算各个连通区域的面积，忽略面积过大的部分（背景）和面积过小的部分，从而正确地识别出两个肺。为了满足题目要求，把肺全部白色显示，还需要进行形态学操作的闭运算，过滤掉肺中的噪点与纹路。
实验结果
随意选取一个数值作为阈值，对图像进行二值化。
调用graythresh，所得结果为阈值，对图像进行二值化。
使用全局阈值法，对图像进行二值化。
使用全局阈值法，对图像进行二值化。
使用系统函数实现kmeans进行图像分割。
自行编写函数，实现kmeans进行图像分割。k=3
自行编写函数实现Roberts算子的边缘检测。
自行编写函数实现Sobel算子的边缘检测。
自行编写函数实现各种算子的边缘检测。
实现肺的分割。首先对图像进行二值化，然后过滤掉背景，最后对肺进行形态学操作的闭运算。
实验结果分析
实验需要比较分析不同边缘检测算子的效果与不同之处。其中，Roberts算子不自带平滑效果，对噪声较为敏感，但对边缘定位准确。Sobel算子检测的图像边缘宽度可能大于两个像素。对灰度渐变的低噪声图像存在良好的检测效果，但对混合的多复杂噪声图像效果不明显。Log算子自带平滑效果，但在平滑效果较好时，图像细节损失也大，边缘精度也会降低。
