目录
一、实验目的二、实验环境三、实验内容题目相关知识部分核心代码实验结果实验结果分析 四、实验小结一、实验目的
1、理解图像分割的基本概念。
2、掌握阈值法、K-means聚类方法、边缘提取及区域生长和分裂方法进行图像分割。
二、实验环境
Matlab B
三、实验内容
题目
1、分别采用两种阈值选取方法实现图像分割（如全局阈值、OTSU等），要求根据阈值选取的思想自己写代码。（分割图像可自由选择）
2、采用K-means聚类算法实现图像分割，要求根据K-means的思想自己写代码。（分割图像可自由选择）
3、分别用Roberts,Sobel和拉普拉斯高斯算子对图像进行边缘检测(可使用系统函数)，比较三种算子处理的不同之处。
4、选择适当方法实现肺的分割，结果包括两部分：肺（白色显示）和背景（黑色显示）。
相关知识
实现图像分割时，需要进行阈值选取。阈值选取可以使用全局阈值法、Otsu法这两种阈值选取法。
全局阈值法过程如下：1.为全局阈值选取一个初始估计值 T T T。2.使用 T T T分割图像，产生了两组像素：由所有灰度值大于 T T T的像素组成的 G 1 G_1 G1，由所有灰度值小于等于 T T T的像素组成的 G 2 G_2 G2。3.分别计算区域 G 1 G_1 G1和 G 2 G_2 G2中像素的平均灰度值 m 1 m_1 m1和 m 2 m_2 m2。4.计算一个新的阈值： T = 1 2 ( m 1 + m 2 ) T=\frac 12(m_1+m_2) T=21(m1+m2)。5.重复步骤2到4，直到后续迭代中 T T T的差小于一个预定义的 Δ T \Delta T ΔT为止。此时所得到 T T T的为图像分割所需的阈值。
Otsu方法可以进行最佳全局阈值处理。首先预处理出图像的直方图分量。 p q = n q n , q = 0 , 1 , 2 , . . . , L − 1 p_q=\frac{n_q}n,q=0,1,2,...,L-1 pq=nnq,q=0,1,2,...,L−1，其中 L L L为灰度级数， n n n为图像中像素点的总数， n q n_q nq为灰度级为 q q q的像素点的总数量。对 p q p_q pq进行前缀和运算，求得 P 1 ( k ) = ∑ i = 0 k p i P_1(k)=\sum_{i=0}^{k}p_i P1(k)=∑i=0kpi，则 P 2 ( k ) = ∑ i = k + 1 L − 1 p i = 1 − p 1 ( k ) P_2(k)=\sum_{i=k+1}^{L-1}p_i=1-p_1(k) P2(k)=∑i=k+1L−1pi=1−p1(k)。记整个图像的平均灰度 m G = ∑ i = 0 L − 1 i p i m_G=\sum_{i=0}^{L-1}ip_i mG=∑i=0L−1ipi。取 m ( k ) = ∑ i = 0 k i p i m(k)=\sum_{i=0}^kip_i m(k)=∑i=0kipi，则类间方差可写作 σ B 2 ( k ) = [ m G P 1 ( k ) − m ( k ) ] 2 P 1 ( k ) [ 1 − P 1 ( k ) ] \sigma^2_B(k)=\frac{[m_GP_1(k)-m(k)]^2}{P_1(k)[1-P_1(k)]} σB2(k)=P1(k)[1−P1(k)][mGP1(k)−m(k)]2：使得类间方差 σ B 2 ( k ) \sigma_B^2(k) σB2(k)取得最大值的 k k k即为图像分割的阈值。
采用K-means算法实现图像分割时，流程如下：
1.首先需要给定 k k k个初始聚类中心 c 1 , c 2 , . . . , c k c_1,c_2,...,c_k c1,c2,...,ck。
2.循环的迭代进行步骤2~步骤4，直到聚类中心不再变化或迭代次数达到预设的上限。
3.计算每个点 p p p到每个聚类中心的距离 d i d_i di，将像素点 p p p归入离聚类中心 d j d_j dj最近的点，标记为第 j j j类点。
4.在将所有对于聚类中心 c 1 , c 2 , . . . , c k c_1,c_2,...,c_k c1,c2,...,ck，计算新的聚类中心的值。记 c i ′ c_i' ci′为被标记为第 i i i类点的所有点的平均值，更新 c i ′ c_i' ci′。
5.得到最终的聚类中心 c 1 ′ , c 2 ′ , . . . , c k ′ c_1',c_2',...,c_k' c1′,c2′,...,ck′。
在本实验中，处理图像利用K-means算法时，一种可能的实例化方案如下：“聚类中心”指的就是rgb的值 c 1 , c 2 , . . . , c k c_1,c_2,...,c_k c1,c2,...,ck。“距离”指的就是当前像素点的rgb值与各个“聚类中心”的rgb值欧氏距离、曼哈顿距离或其他形式的距离。在算法运行完成得到的最终“聚类中心” c 1 ′ , c 2 ′ , . . . , c k ′ c_1',c_2',...,c_k' c1′,c2′,...,ck′后，所有像素点都以与“聚类中心”的“距离”被重新分类、重新赋予新的颜色或rgb值，从而完成了图像前景与背景的分离。
在利用Roberts算子进行边缘检测时，实质就是用此掩膜对矩阵每个点进行运算： G x = ( 1 0 0 − 1 ) G_x=\begin{pmatrix}1&0\\0&-1\end{pmatrix} Gx=(100−1)， G y = ( 0 − 1 1 0 ) G_y=\begin{pmatrix}0&-1\\1&0\end{pmatrix} Gy=(01−10)。对于输入像素 f ( i , j ) f(i,j) f(i,j)，输出像素 g ( i , j ) g(i,j) g(i,j)，则： g ( x , y ) = ( f ( x , y ) − f ( x + 1 , y + 1 ) ) 2 + ( f ( x + 1 , y ) − f ( x , y + 1 ) ) 2 g(x,y)=\sqrt {(f(x,y)-f(x+1,y+1))^2+(f(x+1,y)-f(x,y+1))^2} g(x,y)=(f(x,y)−f(x+1,y+1))2+(f(x+1,y)−f(x,y+1))2
在利用Sobel算子进行边缘检测时，使用到的掩膜如下：
( − 1 0 1 − 2 0 2 − 1 0 1 ) \begin{pmatrix}-1&0&1\\-2&0&2\\-1&0&1\end{pmatrix} −1−2−1000121 用于计算在 x x x轴方向上的梯度 ∂ P ∂ x \frac{\partial P}{\partial x} ∂x∂P， ( − 1 − 2 − 1 0 0 0 1 2 1 ) \begin{pmatrix}-1&-2&-1\\0&0&0\\1&2&1\end{pmatrix} −101−202−101 用于计算在 y y y轴方向上的梯度 ∂ P ∂ y \frac{\partial P}{\partial y} ∂y∂P。
部分核心代码
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:hfor j=1:wif I(i,j)>T1G2=G2+I(i,j);G2cnt=G2cnt+1;elseG1=G1+I(i,j);G1cnt=G1cnt+1;endendendm1=G1/G1cnt;m2=G2/G2cnt;T1=1/2*(m1+m2);endI=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:hfor j=1:wnq(I(i,j)+1)=nq(I(i,j)+1)+1;endendp=zeros(1,256);for i=1:256p(i)=nq(i)/(h*w);endmaxpos=1;maxsigma=0;P=zeros(1,256);P(1)=p(1);for i=2:256P(i)=P(i-1)+p(i);endmg=sum(sum(I))/(h*w);for k=1:256mk=0;for i=1:kmk=mk+(i-1)*p(i);endif (P(k)==0||P(k)==1)sigma=0;elsesigma=(mg*P(k)-mk).^2/(P(k)*(1-P(k)));endif sigma>maxsigmamaxsigma=sigma;maxpos=k;endendI=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:5old_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:kinvCov=inv(Cov(:,:,j));C(:,:,j)=invCov/norm(invCov);%使用Mahalanobis距离(即马氏距离)度量endfor i=1:m*ndis=zeros(k,1);for j=1:kdis(j)=(M(i,:)-Mu(j,:))*C(:,:,j)*(M(i,:)-Mu(j,:))'; endmaxIdx=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除的情况endw=w./repmat(sum(w),m*n,1); %权重归一化% 更新协方差Covfor j=1:kPixsThisClass=ClassedPixs(1:Num(j),:,j);if Num(j)~=0Cov(:,:,j)=cov(PixsThisClass)+0.01*rand(size(I,3)); %加随机量以防止矩阵奇异Mu(j,:)=sum(repmat(w(1:Num(j),j),1,size(I,3)).*PixsThisClass);endendnormDiff=norm(old_Mu-Mu);if normDiff<epsflag=0;endnormMuDiff=[normMuDiff;normDiff];sg=reshape(Label,m,n);endimshow(mat2gray(sg))
实现了K-means算法，将其展开实现。先随机选取 k k k个聚类中心，度量每个像素点到每个聚类中心的“距离”并按“距离”分类。对所有像素点的分类完成后，根据每个像素点的类型计算更新其“聚类中心”的灰度级值。
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-1for j=1:w-1temp=(I(i,j)-I(i+1,j+1)).^2+(I(i+1,j)-I(i,j+1)).^2;newimg(i,j)=sqrt(temp);endend 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-1for j=2:w-1gx=-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);endendres=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:hfor j=1:wif I(i,j)>T1G2=G2+I(i,j);G2cnt=G2cnt+1;elseG1=G1+I(i,j);G1cnt=G1cnt+1;endendendm1=G1/G1cnt;m2=G2/G2cnt;T1=1/2*(m1+m2);endI=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算子自带平滑效果，但在平滑效果较好时，图像细节损失也大，边缘精度也会降低。
