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