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
