分水岭图像分割和区域合并的matlab实现（毕业设计）

%% 基于形态学分水岭和区域合并的图像分割Matlab实现 % Date:2017/2/25 % 1.读取待处理图像 Img = imread('horse.jpg'); Img_gray = rgb2gray(Img); [width,height] = size(Img_gray); figure,imshow(Img),title('原始图像') figure,imshow(Img_gray),title('灰度图像') % 2.利用Sobel求取梯度图像 hy=fspecial('sobel');%生成sobel垂直梯度模板 hx=hy'; %生成sobel水平梯度模板 gradx=imfilter(double(Img_gray),hx,'replicate'); figure,imshow(abs(gradx),[]),title('图像的sobel水平梯度'); grady=imfilter(double(Img_gray),hy,'replicate'); figure, imshow(abs(grady),[]),title('图像的sobel垂直梯度'); grad=sqrt(gradx.^2+grady.^2); %得到图像的sobel梯度 figure,imshow(grad,[]),title('图像的sobel梯度图像'); L1 = watershed(grad); Lrgb1 = label2rgb(L1); figure,imshow(Lrgb1),title('梯度幅值直接做分水岭变换(传统的分水岭变换)') % 3.对灰度图像进行形态学开闭重建 %基于重建（reconstruct）的开—闭运算相对基本的开—闭运算来说，在不影响目标图象的前提下，去除图象中的碎片效果更好 %对灰度图象进行“opening-by-reconstruction”运算 %I = imcomplement(Img_gray); %对图像进行求反运算 se =strel('disk',20); %创建由指定形状shape对应的结构元素 Ie = imerode(Img_gray, se); %腐蚀图像 Iobr = imreconstruct(Ie,Img_gray); %数学形态学重建图像 figure,imshow(Iobr);title('基于重建的形态学开运算'); %对进行开运算后的图象进行“closing-by-reconstruction”运算 %对重建前后的图象都要进行求补运算 Iobrd = imdilate(Iobr,se); %膨胀图像 Iobrcbr = imreconstruct(imcomplement(Iobrd),imcomplement(Iobr)); Iobrcbr = imcomplement(Iobrcbr); figure,imshow(Iobrcbr);title('基于重建的形态学闭运算'); % 4.计算Iobrcbr的局部最大值二值图像，以得到更好的前景目标 fgm = imregionalmax(Iobr); %确定所有局部最大值，fgm是一个二值图像 figure,imshow(fgm);title('局部最大值二值图像'); g2= imcomplement(Iobrcbr); hs=round(255*double(g2)); g22=uint8(hs); %将前景目标叠加到原图上 I2 = Img; I21 = I2(:,:,1); I22 = I2(:,:,2); I23 = I2(:,:,3); I21(fgm) = 255; %I2与fgm重叠部分区fgm，不重叠的部分取255 I22(fgm) = 0;I23(fgm) = 0; I2 = cat(3,I21,I22,I23); figure,imshow(I2);title('前景目标叠加到原图上'); % 5.进行开闭运算平滑边缘，去掉图像中局部极小值区域，得到前景标记图 %清除标记边缘的模糊点 se2 = strel(ones(3,3)); fgm2 = imclose(fgm, se2); fgm3 = imerode(fgm2, se2); %去除所有小于特定像素数（此处设定为20）的斑点污渍 fgm4 = bwareaopen(fgm3, 20); %用于从对象中移除小对象，此处指灰度值小于20的对象,fgm4是二值图 figure,imshow(fgm4),title('前景标记图') I3 = Img; I31 = I3(:,:,1); I32 = I3(:,:,2); I33 = I3(:,:,3); I31(fgm4) = 255; %I3与fgm4重叠部分区fgm4，不重叠的部分取255 I32(fgm4) = 0;I33(fgm4) = 0; I3 = cat(3,I31,I32,I33); figure,imshow(I3);title('修改局部最大值叠加在原图上'); % 6.计算背景标记(外部约束) %将基于重建的开闭运算图象转化为二值图象（二值化） bw = im2bw(Iobrcbr,graythresh(Iobrcbr)); %grayshresh计算一个合适的阈值（OTSU阈值） figure,imshow(bw);title('基于重建的阈值化二值图'); % 7.计算距离函数，产生外部约束 D = bwdist(bw);%计算图像的欧几里得矩阵 %将不属于目标的像素置为—inf D = -D; D(bw) = -Inf; %D与bw重叠部分区bw，不重叠的部分取-Inf DL = watershed(D); bgm = DL == 0; %DL不变,bgm中非0元素变成0，0变成1,此处用来求取分割边界 figure,imshow(bgm);title('分水岭脊线（背景标记）'); % 8.采用强制最小技术重构梯度图 gradmag2 = imimposemin(grad, bgm | fgm4); figure,imshow(gradmag2),title('重构梯度图') % 9.分水岭算法图像分割 L = watershed(gradmag2); M = label2rgb(L); figure,imshow(L,[]),title('分水岭分割') figure,imshow(M,[]);title('分水岭分割的彩色显示'); % 10.区域合并 wat = L; stats=regionprops(wat,'Area'); %获取各区域的面积大小 area=[stats.Area]; [u,v]=sort(area); %对各个区域的面积从大到小进行排序 r=length(area); %r为区域的个数 [mm,nn]=size(Img); d=max(max(wat)); %d表示分水岭分割区域的个数 for i=1:d [a,b]=find(wat==i); %找到指定索引区域的像素点的横纵坐标 c=mean(g22(mm*(b-1)+a)); %求出指定索引区域像素均值 ag(i)=c; end I1=wat; for i=1:r a=v(i); [x,y]=find(I1==a); %找出指定区域的像素点 %下面根据区域像素相似性进行区域像素合并 w1=[x y]; w1(w1(:,1)==1|w1(:,1)==2|w1(:,1)==mm-1|w1(:,1)==mm|w1(:,2)==1|w1(:,2)==2|w1(:,2)==nn-1|w1(:,2)==nn,:)=[]; if isempty(w1) I1(I1==a)=v(r); elseif u(i)<5000 w2=w1; a1=w1(:,1); b1=w1(:,2); w1(find(I1(mm*(b1-2)+a1-1)==0|I1(mm*(b1-2)+a1)==0|I1(mm*(b1-2)+a1+1)==0|I1(mm*(b1-1)+a1-1)==0|I1(mm*(b1-1)+a1+1)==0|I1(mm*b1+a1-1)==0|I1(mm*b1+a1)==0|I1(mm*b1+a1+1)==0),:)=[]; w3=setdiff(w2,w1,'rows'); a2=w3(:,1); b2=w3(:,2); w4=w3; w4(find(I1(mm*(b2-3)+a2)==v(i)|I1(mm*(b2-3)+a2)==0),:)=[]; a3=w4(:,1); b3=w4(:,2); h1=I1(mm*(b3-3)+a3); h1=h1'; w5=w3; w5(find(I1(mm*(b2-1)+a2-2)==v(i)|I1(mm*(b2-1)+a2-2)==0),:)=[]; a4=w5(:,1); b4=w5(:,2); h2=I1(mm*(b4-1)+a4-2); h2=h2'; w6=w3; w6(find(I1(mm*(b2-1)+a2+2)==v(i)|I1(mm*(b2-1)+a2+2)==0),:)=[]; a5=w6(:,1); b5=w6(:,2); h3=I1(mm*(b5-1)+a5+2); h3=h3'; w7=w3; w7(find(I1(mm*(b2+1)+a2)==v(i)|I1(mm*(b2+1)+a2)==0),:)=[]; a6=w7(:,1); b6=w7(:,2); h4=I1(mm*(b6+1)+a6); h4=h4'; h=[h1 h2 h3 h4]; h=unique(h); if isempty(h) h=[v(r)]; end x=h(find(h~=0)); y=x'; [q1,q2]=min(abs(ag(i)-ag(y))); zgr=y(q2); I1(I1==v(i))=zgr; %将同一区域的像素值赋同一个值 end end figure,imshow(I1,[]),title('分水岭分割后的区域合并图像'); rgb1=label2rgb(I1); figure,imshow(rgb1),title('区域合并图像的彩色显示');