作者: 凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/
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);
- for k=1:K
- if mark(i,j)==k
- read_new(i,j)=floor(255/K)*(k-1);
- end
- end
- end
- end
- 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)
- fitness=zeros(T,1);
- [data_num,~]=size(data);
- responsivity=zeros(data_num,K);
- R_up=zeros(data_num,K);
- %----------------------------------------------------------------------------------------------------
- X=(data-ones(data_num,1)*min(data))./(ones(data_num,1)*(max(data)-min(data)));
- [X_num,X_dim]=size(X);
- %----------------------------------------------------------------------------------------------------
- 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
- responsivity(i,:)= R_up(i,:)./sum( R_up(i,:),2);
- end
- end
- 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;
- [~,label_1]=max(responsivity,[],2);
- succeed.m
- function [label_new,accuracy]=succeed(real_label,K,id)
- for i=1:N
- for j=1:p_col
- for k=1:K
- if id(i)==k
- 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)
- [m, n]=size(label);
- read_new=zeros(m,n);
- for k=1:K
- if label(i,j)==k
- read_new(i,j)=floor(255/K)*(k-1);
- end
- end
- end
- end
- figure(3);
- cluster_image=imrotate(read_new, 90);
- imshow(uint8(cluster_image));
- title('分割后');
- >> [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
来源: https://www.cnblogs.com/kailugaji/p/10756673.html