FCM 图像分割

%%%%%%%%%%%%%%% FCM算法分割图像 %%%%%%%%%%%%%%
function clusterResult = FCM(imagePath, C, V, M, iter, epsm)
% 模糊C均值(FCM)聚类算法分割图像
% 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 < 3
    V = [0 255];
end

if nargin < 2
    C = 2;
end

% 读入图像及其信息
I = imread(imagePath);
figure, imshow(I);
title('原图像');
[row col] = size(I);
grayHist = imhist(I);
figure, imhist(I);
title('直方图');
histProb = grayHist / (row * col);
len = length(histProb);

tic
% FCM迭代过程
cnt = 0;
while(cnt < iter)
% 计算隶属度函数(注意要特殊考虑某个像素点和聚类中心一样的情况)
    for i = 1 : len
        flag = 0;
        for j = 1 : C
            if i == V(j)
                U(j, i) = 1.0;
                if j == 1
                    U(j + 1 : C, i) = 0.0;
                elseif j == C
                    U(1 : C - 1, i) = 0.0;
                else
                    U(1 : j - 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
toc

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