之前的文章讲述了肺结节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\MATLABb\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\MATLABb\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\MATLABb\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影像特征提取系统主要部分代码，部分工具包代码过长，无法贴出。
