图像分割

摘 要：本实验采用Sobel算子对图像进行边缘检测,编写一个自动估计阈值的全局阈值化程序，对输入图像进行阈值分割，输出为分割后的二值图。另外，本实验用光渐变算法对图像进行污染，并设计了一个阈值化算法对图像进行阈值分割。由于采用自动估计阈值的全局阈值化算法以及Otsu算法都没有得到较好的分割结果，此实验便采用了模糊C均值（FCM）算法来实现分割，得到了相对
好一些的效果。

function [imt,G] = Sobelgrads(ima,T)
%功能：计算Sobel梯度
%ima为输入图像，T为阈值
%imt为输出的二值图像，A为Sobel梯度数据
%M1,M2为Sobel模板
M1=[-1 -2 -1;0 0 0;1 2 1];
M2=M1';
imt=zeros(size(ima));
gx=imfilter(ima,M1);
gy=imfilter(ima,M2);
G=abs(gx)+abs(gy);
imt(G>=T)=1;
end
%-----------------------------------------
function pro1001
%计算Ｓｏｂｅｌ梯度
ima=imread('Fig9.39(a).jpg');
T=160;
[imt,G1] = Sobelgrads(ima,T);
figure(1);
subplot(231); imshow(ima); title('原图像');
subplot(232); imshow(uint8(G1)); title('Sobel梯度图');
subplot(233); imshow(imt); title('二值图像');
figure(2);imhist(G1);
%先平滑处理，在计算Sobel梯度
H = [1,2,1;2,4,2;1,2,1]/16;
imb = imfilter(ima,H);
figure(1);
subplot(234); imshow(imb); title('平滑处理后图像');
[imt2,G2] = Sobelgrads(imb,T);
subplot(235); imshow(uint8(G2)); title('平滑处理后Sobel梯度图');
subplot(236); imshow(imt2); title('平滑处理后二值图像');
end
%-----------------------------------------------------
function [imt] = BGT(ima,t)
%功能：自动估值阈值的全阈值化程序
%ima为输入图像
%ｔ用来控制迭代的次数
%imt为分割后的二值图
[m,n]=size(ima);
%用图像的平均灰度初始T
T=sum(sum(ima))/(m*n);
imt=zeros(m,n);
his=imhist(ima);
cnt=1;
while(cnt<=100)
    r1=(find(his>T));
    r2=(find(his<=T));
    m1=sum(r1)/length(r1);
    m2=sum(r2)/length(r2);
    Tresult=0.5*(m1+m2);
    if (abs(Tresult-T)<t)
        break;
    else
        T=Tresult;    
    end
    cnt=cnt+1;
end
display(cnt);
display(Tresult);
for i=1:m
    for j=1:n
     if ima(i,j)>T
         imt(i,j)=1;
     else
         imt(i,j)=0;
     end
    end
end
imt=logical(imt);

end
%---------------------------------------------
function pro1002
%用阈值化程序对图像中的目标进行分割
ima=imread('Fig10.36(b).jpg');
[m,n]=size(ima);
imt=BGT(ima,1.0e-2);
figure;subplot(121);imshow(ima);title('原图');
subplot(122);imshow(imt);title('阈值处理后的图像');
%用光渐变算法污染图像
imb=zeros(m,n);
r=linspace(0.2,0.6,n);
for i=1:n
    imb(:,i)=r(i);
end
figure;subplot(311);imshow(ima);title('带噪声的图像');
subplot(312);imshow(imb);
title('在【0.2,0.6】范围内的灰度斜坡图像');
imc=immultiply(double(ima),double(imb));
subplot(313);imshow(uint8(imc));title('图（a）和图（b）的乘积');
figure;
subplot(311);imhist(ima);
subplot(312);imhist(imb);
subplot(313);imhist(uint8(imc));imwrite(uint8(imc),'imc.jpg');
%设计一种阈值化算法对图像进行阈值分割
%先平滑
H = [1,2,1;2,4,2;1,2,1]/16;
imd = imfilter(imc,H);
imwrite(uint8(imd),'imd.jpg')
%采用模糊C均值算法对图像进行分割（FCM）
%调用编写的FCM算法，聚类中心设为两个，初始聚类中心设为【0.255】
FCM('imd.jpg',2,[0 277]);
end
%----------------------------------------------------

function [ clusterResult,im ] = FCM( imagePath,C,V,M,iter,epsm )
% clusterResult = FCM(imagePath, C, V, M, iter, epsm)
% Example: clusterResult =? FCM('E:\Image\lena.bmp')
% clusterResult =? FCM('E:\Image\lena.bmp',3,[0 127 255])
% Input:
% imagePath: 图像路径
% C: 类别数，缺省值为2
% V: 初始化聚类中心,缺省值为[0 255]
% M: 加权指数,缺省值为2
% iter: 迭代次数,缺省值为100
%epsm: 迭代停止阈值,缺省值为1.0e-2
% Output:
% clusterResult: 聚类中心结果
% Note:
% C的值要与V的初始化聚类中心个数相同
%设定缺省值
if nargin<6
    epsm=1.0e-2;
end
if nargin<5
    iter=100;
end
if nargin<4
    M=2;
end
if nargin<2
    C=2;
end

%读入图像
I=imread(imagePath);
figure;subplot(121);imshow(I);title('原图像');
[m, n]=size(I);
im=zeros(m,n);
grayHist=imhist(I);
%subplot(132);imhist(I);title('直方图');
%直方图归一化
histProb=grayHist/(m*n);
len=length(histProb);

cnt=0;
while(cnt<iter)
    %计算隶属度函数（注意要特殊考虑某个像素点和聚类中心一样的情况）
    for i=1:len
        flag=0;
        for j=1:C
            if i==V(j)
                U(i,j)=1.0;
                if j==1
                    U(j+1:C,i)=0.0;
                elseif j==C
                    U(1:C-1,i)=0.0;
                    U(j+1:C,i)=0.0;
                end
                flag=1;
                break;
            end
        end

        if flag==0;
            u=(1.0./((i-V).^2)).^(1.0/(M-1));
            uSum=sum(u);
            U(1:C,i)=u'/uSum;
        end
    end
    %计算更新各类聚类中心
    for j=1:C
        i=linspace(1,len,len);
        v=sum(histProb'.*i.*(U(j,:).^M));
        vSum=sum(histProb'.*(U(j,:).^M));
        if vSum==0
            clusterResult(j)=0;
        else
            clusterResult(j)=v/vSum;
        end
    end
    %计算误差并判断算法迭代是否停止
    diff=sum((clusterResult-V).^2);
    if diff<=epsm
        break;
    else
    V=clusterResult;
    end
    cnt=cnt+1;
end

%分割图像
for i=1:m
    for j=1:n
        temp=(double(I(i,j))-clusterResult).^2;
        [~, pos]=min(temp);
        im(i,j)=pos*255/C;
    end
end
subplot(122);imshow(uint8(im));
title('分割后的图像');
disp('迭代次数：iterTimes=');
disp(cnt);
disp('clusterResult');
disp(clusterResult);
end