MATLAB程序：用FCM分割脑图像

MATLAB程序：用FCM分割脑图像

作者：凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/

脑图像基础知识请看：脑图像；FCM算法介绍请看：聚类——FCM；数据来源：BrainWeb: Simulated Brain Database，只选取脑图像中的0、1、2、3类，其余类别设为0。本文用到的数据：Simulated Brain Database

1. MATLAB程序

FCM_image_main.m

function [accuracy,iter_FCM,run_time]=FCM_image_main(filename, num, K)
%num：第几层，K：聚类数
%[accuracy,iter_FCM,run_time]=FCM_image_main('t1_icbm_normal_1mm_pn0_rf0.rawb', 100, 4)
[data_load, label_load]=main(filename, num);  %原图像
[m,n]=size(data_load);
X=reshape(data_load,m*n,1);   %(m*n)*1
real_label=reshape(label_load,m*n,1)+ones(m*n,1);
Ground_truth(num, K);  %标准分割结果，进行渲染
t0=cputime;
[label_1,~,iter_FCM]=My_FCM(X,K);
[label_new,accuracy]=succeed(real_label,K,label_1);
run_time=cputime-t0;
label_2=reshape(label_new,m, n); 
rendering_image(label_2, K);  %聚类结果

main.m

function [read_new, mark]=main(filename, num)
%将真实脑图像中的0、1、2、3拿出来，其余像素为0.
%函数main(filename, num)中的第一个参数filename是欲读取的rawb文件的文件名，第二个参数num就是第多少张。
%例如：main('t1_icbm_normal_1mm_pn0_rf0.rawb',100)
mark=Mark('phantom_1.0mm_normal_crisp.rawb',num);
read=readrawb(filename, num);
[row, col]=size(read);
read_new=zeros(row, col);
for i=1:row   %行
    for j=1:col    %列
        if mark(i,j)==0
            read_new(i,j)=0;
        else
            read_new(i,j)=read(i,j);   %将第0、1、2、3类拿出来，其余类为0
        end
    end
end
%旋转90°并显示出来
figure(1);  
init_image=imrotate(read_new, 90); 
imshow(uint8(init_image)); 
title('原图像');

Mark.m

function mark=Mark(filename,num)
%将标签为1、2、3类分出来，其余为0，mark取值：0、1、2、3
%[mark_new,mark]=Mark('phantom_1.0mm_normal_crisp.rawb',90);
fp=fopen(filename);
temp=fread(fp, 181 * 217 * 181);
image=reshape(temp, 181 * 217, 181);   
images=image(:, num);
images=reshape(images, 181, 217);
mark_data=images;
fclose(fp);
mark=zeros(181,217);
%将第0、1、2、3类标签所在的坐标点拿出来，其余置0
for i=1:181
    for j=1:217
        if (mark_data(i,j)==1)||(mark_data(i,j)==2)||(mark_data(i,j)==3)
            mark(i,j)=mark_data(i,j);
        else
            mark(i,j)=0;
        end
    end
end

readrawb.m

function g = readrawb(filename, num)
% 函数readrawb(filename, num)中的第一个参数filename是欲读取的rawb文件的文件名，第二个参数num就是第多少张。
fid = fopen(filename);
% 连续读取181*217*181个数据，这时候temp是一个长度为181*217*181的向量。
% 先将rawb中的所有数据传递给temp数组，然后将tempreshape成图片集。
temp = fread(fid, 181 * 217 * 181);
% 所以把它变成了一个181*217行，181列的数组，按照它的代码，这就是181张图片的数据，每一列对应一张图。
% 生成图片集数组。图片集images数组中每一列表示一张图片。
images = reshape(temp, 181 * 217, 181);   
% 读取数组中的第num行，得到数组再reshape成图片原来的行数和列数：181*217。
image = images(:, num);
image = reshape(image, 181, 217);
g = image;
fclose(fid);
end

Ground_truth.m

function Ground_truth(num, K)
%标准分割结果
%Ground_truth(100, 4)
mark=Mark('phantom_1.0mm_normal_crisp.rawb',num);  %0、1、2、3
m=181;
n=217;
read_new=zeros(m,n);
mark=mark+ones(m, n);  %标签：1、2、3、4
for i=1:m   %行
    for j=1:n    %列
        for k=1:K
            if mark(i,j)==k
                read_new(i,j)=floor(255/K)*(k-1);               
            end
        end
    end
end
% 旋转90°并显示出来
figure(2)
truth_image=imrotate(read_new, 90); 
imshow(uint8(truth_image)); 
title('标准分割结果');

My_FCM.m

function [label_1,para_miu_new,iter]=My_FCM(data,K)
%输入K：聚类数
%输出：label_1:聚的类, para_miu_new:模糊聚类中心μ，responsivity:模糊隶属度
format long
eps=1e-8;  %定义迭代终止条件的eps
alpha=2;  %模糊加权指数，[1,+无穷)
T=100;  %最大迭代次数
fitness=zeros(T,1);
[data_num,~]=size(data);
count=zeros(data_num,1);  %统计distant中每一行为0的个数
responsivity=zeros(data_num,K);
R_up=zeros(data_num,K);
%----------------------------------------------------------------------------------------------------
%对data做最大-最小归一化处理
X=(data-ones(data_num,1)*min(data))./(ones(data_num,1)*(max(data)-min(data)));
[X_num,X_dim]=size(X);
%----------------------------------------------------------------------------------------------------
%随机初始化K个聚类中心
rand_array=randperm(X_num);  %产生1~X_num之间整数的随机排列
para_miu=X(rand_array(1:K),:);  %随机排列取前K个数，在X矩阵中取这K行作为初始聚类中心
% ----------------------------------------------------------------------------------------------------
% FCM算法
for t=1:T
    %欧氏距离，计算（X-para_miu）^2=X^2+para_miu^2-2*para_miu*X'，矩阵大小为X_num*K
    distant=(sum(X.*X,2))*ones(1,K)+ones(X_num,1)*(sum(para_miu.*para_miu,2))'-2*X*para_miu';
    %更新隶属度矩阵X_num*K
    for i=1:X_num
        count(i)=sum(distant(i,:)==0);
        if count(i)>0
            for k=1:K
                if distant(i,k)==0
                    responsivity(i,k)=1./count(i);
                else
                    responsivity(i,k)=0;
                end
            end
        else
            R_up(i,:)=distant(i,:).^(-1/(alpha-1));  %隶属度矩阵的分子部分
            responsivity(i,:)= R_up(i,:)./sum( R_up(i,:),2);
        end
    end
    %目标函数值
    fitness(t)=sum(sum(distant.*(responsivity.^(alpha))));
     %更新聚类中心K*X_dim
    miu_up=(responsivity'.^(alpha))*X;  %μ的分子部分
    para_miu=miu_up./((sum(responsivity.^(alpha)))'*ones(1,X_dim));
    if t>1  
        if abs(fitness(t)-fitness(t-1))<eps
            break;
        end
    end
end
para_miu_new=para_miu;
iter=t;  %实际迭代次数
[~,label_1]=max(responsivity,[],2);

succeed.m

function [label_new,accuracy]=succeed(real_label,K,id)
%输入K：聚的类，id：训练后的聚类结果，N*1的矩阵
N=size(id,1);   %样本个数
p=perms(1:K);   %全排列矩阵
p_col=size(p,1);   %全排列的行数
new_label=zeros(N,p_col);   %聚类结果的所有可能取值，N*p_col
num=zeros(1,p_col);  %与真实聚类结果一样的个数
%将训练结果全排列为N*p_col的矩阵，每一列为一种可能性
for i=1:N
    for j=1:p_col
        for k=1:K
            if id(i)==k
                new_label(i,j)=p(j,k);  %iris数据库，1 2 3
            end
        end
    end
end
%与真实结果比对，计算精确度
for j=1:p_col
    for i=1:N
        if new_label(i,j)==real_label(i)
                num(j)=num(j)+1;
        end
    end
end
[M,I]=max(num);
accuracy=M/N;
label_new=new_label(:,I);

rendering_image.m

function rendering_image(label,K)
%对分割结果进行渲染,4类,label:1、2、3、4
[m, n]=size(label);
read_new=zeros(m,n);
for i=1:m   %行
    for j=1:n    %列
        for k=1:K
            if label(i,j)==k
                read_new(i,j)=floor(255/K)*(k-1);               
            end
        end
    end
end
% 旋转90°并显示出来 
figure(3); 
cluster_image=imrotate(read_new, 90); 
imshow(uint8(cluster_image)); 
title('分割后');

2. 实验及结果

对T1模态、icmb协议下，切片厚度为1mm，噪声水平为7%，灰度不均匀水平为40%的第90层脑图像进行分割。因为FCM随机初始化，所以聚类结果会有偏差，结果受初始化影响比较大。

>> [accuracy,iter_FCM,run_time]=FCM_image_main('t1_icbm_normal_1mm_pn7_rf40.rawb', 90, 4)

accuracy =

   0.943783893881916


iter_FCM =

    25


run_time =

   1.937500000000000