一、引言

皮肤癌的发生在世界范围内都比较普遍，一直为皮肤科学等相关领域的学者所重视，而皮肤癌的发病频率最近这些年有逐年上升的趋势。临床试验表明，只要尽可能早地对色素性皮肤病进行早期精确诊断，大多数病人都能得到根治。因此，进行皮肤癌前期检测是皮肤癌防治的重要方法。当前皮肤癌分析的金标准是使用活体组织切片检测，但这种有创伤性的技术容易造成病人生理伤害。基于计算机辅助图像分析的方法是当前皮肤癌无创伤性检测的一个新技术，其主要的流程包括：皮肤图像获取、皮损区域的自动分割、特征提取与计算、分类识别并诊断。这里面的一个关键步骤就是皮损区域的分割。在数字图像处理和分析领域，分割是将图像细分为构成它的子区域或者对象，从而将感兴趣的区域或者对象检测出来。

Matlab是Matrix Laboratory的缩写，对于科学计算来说是一种高性能的语言。Matlab 发展到现在,已经成为世界上最流行的科学与工程计算中面向对象的交互编程工具,以可靠的矩阵运算为基础,具有图形绘制、数据处理、系统分析、信号和图像处理、符号运算功能,并且针对不同专业领域,研发了各种工具箱,其中包含了当今流行的新算法,广泛应用于自动控制、信号与图像处理。Matlab提供的功能强大的图像处理工具箱，使得图像处理和分析的工程技术人员可以方便快速地用Matlab来实践自己的算法。在本文中，我们就将使用Matlab图像处理工具箱来处理色素性皮损区域的自动分割问题。

二、皮损区域自动分割的方法

1.皮肤图像的读入和预处理

为了得到类似于VC++中的交互式文件打开的效果，可以调用Matlab提供的标准文件打开对话框函数uigetfile，该函数有8种调用形式，常用的一种是：[fileName,pathName] = uigetfile(‘filterSpec',‘dialogTitle')。其中fileName和pathName是返回的文件名和路径名；filterSpec指定出现在文件打开对话框中文件类型过滤器字符串，而dialogTitle是出现在文件打开对话框窗口的标题，用于给出一些必要的提示信息。将皮肤图像从硬盘读入到内存的实现代码如下：

[fileName,pathName] = uigetfile(‘.jpg',‘Please select a jpg format file'); %打开用于选取文件的标准对话框，此处默认打开的文件类型是jpg格式的图像。

fileName = strcat(pathName,fileName); %通过字符串连接函数获得完整的文件名

img = imread(fileName);%读入选择的图像文件到内存

由于后面的分割算法是针对灰度图像而不是彩色图像，所以在预处理阶段应该调用rgb2gray函数实现彩色图像到8位灰度图像的转换。这时，需要考虑读入的图像是不是彩色图像，如果是，那么我们就将其转换到8位灰度图像，否则就不进行处理。在Matlab中，只要原来的图像是彩色图像（注意，不是8位彩色位图图像），无论是JPG格式，还是16、24或者32位的BMP格式图像，通过imread读入到内存中，都是MxNx3形式的三维数组。因此，在Matlab中，通过条件语句if(size(img,3) == 3)是判断一副图像是否为16位以上彩色图像的常用技巧。整个转换过程的实现代码如下：

imgBak = img;%备份原始图像数据，以备后面构造彩色叠加效果图

if (size(img,3) == 3)%如果是彩色图像，则将其转化为灰度图像进行处理

    img = rgb2gray(img);%rgb2gray只有一个输入参数，表示待转换的彩色图像

end

此外，由于成像设备内在原因，获取的皮肤图像中往往会含有一些不可预测的随机信号，即通常所说的噪声。噪声的存在很大程度上影响到后续步骤的分析结果，因此对噪声的抑制显得尤为重要。Matlab中提供了很多噪声滤波器，比如均值滤波、加权均值滤波以及中值滤波等。考虑到中值滤波器既有抑制噪声的能力，同时又能较好地保留边缘信息，在我们的实现中，选用了中值滤波器，其调用形式的实现代码如下：

img = medfilt2(img,[7 7]);%使用中值滤波器进行噪声的滤除，第一个输入参数表示待滤波的图像，第二个参数表%示滤波窗口的大小，一般窗口越大，噪声抑制的效果越好，但是对边缘的破坏也会比较厉害，这里通过大量的%试验发现，滤波窗口取[7 7]对我们获取的皮肤图像噪声抑制效果比较理想。

2. 基于Otsu方法的初步分割

由于实现的直观性和简单性，图像阈值处理在图像分割应用中占有重要地位。对图像f(x,y)确定一个阈值T，任意满足f(x,y) >= T的点称为对象点，其他点成为背景点。换言之，阈值处理后的二值图像g(x,y)定义为：

                                                  （1）

这种方法称为全局阈值选取。在这里，阈值T的选取将是关键。工具箱提供了一个函数graythresh，该函数使用经典的Otsu的算法法来计算阈值T。函数graythresh的调用方法是：

T = graythresh(f),其中f是输入的图像，T是产生的阈值。

为了分割图像，我们使用函数BW = im2bw(f,T)，T就是graythresh函数的输出结果。函数im2bw使用式（1）的原理来对输入的图像进行二值化处理。基于Otsu的阈值初步分割的代码如下所示，其分割效果参考图1：

imgLevel = graythresh(img);%使用Otsu方法获得阈值，该值取值已经归一化到[0，1]之间

imgMask = im2bw(img,imgLevel);%根据全局阈值分割的思想进行分割，得到二进制的图像imgMask 

图1 (a)是原始rgb彩色图像；(b)是对应的灰度图像，通过函数rgb2gray来实现彩色图像到灰度的图像的转化；(c)使用使用全局阈值化分割后的结果（手动标识的红色线框内含有空洞）

3. 基于数学形态学的皮损区域边缘精定位

从图1(c)中我们可以发现分割所得的色素性皮损区域含有空洞，为了得到没有空洞的皮损区域，我们可以使用Matlab中的数学形态学填充函数imfill来实现，然后为了去除边界上一些毛刺，可以对图1（c）反相变换操作后再进行形态学闭运算。图2（a）是填充空洞后再进行形态学闭运算的结果。整个形态学变换过程的实现代码如下：

imgMask = ~imgMask;%对二进制图像进行去反操作，这里对图像取反的原因是，形态学填充函数imfill认为白%色像素包围的黑色区域才是空洞，因此为了正确进行空洞的填充，需要对图像进行进行取反

imgMask = imfill(imgMask,'holes');%填充图像中的空洞

se = strel('disk',6);%创建形态学圆盘结构元素

BWc = imclose(imgMask,se);%进行形态学闭操作以光滑部分凹陷的边界

图2 （a）填充空洞后再进行形态学闭运算的结果；(b)中红色标示出（a）中的非皮损区域；（c）检测到的正确的皮损区域；(d)最终提取的皮损区域边界

4.去除小的干扰区域并最终提取出皮损区域边界

从图2（a）我们可以发现，图中除了分割出的皮损区域外，还有一些小的非皮损区域（参考图2（b））。为了去除这些小的区域，我们只需要保留最大的白色区域块,使用下面5个步骤就可以达到目的：

Step1 使用bwlabel函数对4连通区域进行顺序标识，函数的调用方式如下：

[labeled, numObjects] = bwlabel(finalMask,4)；其中finalMask表示图2(a)对应的二值图像（即上段代码中的BWc）。其作用是对finalMask中4连通的区域进行顺序编号，编号从0开始到，并且finalMask中所有的背景像素（值为0）都统一标识为0；而前景像素（值非0）按区域连通性顺序编号。比如，下面的一个示例矩阵demoMatrix（关于这个例子的每步结果，可以通过运行源代码目录下bwlabelDemo.m函数来查看）:

demoMatrix = [ 1 0 0 0 0 1;

               1 1 1 0 0 1;

               1 1 1 0 0 1;

               1 0 0 1 1 0;

               0 0 0 0 0 0];

使用上面的语句得到的labeled将是：

    labeled = [ 1 0 0 0 0 3;

                1 1 1 0 0 3;

    1 1 1 0 0 3;

    1 0 0 2 2 0;

    0 0 0 0 0 0];

Step2 使用regionprops函数对labeled矩阵进行区域属性的计算，函数的调用方式如下：

data = regionprops(labeled,‘basic')；其返回值data是一个结构体；

Step3 获取结构体data中所有4连通区域的面积

dataArea = [data.Area]；面积用4连通区域的像素数来刻画，其中labeled中值为0的背景区域不参与计算。对上面的示例矩阵demoMatrix所对应的labeled，其dataArea = [ 8 2 3]。

Step4 获取具有最大面积的4连通区域在labeled中的标识：

bigRegion = find(dataArea == max(dataArea))；其中find查找到dataArea向量中最大值在该向量中的索引值，对示例矩阵demoMatrix，bigRegion等于1，因为标识为1的区域具有最大的4连通面积值。

Step5 获得原始finalMask中具有最大面积的区域

mask = (labeled == bigRegion)；对示例的demoMatrix，其对应的mask是：

mask = [ 1 0 0 0 0 0;

             1 1 1 0 0 0;

 1 1 1 0 0 0;

 1 0 0 0 0 0;

         0 0 0 0 0 0];

可以看出，最后的mask只保留有最大面积的4连通区域，而别的地方均为0。

由于分割一般只需要获得区域的边界像素，而获取二值图像的边界可以使用二值图像的形态学函数：finalMask = bwmorph(mask,‘remove')；该函数判断一个像素其4连通区域是否都是值1，如果是，则将该像素的值置为0，从而保留皮损区域的最外层边界。

上述步骤的Matlab实现代码如下：

finalMask = BWc;

[labeled, numObjects] = bwlabel(finalMask,4);%对二进制图像的4连通区域进行顺序标识

data = regionprops(labeled,'basic');%获得标识后的图像的区域基本属性

dataArea = [data.Area];%获得所有4连通区域的前景区域的面积向量

bigRegion = find(dataArea == max(dataArea));%获得最大面积的区域的标识号

finalMask = (labeled == bigRegion);%原始finalMask中具有最大面积的区域，即皮损区域

finalMask = bwmorph(finalMask,'remove');%获得皮损区域的边界像素

对我们的色素性皮损图像而言，采用上述步骤去除小的干扰区域所得结果如图2(c)所示，最终提取的皮损区域边界参考图2(d)。
   5.叠加边界至原始彩图

为了便于查看上述过程提取的边界和真实的皮损区域边缘的吻合程度，可以通过下面的代码来实现叠加边界至原始彩图的功能。

rChannel = imgBak(:,:,1); 

%获取原RGB图像中R通道数据，在代码段2中用imgBak备份原始彩色图像的数据

gChannel = imgBak(:,:,2); %获取原RGB图像中G通道数据

bChannel = imgBak(:,:,3); %获取原RGB图像中B通道数据

%下面三行代码的功能：使用逻辑索引的方式方便地访问对应于皮损区域边界位置的原始图像数据，在原始图像

%中将皮损区域边界所对应的像素点设置为纯蓝色rgb = （0，0，255）

rChannel(finalMask) = 0; 

gChannel(finalMask) = 0;

bChannel(finalMask) = 255;

imgBak = cat(3,rChannel,gChannel,bChannel);%重构得到带蓝色边界标识的rgb彩色图像

figure,imshow(imgBak),title('Original Image with Mask');%显示皮损区域边界和原始RGB图像的彩色叠加效果

%保存彩色叠加结果图像

location = findstr(fileName,'.jpg');%定位’.jpg’字符串在fileName中的起始位置

fileName = fileName(1:location - 1);%和上一步一起获得打开文件所在的路径

fileName = strcat(fileName,'Result.jpg');%用于保持彩色叠加结果图像的文件名

imwrite(imgBak,fileName);%保存彩色叠加结果图像，第一个参数表示要保存的图像数据文件，第二个参数表示存%储的文件名

在上述代码段中，有一个小的编程技巧值得一提，就是逻辑索引（Logical Index）。如果一个逻辑数组mask的大小等于数组A（即mask的行列数分别等于A的行列数），那么A（mask）= value意味着：如果mask(i,j)为true（即1）,那么A（i,j）= value。尽管我们也可以通过二层for循环间套来实现同样的功能，但是显然没有通过逻辑索引的方式来地简练。

将分割所得皮损区域的边界叠加到RGB原始图像后所得结果可以参考图3，从分割结果看，使用上面的方法能较好地对色素性皮损区域进行定位。

图3一些色素性皮损图像的自动分割结果

三、结语

Matlab内含的大量工具箱函数以及内部函数在性能方面往往都做了优化，从而为科学工作者以及工程开发人员提供了一个良好的平台。在熟悉工具箱及内部函数的功能后，可以进行系统有机的组合以快速地验证我们对问题的解决思路是否合理。如果经过验证发现我们的思路是正确的，那么我们就可以根据项目开发的需求采用其他编程语言进行系统开发，或者将其他语言和Matlab进行混合编程。因此，从这个角度看，Matlab在降低开发人员的时间和精力的消耗，以及推进项目进程方面具有重要的意义。