function [IX2] = fcm(IM);
%IM是輸入的原圖像
%IX2是分類結果
[maxX,maxY]=size(IM);

IMM=cat(3,IM,IM,IM);
%初始化聚類中心（3類）
cc1=8;
cc2=100;
cc3=200;

ttFcm=0;

while(ttFcm<15)
    
    ttFcm=ttFcm+1
    
    c1=repmat(cc1,maxX,maxY);
    c2=repmat(cc2,maxX,maxY);
    c3=repmat(cc3,maxX,maxY);
    c=cat(3,c1,c2,c3);
    
    ree=repmat(0.000001,maxX,maxY);
    ree1=cat(3,ree,ree,ree);
    
    distance=IMM-c;
    distance=distance.*distance+ree1;
    
    daoShu=1./distance;
    daoShu2=daoShu(:,:,1)+daoShu(:,:,2)+daoShu(:,:,3);
%計算隸屬度u    
    distance1=distance(:,:,1).*daoShu2;
    u1=1./distance1;
    distance2=distance(:,:,2).*daoShu2;
    u2=1./distance2;
    distance3=distance(:,:,3).*daoShu2;
    u3=1./distance3;
%計算聚類中心z    
    ccc1=sum(sum(u1.*u1.*IM))/sum(sum(u1.*u1));
    ccc2=sum(sum(u2.*u2.*IM))/sum(sum(u2.*u2));
    ccc3=sum(sum(u3.*u3.*IM))/sum(sum(u3.*u3));
    
    tmpMatrix=[abs(cc1-ccc1)/cc1,abs(cc2-ccc2)/cc2,abs(cc3-ccc3)/cc3];
    
    pp=cat(3,u1,u2,u3);
    
    for i=1:maxX
        for j=1:maxY
            if max(pp(i,j,:))==u1(i,j)
                IX2(i,j)=1;
            elseif max(pp(i,j,:))==u2(i,j)
                IX2(i,j)=2;
            else
                IX2(i,j)=3;
            end
        end
    end
%判結束條件    
    if max(tmpMatrix)<0.0001
        break;
    else
        cc1=ccc1;
        cc2=ccc2;
        cc3=ccc3;
    end
    
    for i=1:maxX
        for j=1:maxY
            if IX2(i,j)==3
                IMMM(i,j)=240;
            elseif IX2(i,j)==2
                IMMM(i,j)=130;
            else
                IMMM(i,j)=20;
            end
        end
    end
%顯示每次聚類分割結果    
    figure(1);
    imshow(uint8(IMMM));
    
end

for i=1:maxX
    for j=1:maxY
        if IX2(i,j)==3
            IMMM(i,j)=240;
        elseif IX2(i,j)==2
            IMMM(i,j)=130;
        else
            IMMM(i,j)=20;
        end
    end
end
%顯示最終分類結果
IMMM=uint8(IMMM);
figure(2);
imshow(IMMM);

end