MATLAB程序：用FCM分割腦圖像 聚類——FCM

時間  2019-12-10
標籤 matlab 程序 fcm 分割 圖像 欄目 MATLAB 简体版
原文   原文鏈接

MATLAB程序：用FCM分割腦圖像

做者：凱魯嘎吉 - 博客園 http://www.cnblogs.com/kailugaji/html

腦圖像基礎知識請看：腦圖像；FCM算法介紹請看：聚類——FCM；數據來源：BrainWeb: Simulated Brain Database，只選取腦圖像中的0、一、二、3類，其他類別設爲0。本文用到的數據：Simulated Brain Databaseweb

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、一、二、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、一、二、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)
%將標籤爲一、二、3類分出來，其他爲0，mark取值：0、一、二、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、一、二、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、一、二、3
m=181;
n=217;
read_new=zeros(m,n);
mark=mark+ones(m, n);  %標籤：一、二、三、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:一、二、三、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

                                  

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程