1 引言

近年来, 随着科学技术的高速发展, 颜色作为与物质化学信息紧密联系且能够进行客观测量的光学属性, 越来越受到各国研究者的关注, 从而引发了计算机视觉在分析化学中的广泛应用. 当前, 最典型的例子就是Suslick等人基于比色传感器阵列所构建的可视化仿生鼻技术^[1]. 作为可视化仿生鼻的核心关键, 比色传感器阵列(如图1所示)由多个具有非特异性的化学传感器单元, 以5×5或6×6的等阵列形式组成. 这些传感器单元与不同目标的分析物发生分子间相互作用时会导致吸收光谱发生改变, 并呈现出具有特异性的颜色的改变^[2]. 可视化仿生鼻技术是借助扫描仪、摄像头等设备提取传感器单元反应前后所产生颜色差异特征, 从而实现物质定性甚至定量的分析^[3]的先进技术. 目前该技术已经能够实现对挥发性气体、生物样品、离子和有机小分子等多种微量级物质的检测^[3,4], 具有了响应速度快、高选择性和高特异性等优点.

在基于比色传感器阵列的可视化仿生鼻技术中, 传感器单元的有效分割^[4,5]决定了特征提取信息的可靠性, 它是后期计算机视觉识别的前提和关键. 当前比色传感器阵列在特征提取中主要以人工手动分割方式为主^[5], 这种方式效率低下、可重复性较差, 且容易引入人为误差等. 为了克服上述不足, 基于图像阈值^[6,7]、区域生长^[8]、人工神经网络^[9]、聚类分析^[10]的分割方法相继被提出并应用于阵列图像分割, 取得了不错的效果. 这些方法中, 基于图像阈值的Otsu算法是最常采用的方法^[5], 该算法便捷简单且应用广泛, 但其在主要以灰度阈值分割为主, 考虑因素单一, 在分割中很容易受到图像噪声以及光照环境所影响而产生错误的分割结果. 区域生长的图像分割算法在比色传感器阵列分割中也展示了不错的整体性能, 但该方法很容易引入错误的边界, 对于高维且富含不同特征的彩色图像, 该算法的时间复杂度是需要克服的难题. 人工神经网络的方法可以实现具有复杂信息图像的有效分割, 但其网络层的设计缺乏理论基础, 同时在分割前需要时间进行大样本的训练和学习^[11], 并且它也离不开算法复杂度的限值. 对于聚类分析的分割算法来说, 它的基本原理是将图像中特征空间内相似的像素点归为一类, 特别是基于特征空间聚类的模糊C均值(此后简称FCM)图像分割算法, 其不直接表明某像素点是否属于某一类, 而是通过值范围为[0, 1]的模糊隶属度来描述该点属于某一类的程度, 最后通过不断迭代隶属度和聚类中心使其目标函数值达到最小的情况下, 实现最优的图像分割^[12,13]. 因该算法在图像分割领域应用广泛, 其存在的主要难点是如何有效地确定初始聚类条件.

由于FCM算法所提出的模糊隶属度理论与图像信息模糊性的特点非常匹配^[14], 且考虑到阵列图像为颜色信息各异的二维点阵图像, 借助FCM算法的目标函数能够在分割中利用到更多的图像信息. 因此, 本文提出了一种基于模糊C均值聚类的比色传感器阵列图像分割算法. 首先通过网格划分下直方图信息分析实现初始聚类条件的确定, 然后通过在目标函数引入色彩信息实现了比色传感器阵列的精确分割. 通过不同类型阵列点图像分割测试, 该算法相比于Otsu算法、以及传统的FCM算法具有更优的整体分割精度, 为比色传感器阵列特征提取提供了一种行之有效的图像分割方法.

2 HSI颜色空间

颜色空间不仅影响着图像特征提取的有效性, 同时也决定了颜色信息的量化基础以及后期分析结果的可靠性等问题^[15,16]. 此外, 比色传感器阵列图像包含了大量的颜色信息, 属于彩色图像分割的范畴. 因此, 选择合适的颜色空间是必须首要考虑的问题.

当前分析化学广泛应用的颜色空间中, HSI颜色空间巧妙地将色彩信息和亮度分离开来, 各分量间具有高度的相对独立性, 而且I分量可以看作RGB颜色空间下的均值灰度化结果, 富含了大量的图像信息, 仅仅利用I分量就可以获得较好地分割效果^[17,18]. 除此之外, 表征最主要颜色信息的H分量基本不受光照环境影响. 因此, 图像分割中利用H分量可以有效避免提光照环境的影响. 综合考虑, 该研究选择了在HSI颜色空间下实现比色传感器阵列图像分割算法设计, 而实现后续图像分割的首要任务则是将采集到的RGB图像根据式(1)转换到HSI颜色空间.

$\left\{ \begin{array}{l} S = 1 - \dfrac{{3 \cdot \min (R,G,B)}}{{R + G + B}}\\ H = \left\{ \begin{array}{l} \theta \begin{array}{*{20}{c}} {}&{} \end{array}\begin{array}{*{20}{c}} {}&{if\begin{array}{*{20}{c}} {B \le G}&{} \end{array}} \end{array}\\ 2\pi - \theta \begin{array}{*{20}{c}} {}&{if\begin{array}{*{20}{c}} {B > G}&{} \end{array}} \end{array} \end{array} \right.\\ \theta = {\cos ^{ - 1}}\left( {\dfrac{{(R - G) + (R - B)}}{{2\sqrt {{{(R - G)}^2} + (R - B)(G - B)} }}} \right)\\ I = \dfrac{{R + G + B}}{3} \end{array} \right.$

(1)

3 算法设计 3.1 网格化分下的聚类条件初始化

首先根据Deng等人^[19]在生物芯片类等微阵列图像中所提出的图像投影方法基础上. 如图2所示, 将阵列图像的I分量, 根据式(2)实现水平方向和垂直方向投影, 并得到了对应行、列方向的投影曲线.

$ \left\{ \begin{gathered} pro\_row(i) = \sum\limits_{j = 1}^N {\left( {{\text{1 - }}I(i,j)} \right)} \begin{array}{*{20}{c}} {}&{i = 1,2,\cdots,N} \end{array} \\ pro\_col(j) = \sum\limits_{i = 1}^M {\left( {{\text{1}} - I(i,j)} \right)} \begin{array}{*{20}{c}} {}&{j = 1,2,\cdots,M} \end{array} \\ \end{gathered} \right. $

(2)

其中, I(i, j)为对应像素位置的I分量值, 经归一化处理后其值变化范围为[0, 1], M为行数, N为列数.

如图2(b)、图2(c)所示的行、列投影曲线, 其中红色菱形所标记的“波峰”对应着阵列图像的传感器单元区域, “波谷”则为对应阵列图像的背景区域. 通过计算相邻波峰坐标的中值即获得对应的网格划分坐标, 并通过该坐标将阵列图像划分成如图所示的多个子区域.

3.1.1 聚类数目的确定

通过网格划分, 图像分割范围从全局分割缩小为局部的阵列点图像分割(如图3所示), 而阵列点图像中只存在需要实现特征提取的目标和背景两部分, 因此聚类数目可以很容易的确定为2.

3.1.2 初始聚类中心的确定

对于单通道图像的直方图, 特征一致性区域会以显著性波峰存在, 该波峰位置所对应的分量值可以作为FCM图像分割算法的初始聚类中心^[20]. 以I分量为例进行直方图信息统计, 根据式(3)进行平滑操作以避免噪声尖峰的干扰, 最终得到了如图4(d)所示的平滑直方图.

$ {H_{smooth}}(i) = \dfrac{{\sum\limits_{n = - N}^N {H(i + n)} }}{{2N + 1}}\begin{array}{*{20}{c}} ,&{N \leqslant i \leqslant L - N} \end{array} - 1 $

(3)

其中, N为平滑窗的半窗尺寸, N值一般设置为2到5之间的整数最为合适; L为图像颜色分量的灰度级范围, 在研究中统一规定其值变化范围在[0, 255]之间.

结合文献[20]所提出的波峰筛选算法, 并根据多次试验, 设置峰值高度筛选阈值T_peak=0.006×S(S为该图像的总像素数目), 峰峰最小间距T_distance=30, 可以得到如图4(d)所示的两个显著性波峰P₁和P₂, 则这两个波峰所对应的分量值I₁, I₂即为初始聚类中心. 在HSI图像中, H分量存在不稳定和易受噪声影响等特点, 它的平滑直方图经过波峰筛选后可能会存在3个或者3个以上的波峰, 则可以选择峰值最大的两个波峰对应的分量值作为初始聚类中心.

3.2 基于加权HI分量的FCM图像分割算法

假设将一个包含n个样本的数据集X={x₁, x₂, …, x_n}划分为c个子类(2≤c≤n), 其中x_i表示数据集中的某一个样本(1≤i≤n), FCM算法实质就是在式(4)的条件限制下, 通过不断迭代使式(5)中的目标函数值最优化的一种无监督分类方法.

$\sum\limits_{k = 1}^c {{u_{ik}}} = 1,0 \leqslant {u_{ik}} \leqslant 1$

(4)

${J_{\rm{FCM}}} = \sum\limits_{i = 1}^n {\sum\limits_{k = 1}^c {u_{ik}^m{{\left\| {{x_i} - {v_k}} \right\|}^2}} } $

(5)

式(4)、(5)中, v_k表示c个子类中第k个子类的聚类中心; u_ik为数据集中第i个样本对第k个子类的模糊隶属度; m则表征了FCM算法的模糊程度, 其值越小, FCM算法的性质接近硬聚类算法, 在本文中m取值为2; $\left\| \bullet \right\|$ 表示某一样本到其聚类中心的范数, 当前的FCM算法主流研究中多采用2维范数, 即欧式距离, 本文也不例外.

结合传统FCM算法的目标函数, 重新定义其目标函数为:

${J_{{\rm{HI\_FCM}}}} = \sum\limits_{i = 1}^n {\sum\limits_{k = 1}^c {u_{ik}^m\left( {{W_H} \cdot {D^2}\left( {{H_i},{v_{Hk}}} \right) + {W_I} \cdot {D^2}\left( {{I_i},{v_{Ik}}} \right)} \right)} } $

(6)

其中, W_H、W_I分别为H分量和I分量的权重系数, H_i、I_i分别为H分量和I分量中的第i个待分类的样本值, v_Hk、v_Ik分别为H分量和I分量的第k个聚类中心值, D(H_i, v_Hk)表示H分量中第i个待分类样本与聚类中心v_Hk的相似度距离, 由于HSI颜色空间下, 图像的H分量具有圆循环特性^[11,21], 其计算公式为:

$D\left( {{H_i},{v_{Hk}}} \right) = \left( {\frac{{{H_{\max }}}}{2} - \left| {\left| {{H_i} - {v_{Hk}}} \right| - \frac{{{H_{\max }}}}{2}} \right|} \right)$

(7)

式(7)中, H_max为H分量的最大值, H分量为8位数据格式时, 其最大值为255. D(I_i, v_Ik)则为I分量中第i个待分类样本与聚类中心v_Ik的相似度距离, 计算公式为:

$D\left( {{I_i},{v_{Ik}}} \right) = \left| {{I_i} - {v_{Ik}}} \right|$

(8)

结合式(7)、(8), 运用拉格朗日乘子法求解式(6), 可以得其模糊隶属度矩阵U=[u_ik]:

${u_{ik}} = \frac{{{{\left( {{W_H} \cdot {D^2}\left( {{H_i},{v_{Hk}}} \right) + {W_I} \cdot {D^2}\left( {I{}_i,{v_{Ik}}} \right)} \right)}^{\frac{{ - 1}}{{m - 1}}}}}}{{\displaystyle\sum\limits_{j = 1}^c {{{\left( {{W_H} \cdot {D^2}\left( {{H_i},{v_{Hj}}} \right) + {W_I} \cdot {D^2}\left( {{I_i},{v_{Ij}}} \right)} \right)}^{\frac{{ - 1}}{{m - 1}}}}} }}$

(9)

因此, H分量和I分量聚类中心向量V={v_Hk, v_Ik}也可以根据式(10)而确定.

${v_{Hk}} = \dfrac{{\displaystyle\sum\limits_{i = 1}^n {u_{ik}^m{H_i}} }}{{\displaystyle\sum\limits_{i = 1}^n {u_{ik}^m} }}\begin{array}{*{20}{c}} ,&{{v_{Ik}} = \dfrac{{\displaystyle\sum\limits_{i = 1}^n {u_{ik}^m{I_i}} }}{{\displaystyle\sum\limits_{i = 1}^n {u_{ik}^m} }}} \end{array}$

(10)

目标函数(6)重新代入式(9)所计算的最新模糊隶属度矩阵U和聚类中心向量V, 不停地迭代直到满足式(11)的迭代终止条件, 则可以实现HSI图像的最优聚类, 并最终实现图像分割.

$\left| {U(t + 1) - U(t)} \right| < \varepsilon $

(11)

3.3 权重系数W_H和W_I的确定

为了在图像分割中能够引入色彩信息, 研究所设计的算法在目标函数中同时引入了H分量和I分量, 并加入了W_H和W_I两个权重系数, 其值大小也暗示了对应分量在图像分割中的重要程度.

3.3.1 峰值熵

熵在信息论中具有非常重要的地位, 它是对系统本身的不确定性的一种度量方法. 由于HSI图像中H分量和I分量均来自同一图像, 它们各自的图像熵差异不大. 因此, 为了表征各分量的噪声程度和不稳定性, 因此引入了峰值熵这一概念. 峰值熵, 即仅仅统计在直方图中所有符合条件峰值的熵值, 结合式(12)可分别得到H分量的峰值熵H_peak_E和I分量的峰值熵I_peak_E的计算公式:

$\left\{ \begin{gathered} {H_{\rm{peak}}}\_E = - \sum\limits_{i = 1}^{{m_H}} {{p_{Hi}}\log ({p_{Hi}})} \\ {I_{\rm{peak}}}\_E = - \sum\limits_{i = 1}^{{m_I}} {{p_{Ii}}\log ({p_{Ii}})} \\ \end{gathered} \right.$

(12)

式(12)中, m_H、m_I分别为经过直方图峰值算法筛选后的峰值数目, p_Hi和p_Ii分别为H分量、I分量中第i个峰值的均一化频率值.

3.3.1 图像对比度

图像对比度是图像质量评价的重要参数, 它是对图像中最亮的白和最暗的黑之间不同亮度层级的测量, 反映了图像的清晰度和纹理沟纹深浅的程度. 图像对比度越大, 图像越清晰, 纹理沟纹越深, 同时也有可能暗示其噪声干扰比较严重, 反之其值越小, 图像越模糊, 纹理沟纹越浅, 图像噪声干扰较小. 因此, 对于阵列点图像H分量的对比度H_C和I分量的对比度I_C的计算公式分别为:

$\left\{ \begin{array}{l} H\_C = \sum\limits_{i = 1}^k {\sum\limits_{j = 1}^k {{{\left( {\dfrac{{{H_{\max }}}}{2} - \left| {\left| {i - j} \right| - \dfrac{{{H_{\max }}}}{2}} \right|} \right)}^2}} } {P_H}(i,j)\\ I\_C = \sum\limits_{i = 1}^k {\sum\limits_{j = 1}^k {{{(i - j)}^2}} } {P_I}(i,j) \end{array} \right.$

(13)

式(13)中, k为H分量、I分量进行量化后的灰度级数, 为了降低整体算法的时间复杂度, k值设置为8; P_H(i, j)和P_I(i, j)为H分量、I分量的灰度共生矩阵; H_max为量化后H分量的最大灰度值, 由于k为8, 因此H_max设置为7.

3.3.3 权重系数确定

由于在高I值和低I值时, H值非常不稳定, 而I_C受图像的整体亮度所影响, 并且一定程度上反映了图像的光照情况, 因此H_peak_E和I_C组成了影响W_H的主要因素. 另一方面, 相对于I_C来说, H_C几乎都处于高值状态, 它可以影响 I分量在分割中的重要程度, 因此I_peak_E和H_C则组成了影响W_I的主要因素. 综上所述, W_H和W_I可以通过经验公式(14)来确定.

$\left\{ \begin{gathered} {W_H} = \dfrac{{{H_{\rm{peak}}}\_{E^2} \cdot I\_{C^2}}}{{{H_{\rm{peak}}}\_{E^2} \cdot I\_{C^2} + {I_{\rm{peak}}}\_{E^2} \cdot H\_{C^2}}} \\ {W_I} = \dfrac{{{I_{\rm{peak}}}\_{E^2} \cdot H\_{C^2}}}{{{H_{\rm{peak}}}\_{E^2} \cdot I\_{C^2} + {I_{\rm{peak}}}\_{E^2} \cdot H\_{C^2}}} \\ \end{gathered} \right.$

(14)

3.4 算法流程图

综上所述, 本文基于传统FCM算法, 结合直方图信息和色彩信息提出了一种加权HI分量的比色传感器阵列图像分割算法, 它最终的算法流程图如图5所示.

4 图像分割效果测试

通过1024×768分辨率摄像头(奥尼酷克生产)采集了5张6×6阵列尺寸的24位比色传感器阵列真彩图像, 经过网格划分后得到共计180张阵列点图像(分辨率为60×60). 这180张阵列点图像根据光照、噪声情况等分为三类: 正常图像127张, 光照不均匀或者模糊的图像36张, 带有污染斑点的图像17张.

为了测试所提出的图像分割算法的的性能, 针对这180张阵列点图像进行了图像分割效果测试. 在测试中, 除了所提出的加权HI分量的FCM图像分割算法(简称为H&I_FCM算法), 还使用了阵列图像在RGB颜色空间下常用的Otsu阈值分割算法(简称RGB_Otsu算法), 以及只包含I分量的传统FCM算法(简称I_FCM算法)等. 所有的图像分割测试都是基于1.7 GHz CPU频率, 4 GB内存计算机, 64位Windows 8操作系统下的Matlab 2014a平台实现的.

4.1 分割精度评价

首先通过比色传感器阵列专业研究人员对这180张60×60阵列点图像进行了手工分割, 其处理结果作为标准图像, 然后分别使用H&I_FCM、I_FCM、RGB_Otsu这三种算法实现了所有阵列点图像的分割, 并根据不同类型图像统计各自的平均分割精度SA^[12,22]来评价分割结果, 其中SA的计算公式如下:

$SA = \frac{{\left| {{O_S} \cap O} \right| + \left| {{B_S} \cap B} \right|}}{{\left| {{O_S}} \right| + \left| {{B_S}} \right|}}$

(15)

式(15)中, O_S和B_S分别为标准图像的目标和背景, O和B为测试算法分割结果的目标和背景, $\left| {{O_S}} \right| + \left| {{B_S}} \right|$ 为图像的总像素数目, $\left| {{O_S} \cap O} \right|$ 为正确分割的目标像素数目, $\left| {{B_S} \cap B} \right|$ 为正确分割的背景像素数目.

4.1.1 不同算法分割精度评价

如图6所示为三种不同类型阵列点图像在不同分割算法下的典型分割结果, 表1为三种算法在所有阵列点图像分割中的分割精度统计结果.

通过该结果可以得知: 在正常阵列点图像分割中, 三种算法的分割精度不相上下; 在光照不均或模糊的阵列点图像分割中, 由于H&I_FCM考虑了色调信息, 其90.75%的分割精度在三者中最高, I_FCM和RGB_Otsu分割精度比较接近; 在含有污染斑点的阵列点图像分割中, I_FCM的分割精度最高96.52%, 但H&I_FCM分割精度也高达96.30%, RGB_Otsu的分割精度最低. 总体来说, 相比于其余两种图像分割算法, H&I_FCM在所有阵列点图像分割中展示了96.54%的总体最优分割精度, 特别是在光照环境和污染斑点等因素影响下, 仍然保持了良好的分割效果.

4.1.2 不同权重系统分割分割精度评价

为了评价H分量的引入对于分割精度的影响, 统计了H&I_FCM算法下, 三种不同类型阵列点在不同权重系数H分量下的分割精度. H分量的权重系数取值范围为0–1, 取值间隔为0.1.

如图7所示, 三种不同类型阵列点图像分割中, 随着H分量在目标函数中的权重增加, 在初期分割精度有所提高, 但是在后期随着H分量比重的不断增大, 分割精度会下降, 甚至低于单纯I量的FCM算法分割精度. 由于H分量的不连续性, 其在图像分割中所占的比重不宜过大, 因此需要通过合理的权重系数公式对H分量所占比重进行分配. 式(14)通过借鉴Rajaby^[11]所提出的权重分配公式, 通过平方和方式进一步增大了I分量所占的比重, 将W_H控制在合理的水平范围内.

4.2 运行时间评价

时间复杂度是评价算法好坏的重要指标, 对于一种优良的图像分割算法来说, 除了具备优良的分割结果, 其运性时间还必须控制在可以接受水平范围内. 因此, 本文在图像分割测试中分别统计了三种算法在5张阵列图像以及网格划分后的180张阵列点图像的平均运行时间 $\overline t $ , 其统计结果如表2所示.

通过表2所示平均运行时间可以得知: 在比色传感器阵列图像分割中, H&I_FCM的运行时间最长, I_FCM的运行时间最短, RGB_Otsu适中. 实际上, I_FCM和RGB_Otsu都是单通道的图像分割, 但是I_FCM的运行时间明显优于RGB_Otsu. 这也表明单通道比色传感器阵列图像分割中, 传统FCM图像分割算法在效率上是优于RGB_Otsu算法的. H&I_FCM的运行时间由于在分割中使用双通道图像信息, 并且涉及W_H和W_I的计算, 因此它的运行时间在三者之中最为耗时是合乎情理, 但其运行时间处于可以接受的水平范围.

5 结论

在比色传感器阵列图像分割中, 传感器单元的颜色信息是实现后续特征分析的基础, 而优良的分割精度保证了信息提取的正确性, 因此分割精度应该是本文算法设计中更看重的方面. 结合本文的图像分割效果测试结果, 本文可以得出以下结论: 相比于传统的FCM算法和当前常用的RGB_Otsu算法, 本文所设计的H&I_FCM图像分割算法在所有阵列点图像分割中展示了96.54%的总体最优分割精度, 更能够适应不同光照环境和斑点污染的影响, 可以有效、准确地实现比色传感器阵列图像的分割. 当然本文所设计的算法在运行时间上还有待改进, 该内容也是后续的研究重点.