MATLAB程序：用FCM分割脑图像

成长路上除了

1. MATLAB程序
FCM_image_main.m

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

复制代码

main.m

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

复制代码

Mark.m

1. function mark=Mark(filename,num)
   
2. %将标签为1、2、3类分出来，其余为0，mark取值：0、1、2、3
   
3. %[mark_new,mark]=Mark('phantom_1.0mm_normal_crisp.rawb',90);
   
4. fp=fopen(filename);
   
5. temp=fread(fp, 181 * 217 * 181);
   
6. image=reshape(temp, 181 * 217, 181);  
   
7. images=image(:, num);
   
8. images=reshape(images, 181, 217);
   
9. mark_data=images;
   
10. fclose(fp);
    
11. mark=zeros(181,217);
    
12. %将第0、1、2、3类标签所在的坐标点拿出来，其余置0
    
13. for i=1:181
    
14.     for j=1:217
    
15.         if (mark_data(i,j)==1)||(mark_data(i,j)==2)||(mark_data(i,j)==3)
    
16.             mark(i,j)=mark_data(i,j);
    
17.         else
    
18.             mark(i,j)=0;
    
19.         end
    
20.     end
    
21. end

复制代码

readrawb.m

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

复制代码

Ground_truth.m

1. function Ground_truth(num, K)
   
2. %标准分割结果
   
3. %Ground_truth(100, 4)
   
4. mark=Mark('phantom_1.0mm_normal_crisp.rawb',num);  %0、1、2、3
   
5. m=181;
   
6. n=217;
   
7. read_new=zeros(m,n);
   
8. mark=mark+ones(m, n);  %标签：1、2、3、4
   
9. for i=1:m   %行
   
10.     for j=1:n    %列
    
11.         for k=1:K
    
12.             if mark(i,j)==k
    
13.                 read_new(i,j)=floor(255/K)*(k-1);              
    
14.             end
    
15.         end
    
16.     end
    
17. end
    
18. % 旋转90°并显示出来
    
19. figure(2)
    
20. truth_image=imrotate(read_new, 90);
    
21. imshow(uint8(truth_image));
    
22. title('标准分割结果');

复制代码

My_FCM.m

1. function [label_1,para_miu_new,iter]=My_FCM(data,K)
   
2. %输入K：聚类数
   
3. %输出：label_1:聚的类, para_miu_new:模糊聚类中心μ，responsivity:模糊隶属度
   
4. format long
   
5. eps=1e-8;  %定义迭代终止条件的eps
   
6. alpha=2;  %模糊加权指数，[1,+无穷)
   
7. T=100;  %最大迭代次数
   
8. fitness=zeros(T,1);
   
9. [data_num,~]=size(data);
   
10. count=zeros(data_num,1);  %统计distant中每一行为0的个数
    
11. responsivity=zeros(data_num,K);
    
12. R_up=zeros(data_num,K);
    
13. %----------------------------------------------------------------------------------------------------
    
14. %对data做最大-最小归一化处理
    
15. X=(data-ones(data_num,1)*min(data))./(ones(data_num,1)*(max(data)-min(data)));
    
16. [X_num,X_dim]=size(X);
    
17. %----------------------------------------------------------------------------------------------------
    
18. %随机初始化K个聚类中心
    
19. rand_array=randperm(X_num);  %产生1~X_num之间整数的随机排列
    
20. para_miu=X(rand_array(1:K),:);  %随机排列取前K个数，在X矩阵中取这K行作为初始聚类中心
    
21. % ----------------------------------------------------------------------------------------------------
    
22. % FCM算法
    
23. for t=1:T
    
24.     %欧氏距离，计算（X-para_miu）^2=X^2+para_miu^2-2*para_miu*X'，矩阵大小为X_num*K
    
25.     distant=(sum(X.*X,2))*ones(1,K)+ones(X_num,1)*(sum(para_miu.*para_miu,2))'-2*X*para_miu';
    
26.     %更新隶属度矩阵X_num*K
    
27.     for i=1:X_num
    
28.         count(i)=sum(distant(i,:)==0);
    
29.         if count(i)>0
    
30.             for k=1:K
    
31.                 if distant(i,k)==0
    
32.                     responsivity(i,k)=1./count(i);
    
33.                 else
    
34.                     responsivity(i,k)=0;
    
35.                 end
    
36.             end
    
37.         else
    
38.             R_up(i,:)=distant(i,:).^(-1/(alpha-1));  %隶属度矩阵的分子部分
    
39.             responsivity(i,:)= R_up(i,:)./sum( R_up(i,:),2);
    
40.         end
    
41.     end
    
42.     %目标函数值
    
43.     fitness(t)=sum(sum(distant.*(responsivity.^(alpha))));
    
44.      %更新聚类中心K*X_dim
    
45.     miu_up=(responsivity'.^(alpha))*X;  %μ的分子部分
    
46.     para_miu=miu_up./((sum(responsivity.^(alpha)))'*ones(1,X_dim));
    
47.     if t>1
    
48.         if abs(fitness(t)-fitness(t-1))<eps
    
49.             break;
    
50.         end
    
51.     end
    
52. end
    
53. para_miu_new=para_miu;
    
54. iter=t;  %实际迭代次数
    
55. [~,label_1]=max(responsivity,[],2);

复制代码

succeed.m

1. function [label_new,accuracy]=succeed(real_label,K,id)
   
2. %输入K：聚的类，id：训练后的聚类结果，N*1的矩阵
   
3. N=size(id,1);   %样本个数
   
4. p=perms(1:K);   %全排列矩阵
   
5. p_col=size(p,1);   %全排列的行数
   
6. new_label=zeros(N,p_col);   %聚类结果的所有可能取值，N*p_col
   
7. num=zeros(1,p_col);  %与真实聚类结果一样的个数
   
8. %将训练结果全排列为N*p_col的矩阵，每一列为一种可能性
   
9. for i=1:N
   
10.     for j=1:p_col
    
11.         for k=1:K
    
12.             if id(i)==k
    
13.                 new_label(i,j)=p(j,k);  %iris数据库，1 2 3
    
14.             end
    
15.         end
    
16.     end
    
17. end
    
18. %与真实结果比对，计算精确度
    
19. for j=1:p_col
    
20.     for i=1:N
    
21.         if new_label(i,j)==real_label(i)
    
22.                 num(j)=num(j)+1;
    
23.         end
    
24.     end
    
25. end
    
26. [M,I]=max(num);
    
27. accuracy=M/N;
    
28. label_new=new_label(:,I);

复制代码

rendering_image.m

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

复制代码

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

1. >> [accuracy,iter_FCM,run_time]=FCM_image_main('t1_icbm_normal_1mm_pn7_rf40.rawb', 90, 4)
   
2. 
   
3. accuracy =
   
4. 
   
5.    0.943783893881916
   
6. 
   
7. 
   
8. iter_FCM =
   
9. 
   
10.     25
    
11. 
    
12. 
    
13. run_time =
    
14. 
    
15.    1.937500000000000

复制代码

作者：凯鲁嘎吉