基于改进空间约束高斯混合模型的图像分割技术

论文：Fast and Robust Spatially Constrained Gaussian Mixture Model for Image Segmentation

实现：

%%%%%主程序

clear all;
clc;

Img=imread('woman.jpg');
figure(1),imshow(Img),title('原图');
Img = rgb2gray(Img);
figure(2),imshow(Img),title('灰度图');
[height,width]=size(Img);
%Img=imnoise(Img,'gaussian',0,0.01);
Img=imnoise(Img,'salt & pepper',0.02);
figure(3),imshow(Img),title('原图加高斯噪声');
Img = im2double(Img);
y=Img(:);
K=3;

[u,afa,Sigma,fai,n] = initial_cal(K,y);       %输入:维度dim,图像长宽height,width,聚类数K,图像输入y
PI = repmat(afa',n,1);
img_line = zeros(n,1);

for m = 1:10
    [Z,G,u,Sigma,PI,J,fai,IDX] = update_cal(PI,fai,height,width,n,K,y);
    JJ(1,m) = J;
    for j=1:K
        c_label = find(IDX==j);
        img_line(c_label,:) = repmat(u(j,:),[size(c_label),1]);
    end
    img_out = reshape(img_line,height,width); 
    titlename = strcat('第',num2str(m),'次Gmm');
    figure(4),imshow(img_out),title(titlename);
    pause(0.001);
end
figure(5),plot(JJ);

%Function1 初始化

function [ u,afa,Sigma,fai,n ] = initial_cal( K,y )
%INITIAL_CAL Summary of this function goes here
%   Detailed explanation goes here
    n=size(y,1);
    [IDX,u] = kmeans(y,K);
    
    c_sum = zeros(K,1);
    cluster = zeros(n,K);
    Sigma = zeros(1,K);
    afa = zeros(K,1);


    for i = 1:n
        c_sum(IDX(i)) = c_sum(IDX(i)) + 1;
        cluster(c_sum(IDX(i)),IDX(i)) = y(i,:);   %%将每个坐标对对应的群记下来，并把该像素点的记下来
    end
    
    for j = 1:K
        clu=cluster(1:c_sum(j),j);                 
        Sigma(:,j) = cov(clu);                       %获得第j聚类的所有坐标rgb坐标所得到的积累
        afa(j)=c_sum(j)/n;                             %每个点j聚类概率
    end


    for j = 1:K
         fai(:,j) = ((2*pi)^(-0.5))*((Sigma(:,j)+realmin)^(-0.5))*exp(-0.5*(y - repmat(u(j,:),n,1))/(Sigma(:,j)+realmin).*(y - repmat(u(j,:),n,1)));
    end
end

%Function2 迭代更新

function [Z,G,u,Sigma,PI,J,fai,IDX] = update_cal(PI,fai,height,width,n,K,y)
%UPDATE-CAL Summary of this function goes here
%   Detailed explanation goes here



    f = PI.*fai;
    f_sum = sum(f,2);        %每个点 3个聚类的先验概率*高斯密度函数的和 
   
    for j=1:K
        Z(:,j) = f(:,j)./(f_sum); %获得后验概率Z
    end
    

    Z_temp = reshape(Z,height,width,K);
    PI_temp = reshape(PI,height,width,K);
    ZPI = Z_temp+PI_temp;


    H1 = zeros(2,width);
    H2 = zeros(height+4,2);    
 
    for j=1:K
        ZPI_m(:,:,j) = [H1;ZPI(:,:,j);H1];
        ZPI_img(:,:,j) = [H2 ZPI_m(:,:,j) H2];
    end


    for k=1:K
        for i=3:height+2            
            for j=3:width+2
                ZPI_AreaSum(i-2,j-2,k) = sum(ZPI_img(i-2:i+2,j-2,k))+sum(ZPI_img(i-2:i+2,j-1,k))+sum(ZPI_img(i-2:i+2,j,k))+sum(ZPI_img(i-2:i+2,j+1,k))+sum(ZPI_img(i-2:i+2,j+2,k));
            end
        end
    end
    
    b=6;
    for j=1:K
         G_img(:,:,j) = exp(b/50*(ZPI_AreaSum(:,:,j))); 
    end


    G = reshape(G_img,n,1,K);
    G = squeeze(G);                 % 更新G
    
    for j=1:K
         u(j,:) = sum(repmat(Z(:,j),1).*y)/sum(Z(:,j));  %update u

    end

    for j=1:K    
        Sigma(:,:,j) = sum(Z(:,j).*(y-repmat(u(j,:),n,1)).*(y-repmat(u(j,:),n,1)))/sum(Z(:,j));  %更新update Sigma
    end
    
    for j=1:K
        PI(:,j) = (Z(:,j)+G(:,j))./sum(Z+G,2);%更新PI
    end
    
    fai = zeros(n,K);
    
    for j = 1:K
         fai(:,j) = ((2*pi)^(-0.5))*((Sigma(:,j)+realmin)^(-0.5))*exp(-0.5*(y - repmat(u(j,:),n,1))/(Sigma(:,j)+realmin).*(y - repmat(u(j,:),n,1)));%更新fai
    end
    
    for j = 1:K
        jj(:,j) = Z(:,j).*(log(PI(:,j))+log(fai(:,k)))+G(:,j).*log(PI(:,j));
    end

    J = sum(sum(jj));%计算J
    
    for j = 1:K
        IDX(i,:) = find(Z(i,:)==max(max(Z(i,:))));  %更新IDX
    end
end

%%%%%%%%%%结束%%%%%%%%%%