基于混合高斯模型的图像分割方法

clc;clear;close all;

%逐幅打开图像列表中的图像，并计算其特征
[filenamefilepathfilterindex] = uigetfile（...
    {‘*.txt‘‘All Image Files（*.txt）‘;
    ‘*.mat‘‘mat-files（*.mat）‘;...
    
    ‘*.*‘‘All Image Files（*.*）‘}...
    ‘文件名列表‘‘untitled.mat‘...
    ‘multiselect‘‘off‘）;

if isequal（filename0）||isequal（filepath0）
    return;
else
   filefullpath = fullfile（filepathfilename）;
end

fid = fopen（filefullpath‘r‘）;

ffnum = fscanf（fid‘%d‘）;
fwnum = ffnum（1）;
fonum = ffnum（2）;
fprintf（1‘\nwaterNum=%d    otherNum=%d\n‘fwnumfonum）;
fnum = fwnum + fonum;

fpath = fgetl（fid）;
fprintf（1‘\n%s\n‘fpath）;

nxyA = zeros（fwnum1）;
cmA = zeros（33fwnum）;
mvA = zeros（3fwnum）;  
pixelnum = 0;

% learn_zoom = 1/5;
% test_zoom = 1/5;

%B6 pics
% learn_zoom = 1/10;
% test_zoom = 1/10;
learn_zoom_size = 100;
test_zoom_size = 200;

for k = 1:fwnum
    fName = fgetl（fid）; 
    ffName = fullfile（fpathfName）;
    fprintf（1‘\n%s\n‘ffName）;
    img = imread（ffName）;     %读取图像
    [imximyimd] = size（img）;
    max_xy = max（imximy）;
    if max_xy>learn_zoom_size
        learn_zoom = learn_zoom_size/max_xy;
        img = zoomImg（imglearn_zoom）;
    end
    simg = double（img）;    %转换数据格式
    Rd = simg（::1）;
    Gd = simg（::2）;
    Bd = simg（::3）;
    [nx ny] = size（Rd）;     %图像的行、列数
    nxy = nx * ny;           %图像的总象素数
    fprintf（1‘\n第%d测试图像的行数:%d   列数:%d   像素数:%d\n‘knxnynxy）;
%     subplot（221）;  imshow（img）;title（‘background RGB‘）;
    nxyA（k） = nxy;
    nxA（k） = nx;
    nyA（k） = ny;
    pixels = reshape（simgnxy3）;
    pixelsA（pixelnum+1:pixelnum+nxy:） = pixels;  
    pixelnum = pixelnum+nxy;
end
figure（1）;
subplot（221）;imshow（img）;title（‘background‘）;grid on;

x = pixelsA;
nbVar = size（x2）;
nn = 1000;
c = 1;
pixelsI = pixelsA（1:nn:）;
Mu（1:） = sum（pixelsI）/nn;
Sigma（::1） = cov（[pixelsI‘ pixelsI‘]‘）;
Priors（1） = 1;
% [Data_id Priors Mu Sigma] = EM_init_FCM（pixelsI c）;
% [Priors Mu Sigma] = EM（pixelsI Priors Mu Sigma）;
t = zeros（pixelnum1）;
t（1:nn） = 1;
beta = 1;
lamda = 3;
for j=1:pixelnum
    for k=1:c
        covt = reshape（Sigma（::k）nbVarnbVar）;
        incov = inv（covt）;
        dcov = det（covt）;
        logd = log（abs（dcov））;
        logp = log（Priors（k））;
        T（k） = beta*nbVar*lamda^2-2*logp+logd;
        a = x（j:）-Mu（k:）;
        tt = a*incov*a‘-2*logp+logd;
        if tt>T（k）
            t_3std（jk） = -1;
            N_ab（jk） = -10^10;
        else
            t_3std（jk） = 1;
            N_ab（jk） = -（0.5）*a*incov*a‘+logp+（-0.5）*logd;%后验概率值                        
        end 
    end        
    %判断第j个样本是哪一类：新的类别或已知类别中的一类
    [ir ic] = find（t_3std（j:）==1）;
    lnstd = length（ir）; 
    %若在所有类别的3sigma之外，则为新的类别；若在3sigma之内，则判定为后验概率最大的那一类
    if lnstd==0
            c = c+1;
            t（j） = c;
            if j                t（j+1:j+10） = c;
            end
    else
        N_ab_max = max（