肺結節CT影像特征提取（四）——肺結節CT影像特征提取MATLAB代碼實現

之前的文章講述了肺結節CT影像數據特征提取算法及基於MATLAB GUI設計的肺結節CT影像特征提取系統。本文將講述幾個主要部分的代碼實現，分別是預處理、灰度特征提取、紋理特征提取、形態特征提取數據。

一.預處理部分代碼

1、讀取肺結節CT數據和專家標記的mask數據

function [ sData ] = read_dcm_mask( dcmPath,maskPath,Ng )

function [ sData ] = read_dcm_mask( dcmPath,maskPath,Ng )
%read_dcm_mask.m 讀取dcm文件和mask文件為矩陣，為后期使用准備
%第一個程序
% DESCRIPTION:
%此函數處理dcm文件和mask文件
%1.設置dcm和mask文件所在路徑
%2.執行函數即可
%INPUTS:
%dcmPath:dcm文件所在路徑
%maskPath:mask文件所在路徑
%Ng:標准化的CT灰度級數
%OUTPUTS:
%sData:保存了volume和mask以及prepareVolume函數所需參數的一個cell結構

nowPath=cd;
mkdir(nowPath,'feature_extraction');%創建文件夾存放數據
dcmList=dir(dcmPath);               %獲取dcm文件列表
maskList=dir(maskPath);             %獲取mask文件列表
nDcm=size(dcmList,1);               %取得處理數據數目
%nMask=size(maskList,1);

%獲取prepareVolume函數需要的參數，創建cell結構的sData,保存volume和para.
dcm1Path=fullfile(dcmPath,dcmList(3).name);   %獲取第一個dcm文件
info=dicominfo(dcm1Path);                      %取得dcm文件部分信息用於參數設置
data.volume=[];
data.mask=[];
para.scanType='Other';
para.pixelW=info.PixelSpacing(1);
para.sliceS=info.SliceThickness;
para.R=1;
para.scale=info.PixelSpacing(1);
para.textType='Matrix';
para.quantAlgo='Lloyd';
para.Ng=Ng;
%sData={data,para}

%開始讀取數據
for i=1:nDcm-2
    dcmName=dcmList(i+2).name;
    dcmP=fullfile(dcmPath,dcmName);
    maskName=maskList(i+2).name;
    maskP=fullfile(maskPath,maskName);
    volume=dicomread(dcmP);
    data(i).volume=volume;
    mask1=imread(maskP);%醫生標記的ROI區域
    mask=im2bw(mask1);
    data(i).mask=mask;
end
% 
disp('數據讀取完畢！');
sData={data,para};
save sData.mat sData;
save info.mat info;

% file=fullfile(nowPath,'feature_extraction');
% movefile('sData.mat',file);

　2.獲取ROI區域數據

function [ROIdata] = getROI( sDataPath )

function [ROIdata] = getROI( sDataPath )
%function getROI.m 獲得ROI區域
% 第二個函數
%DESCRIPTION:
%讀取sData數據，對每個volume和mask進行預處理，調用prepareVolume函數
%INOUTS:sData數據的路徑
%OUTPUTS:保存ROIonly,levels,ROIbox,maskBox數據的ROIdata.

load(sDataPath);
nFile=size(sData{1},2);

%獲取參數
scanType=sData{1, 2}.scanType;
pixelW=sData{1, 2}.pixelW;
sliceS=sData{1, 2}.sliceS;
R=sData{1, 2}.R;
scale=sData{1, 2}.scale;
textType=sData{1, 2}.textType;
quantAlgo=sData{1, 2}.quantAlgo;
Ng=sData{2}.Ng;

for i=1:nFile
    volume=sData{1}(i).volume;  %獲得dcm數據
    mask=sData{1}(i).mask;      %獲取標記
    %調用prepareVolume函數獲得ROI區域數據
    [ROIonly,levels,ROIbox,maskBox] = prepareVolume(volume,mask,scanType,pixelW,sliceS,R,scale,textType,quantAlgo,Ng);
    ROI(i).ROIonly=ROIonly;
    ROI(i).levels=levels;
    ROI(i).ROIbox=ROIbox;
    ROI(i).maskBox=maskBox;
    fprintf('得到第%d組圖像ROI數據\n',i);
end
ROIdata=ROI;
save ROIdata.mat ROIdata;

二、提取灰度特征

　　肺結節區域對應的灰度直方圖，是表現了肺結節區域每一個像素出現的概率的圖像。對每張CT影像ROI區域進行計算，得到灰度直方圖。然后根據灰度統計的直方圖提取8個灰度特征，用這8個灰度特征來描述肺結節的特點，8個灰度特征分別如下表所示，是方差、標准差、最大像素值、最小像素值、偏離度、峰態、能量和熵。

　　代碼如下：

function [grayFeature] =get_gray_feature(ROIdataPath)
%function get_gray_feature.m  獲取圖像的灰度特征
%第三個函數
%DESCRIPTION:讀取ROIdata數據，得到灰度直方圖。提取灰度特征
%INPUTS:
%ROIdataPath:ROIdata的路徑
%OUTPUTS:
%grayFeature:灰度特征
% grayFeature(j).mean:均值
% grayFeature(j).variance：方差
% grayFeature(i).maxP：最大值
% grayFeature(i).minP：最小值
% grayFeature(j).skewness：偏離度
%grayFeature(j).kurtosis：峰態
%grayFeature(j).energy：能量
%grayFeature(j).entropy：熵
% 

load(ROIdataPath);%打ROIdata數據
mkdir(cd,'histogram');
%featurePath=fullfile('D:\wuProgram\MATLAB2014b\work\test\feature_extraction_box\feature_extration');
n=size(ROIdata,2);        %文件數量
Ng=size(ROIdata(1).levels,2);       % 灰度級數

%得到所有ROI區域的灰度直方圖
for i=1:n
    ROIonly=ROIdata(i).ROIonly;
    histData(i).H=my_hist(ROIonly,Ng);
    %name=strcat(num2str(i),'.jpg');
    %print('name.jpg');
end
save hist.mat histData;
disp('得到灰度直方圖');
%.1計算均值mean
for j=1:n
     mean=0;
    for k=0:Ng-1
        H=histData(j).H;
        mean=mean+k*H(k+1);
        grayFeature(j).mean=mean;
    end
end
disp('得到均值');

%2.計算方差variance
for j=1:n
    variance=0;
    for k=0:Ng-1
        mean=grayFeature(j).mean;
        H=histData(j).H;
        variance=variance+((k-mean).^2)*H(k+1);
        grayFeature(j).variance=variance;
        
    end
    grayFeature(j).deviation=sqrt(variance);
end
disp('得到方差');

%3.計算最大值最小值maxP,minP,
 for i=1:n
     ROIonly=ROIdata(i).ROIonly;
     maxP=max(max(ROIonly));
     grayFeature(i).maxP=maxP;
     minP=min(min(ROIonly));
     grayFeature(i).minP=minP;
 end
disp('得到最大值最小值');

%4.計算偏離度，峰態，能量
for j=1:n
    skewness=0;
    kurtosis=0;
    energy=0;
    for k=0:Ng-1
        H=histData(j).H;
        mean=grayFeature(j).mean;
        deviation= grayFeature(j).deviation;
        skewness=skewness+((k-mean).^3*H(k+1))./(deviation.^3);
        kurtosis=kurtosis+((k-mean).^4*H(k+1))./(deviation.^4);
        energy=energy+H(k+1).^2;
    end
    grayFeature(j).skewness=skewness;
    grayFeature(j).kurtosis=kurtosis-3;
    grayFeature(j).energy=energy;
    
    
end
disp('得到偏離度，峰態，能量');

%5.計算熵
entropy=0;
for j=1:n
    for k=0:Ng-1
        H=histData(j).H;
        if(H(k+1)==0)
            entropy=entropy+0;
        else
            entropy=entropy+H(k+1)*log2(H(k+1));
        end
    end
    grayFeature(j).entropy=entropy;
end
disp('得到熵');
save grayFeature.mat grayFeature;
%movefile('grayFeature.mat',featurePath);

三、提取紋理特征

　　紋理特征是一類人類視覺可以明顯感覺到的特征，同時也是圖像的一類重要特征，主要表現為像素在空間分布模式的描述，可以反映圖像表示的物體表面的粗糙度、光滑性、顆粒度、隨機性等性質。本文采用灰度共生矩陣（GCLM）來提取肺結節CT影像數據的紋理特征。

　　灰度特征是基於圖像矩陣的特點，利用數學方法構造的灰度共生矩陣，從灰度共生矩陣中得到圖像的信息，本文利用設計的肺結節CT影像特征提取系統對選取的9張肺結節CT影像進行特征提取，本文采用灰度共生矩陣方法（GCLM）提取肺結節的紋理特征，這里提取了能量、對比度、相關、熵、差分矩、和平均6個特征。

　　代碼如下(紋理特征利用了GCLM工具包)：

function [textureFeature] = get_texture_feature( ROIdataPath )
%get_texture_feature.m ：獲取紋理特征
%DESCRIPTION:
%調用以下兩個函數得到紋理特征
%[GLCM] = getGLCM(ROIonly,levels); 
%[glcmTextures] = getGLCMtextures(GLCM);
%IMPUTS:
%
%OUTPUTS:


load(ROIdataPath);%打ROIdata數據
% featurePath=fullfile('D:\wuProgram\MATLAB2014b\work\test\feature_extraction_box\feature_extration')
n=size(ROIdata,2);        %文件數量
%Ng=size(ROIdata(1).levels,2);       % 灰度級數

for i=1:n
    ROIonly=ROIdata(i).ROIonly;
    levels=ROIdata(i).levels;
    ROIbox=ROIdata(i).ROIbox;
    maskBox=ROIdata(i).maskBox;
    [GLCM] = getGLCM(ROIonly,levels); %調用getGLCM獲得GCLM矩陣
    [glcmTextures] = getGLCMtextures(GLCM);%調用getGCLMtextures函數獲得GCLM紋理
    textureFeature(i).glcmTextures=glcmTextures;
    fprintf('獲取第%d組數據GCLM紋理特征\n',i);
end
save textureFeature.mat textureFeature;
%featurePath=fullfile('D:\wuProgram\MATLAB2014b\work\test\feature_extraction_box\feature_extration');
%movefile('textureFeature.mat',featurePath);

四、提取形態特征數據

　　提取了肺結節的大小、周長、面積、重心和形狀參數特征，另外根據Hu不變矩算法提取了7個不變矩組作為形態特征。

1、基本形態特征提取代碼

function [geometryFeature] = get_geometry_feature(sDataPath)
%function get_geometry_feature.m
%purpose:
%獲取幾何參數，邊界長度perimeter、直徑diameter、面積area、重心orthocenter、形狀參數shape
%INPUTS:
%sDataPath:存儲有mask的數據文件位置
%OUTPUTS:
%geometryFeature:存儲有幾何特征的數據
load(sDataPath);
n=size(sData{1, 1},2);
for i=1:n
    mask=sData{1, 1}(i).mask;
    perimeter = get_perimeter(mask);
    [diameter,myarea] =get_diameter_area( mask );
    orthocenter = get_orthocenter( mask,myarea );
     shape = get_shape(perimeter ,myarea );
     geometryFeature(i).perimeter=perimeter;
     geometryFeature(i).diameter=diameter;
     geometryFeature(i).myarea=myarea;
     geometryFeature(i).orthocenter=orthocenter;
     geometryFeature(i).shape=shape;
     fprintf('得到第%d組形狀參數\n',i);
end

save geometryFeature.mat geometryFeature

2.形態特征提取子函數

function [perimeter] = get_perimeter(mask)
%function get_perimeter :獲取周長
%purpose：
%獲取mask中的周長
%INPUTS:
%mask：圈出區域
%OUTPUTS:
%perimeter：周長

m_edge=edge(mask);
edgeSpot=find(m_edge==1); %獲取邊界坐標數組
perimeter=size(edgeSpot,1);
end

function [diameter,myArae] =get_diameter_area( mask )
%function get_diameter_area:獲取mask的直徑和面積
%description：
%輸入mask，值不為0的地方為感興趣區域，求其直徑和面積
%INPUTS:
%mask:處理的矩陣
%OUTPUTS:
%diameter:直徑
%area：面積

[m,n]=size(mask);
myArae=0;
diaTemp=[];
k=1;
for i=1:m
    for j=1:n
        if(mask(i,j)~=0)
            myArae=myArae+1;   %計算面積
            diaTemp(k,1)=i;  %區域存儲坐標
            diaTemp(k,2)=j;
            k=k+1;
        else
            continue;
        end
    end
end

n=size(diaTemp);
length=[];
k=1;
for i=1:n
    for j=i+1:n
        length(k)=sqrt((diaTemp(j,1)-diaTemp(i,1)).^2+(diaTemp(j,2)-diaTemp(i,2)).^2);%計算任意兩個坐標間的距離
        k=k+1;
    end
end

diameter=max(length);%得到最大距離，即直徑
end

function [orthocenter] = get_orthocenter( mask,area )
%function get_orthocenter:獲取mask重心
%description：
%獲取mask區域重心
%INPUTS:
%mask:需要計算的的圖像
%OUTPUTS:
%orthocenter:重心

[m,n]=size(mask);
xTemp=0;
yTemp=0;

for i=1:m
    for j=1:n
        if(mask(i,j)~=0)
            xTemp=xTemp+i;
            yTemp=yTemp+j;
        else
            continue;
        end
    end
end
orthocenter.x=ceil(xTemp/area);
orthocenter.y=ceil(yTemp/area);

end

function [shape] = get_shape(perimeter ,area )
%function get_shape:獲取形狀參數
%description：
%獲取mask形狀參數
%INPUTS:
%perimeter ,area區域的mask需要計算的的圖像的周長，面積
%OUTPUTS:
%shape:形狀參數

shape=(perimeter.^2)/(4*pi*area);
end

　　上述為肺結節CT影像特征提取系統主要部分代碼，部分工具包代碼過長，無法貼出。