CN103236065A - 基于主动轮廓模型和细胞神经网络的生物芯片分析方法 - Google Patents

Abstract

本发明公开了一种基于主动轮廓模型和细胞神经网络的生物芯片分析方法，包括以下步骤：对矩形样点采用改进的Hough变换进行倾斜校正，而对于圆形样点则采用改进的Radon变换；基于投影法对样点进行初次定位，并生成优化后的网格；然后基于邻域搜索自适应地调整网格，对样点进行二次精确定位；根据贪婪算法优化主动轮廓模型，利用CNN将样点按信号强弱进行分类；将Multiple snakes与CNN结合，CNN先学习强信号样点snake的收敛行为，再指导弱信号样点snake的收敛，最终实现样点的合理分割；提取并输出微阵列样点信号数据。本发明解决了生物芯片图像倾斜校正、形状不规则样点及弱信号样点分割困难等问题，实现了生物芯片样点的自动识别，适合大规模的生物芯片样点快速分析。

Description

基于主动轮廓模型和细胞神经网络的生物芯片分析方法

技术领域

本发明属于生物芯片技术领域，特别是一种基于主动轮廓模型和细胞神经网络的生物芯片分析方法，该方法可应用于生物芯片扫描图像的处理、自动识别和分析。

背景技术

生物芯片是通过平面微细加工技术在固体芯片表面构建的微流体分析单元和系统，以实现对细胞、蛋白质、核酸以及其他生物组分的准确、快速、大信息量的检测，具有多通道、高通量、并行自动处理等优点。生物芯片技术充分利用了生物科学、信息学等学科的成果，已经广泛应用于医学、生命科学、环境科学等相关领域。

生物芯片样点识别是生物芯片检测系统的关键技术。生物芯片上集成了大量点阵的生物识别分子，每个样点对应一个基因表达或者生物信息，对生物芯片上所包含的信息进行准确检测是一项至关重要的工作。生物芯片信息处理的关键是准确识别每个样点，在样点区域内提取灰度强度作为后续处理的依据。样点识别是生物信息分析的基础，识别准确与否直接影响后续数据提取及生物信息建立，其速度和精度对整个生物芯片信息分析具有决定性作用。

样点识别技术是生物信息提取和分析方法领域的一个经典难题。芯片图像经常受到点样仪误差、扫描仪误差、反应过程污染、清洗过程污染以及环境光污染等因素的综合影响，从而增大了准确、快速地识别样点的难度。而且现在还没有统一的国际标准或行业标准，算法的一致性、精确性、图像依赖性等关键技术指标不明确，用户难以选择合适芯片特点的样点识别方案。目前较多的生物芯片分析软件仍采用人为干预相结合的方法，一般做法是：首先手动输入设置网格参数，自动生成粗对准网格，然后手动调整实现目标样点局部区域化，按照一定规则不断调整分割区域的位置和大小，最后实现目标样点的分割和数据提取。对于生物芯片高密度、大信息量的特点，手动辅助识别的方式显得效率低下，对操作人员的主观依赖性大，精度难以提高。

KASS于1987年提出SNAKE模型，现已广泛应用于数字图像分析和计算机视觉领域，这种主动轮廓模型具有良好的提取和跟踪特定区域内目标轮廓的能力，因此非常适合边缘提取。LEON O.CHUA等于1988年提出了细胞神经网络，这种网络在分类、学习和预测问题上有较好性能，并提供了分布式信息存储和并行处理能力，使它在生物芯片判定方面具有得天独厚的优势。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种基于主动轮廓模型和细胞神经网络的生物芯片分析方法，解决样点信息提取和分析上的难题，准确、快速和自动地识别样点。

为解决上述技术问题，本发明所采用的技术方案是：一种基于主动轮廓模型和细胞神经网络的生物芯片分析方法，该方法为：

1)若原始生物芯片图像未发生倾斜，即图像准直，则进入4）；若图像发生倾斜，即倾斜角度在[-5°,5°]范围内，则进入2）；若图像倾角超出[-5°,5°]范围，则认为芯片图像的质量不符合要求，不进行后续的分析；

2)将生物芯片图像二值化，计算二值的生物芯片图像的倾斜角度，得到其相反数为最佳的倾斜校正角度；

3)利用所述最佳的倾斜校正角度和图像几何空间变换的原理，对倾斜的原始生物芯片图像进行旋转校正并灰度化，得到准直的芯片灰度图像；

4)对所述准直的芯片灰度图像进行投影计算，得到准直图像在水平和竖直两个方向上的投影波形P_H(x),P_V(y)，将投影波形进行开闭级联运算操作和自适应二值化处理，得到二值投影方波；

5)二值投影方波P_H(x)的波峰中心对应着每一行样点的中心，在所述中心划分网格线，确定样点行在水平方向的坐标；二值投影方波P_H(x)的波谷中心对应着样点行间空隙的中心位置，在所述中心位置划分网格线，将每个样点分割到各自的子区域中；运用上述同样的方法处理二值投影方波P_V(y)；

6)读取上述二值投影方波的各个波峰和波谷的宽度，以及波峰中心与波谷中心的位置；运用数学统计学原理，计算出贴近真实的平均波峰宽度和平均波谷宽度，以及平均波峰中心距；

7)结合数学统计学原理，通过判断每条网格线对应的方波的波峰宽度和中心距，发现冗余网格线并予以剔除；对剩余的网格线重新进行计算，添补缺失的网格线；对网格线的位置与数量进行调整，生成优化后的网格线，完成对样点的初次定位；

8)将步骤6）获得的平均波峰宽度用作二次精确定位网格的初始直径，对每个样点生成一个二次精确定位初始网格；对于圆形样点的芯片图像，网格设计成圆形；对于矩形样点的芯片图像，网格设计成矩形；

9)基于邻域搜索和自适应调整的原理，对样点进行二次精确定位：运用邻域搜索算法，对上述产生的二次精确定位网格的中心位置逐个进行调整；通过比较网格内部样点信号数据，确定尺寸调整方向，再按照直径调整准则进行调整；使得网格最终能将样点套住，信号区域完全落在网格内部，同时套住的背景区域像素数目达到最小值；

10)根据贪婪算法优化主动轮廓模型，利用细胞神经网络CNN将芯片样点按信号强弱进行分类；将多主动轮廓模型Multiple snakes与CNN结合，CNN先学习强信号芯片样点snake的收敛行为，并得到强信号样点的snake轮廓；再指导弱信号芯片样点snake的收敛，得到弱信号芯片样点的snake轮廓，最终实现芯片样点信号区域的分割，提取芯片样点信号区域内灰度强度作为后续处理的数据；

11)计算每个样点信号区域内部的灰度总值、平均值和方差，统计得到芯片样点的信号值；计算样点snake轮廓外部环形区域的灰度强度信息，统计出每个样点的背景值，最后输出生物芯片样点分析数据。

以下进一步详细阐述本发明的原理。

（1）基于改进Hough变换和改进Radon变换的倾斜快速校正算法

本发明的目的之一是提供一种基于Hough变换和Radon变换的倾斜快速校正算法，综合考虑了样点形状、校正精确度和运算速度等因素，解决了现有技术中生物芯片图像倾斜检测精度和速度不佳的问题。

为了达到上述目的一，经过比较常见的主流倾斜校正算法，本发明提出对矩形样点采用改进的Hough变换进行倾斜校正，而对于圆形样点则采用改进的Radon变换，该方法解决问题的具体技术方案由如下A、B和C三部分组成。

（A）针对矩形样点改进的Hough变换方法，包括如下步骤：

步骤A1：对生物芯片图像进行二值化

将倾斜角度范围规定为[-5°,5°]，角度检测精度为0.01°。大量统计得到生物芯片图像的倾斜角度范围是[-5°,5°]，如果图像的倾角超过这范围，则可以认为图像质量不符合要求。但是生物芯片要求的倾角检测精度比一般的直线倾角检测要高，故将倾角检测精度规定为0.01°。

在检测倾斜角度之前，先对生物芯片图像进行灰度化，然后进行一次形态学的开闭运算级联操作。接着选择Otsu算法自动计算最佳分割阈值，并使用这阈值将芯片图像二值化。所得二值图像的白色区域（灰度值为255）代表芯片样点，黑色区域（灰度值为0）代表芯片背景。

步骤A2：计算Hough参数空间

在进行Hough变换之前先对像素点进行判别，如果当前像素点符合条件，则进行Hough变换，否则跳转到下一个像素点。

像素点的约束条件：如果当前像素的灰度值为0，且下一行同样位置像素的灰度值为255；或者如果当前像素灰度值为255，且下一行同样位置的像素灰度值为0，则计算当前像素点的Hough参数(ρ,θ)，即仅统计矩形样点的上下边缘像素。其中，ρ为直线到计算原点的距离，θ为直线与轴之间的夹角。因此，这种改进的Hough算法仅适合矩形样点，不适合圆形样点。

以二值图像的中心为计算原点，以对角线长度为ρ轴最大值，步长为1；以5°为θ轴的最大值，步长为0.01°，计算Hough参数空间。统计ρθ参数空间中各曲线交点的信号累加值，对应着各个角度上特征点的累计数量。

步骤A3:提取占主导的直线组

无论是准直图像还是倾斜图像，同一行矩形样点的上（下）边缘都处在同一条直线上，故可以通过累计同一行样点的上（下）边缘长度，来判别该行是否对图像倾斜贡献较大。显然，占主导地位的都是累计长度较长的直线，对应的矩形样点行信号强度好、样点数量多。

在提取占主导地位的样点上（下）边缘之前，选择图像宽度的1/10为最小直线长度A_min。遍历Hough参数空间，如果参数空间坐标(ρ_i,θ_j)的累加值为A(i,j)，代表直线xcosθ_j+ysinθ_j=ρ_i的累计长度为A(i,j)。如果线长度A(i,j)小于A_min，则表明该直线不占主导地位，应舍去不记录。反之，如果线长度A(i,j)不少于A_min，则要进一步判别它是否局部极大值。根据检测精度确定合理的局部峰值半径r，将Hough参数空间中A(i,j)的r×r邻域值与当前累加值A(i,j)比较，如果A(i,j)在局部区域内是极大值，则表明该直线起主导作用，应保存并用于下一步计算。

步骤A4:计算直线组的平均倾斜角度

将提取到的直线组按照累计长度从大到小地排序，从而得到最大累计长度A_max，并分别计算累计长度较长的前10条直线的相对长度

（累计长度最长的直线的相对长度为1）。最后，提取相对长度不少于0.5的直线，计算它们的平均倾斜角度

则最佳的倾斜校正角度为平均倾斜角度的相反数

其计算公式为：

$\stackrel{\‾}{\mathrm{\θ}}=\frac{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\θ}}_m{R}_m}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{R}_m}$

其中，相对长度不少于0.5的直线有M条，R_m为直线的相对长度，θ_m为直线的倾斜角度。

由Hough变换的原理可知，倾角为

的一族直线占主导地位，对应ρθ参数空间中的累计值也是局部极大值。观察Hough正弦图，在

附近正弦曲线呈有规律的聚焦状，图中局部极亮点代表矩形样点上下边缘的所在位置。

（B）针对圆形样点改进的Radon变换方法，包括如下步骤：

步骤B1:提取圆形样点边缘

在检测倾斜角度之前，先对圆形样点图像进行灰度化和二值化，所用方法跟改进的Hough变换中所用方法是一样的。

由于Radon变换的计算量较大，为了减少检测的计算量，在进行Radon变换之前，先提取圆形样点的边缘。使用最小的结构元B对原二值图像A进行腐蚀，再用原二值图像A减去腐蚀图像ε_B(A)，即可得到边缘图像β(A)=A-ε_B(A)，且其边缘宽度为1。

对比提取边缘前后的Radon参数空间可知，提取边缘后得到的正弦图的整体亮度减弱，即正弦曲线的数量减少，但具体位置未变，故仍能很好地代表样点投影信号。通过精简正弦图来减少计算量，保持了准确度而又不影响检测效果，这就是Radon变换的初步改进。

步骤B2：二级Radon变换

Radon变换一般将图像的中心位置作为投影原点，采用其离散形式，投影信号表达公式为：

$g(\mathrm{\ρ},\mathrm{\θ})=\underset{x=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}f(x,y)\mathrm{\δ}(x\cos \mathrm{\θ}+y\sin \mathrm{\θ}-\mathrm{\ρ})$

其中，(x,y)为二值边缘图像中的坐标，(ρ,θ)为Radon参数空间中的坐标，ρ为像素点到投影原点的距离，θ为投影角度，x,y,ρ,θ皆为离散变量。固定θ值为θ_k，通过覆盖图像所要求的ρ的所有值产生一个投影g(ρ,θ_k)。改变θ值并重复上述过程，则产生另外的投影。最终得到以ρ,θ为变量的二维数组g(ρ,θ)，并可以利用这数组合成正弦图，而正弦图中像素的亮点与对应的g值成正比。

针对圆形样点边缘进行Radon变换，如果要求的检测精度较高，计算量的数量级也相应增加。为了减少搜索最佳的倾斜校正角度的次数，采用二级不同尺度的搜索策略，先进行粗略搜索，在确定大致范围后再进行精细搜索。

在倾斜角度[-5°,5°]区间内θ以0.1°为步长，对二值边缘图像进行第一级Radon变换，得到宽度为101的Radon参数空间g(ρ,θ)。固定θ值为θ_k，计算投影g(ρ,θ_k)的累计值

遍历所有累计值，最大累计值A_m所对应的θ_m为初次倾斜角度。记初次倾斜角度为

在

区间θ以0.01°为步长，进行第二级Radon变换，得到宽度为301的更精细的Radon参数空间g(ρ,θ)。同样原理，通过计算投影的累计值并比较出最大值，从而得到对应的二次倾斜角度

二次倾斜角度的相反数

即为精确度更高的最佳的倾斜校正角度。由Radon变换的原理可知，边缘图像在的方向上的投影为g(ρ,θ_k)，其累加值A_k是局部极大值。在

附近正弦曲线有规律地成束分布，Radon正弦图中正弦曲线束代表圆形样点边缘的所在位置。

不同尺度的搜索策略在很多倾斜校正的算法中都有体验，其减少存储空间、提高运算速度的作用是显而易见的。大量实验证明，采用二级Radon变换算法后计算量明显减少，要比同样精度下直接Radon变换的计算量少一半以上。

（C）倾斜图像的旋转校正

在获得最佳的校正角度后，可以利用图像几何空间变换的原理，对倾斜图像进行旋转校正。几何变换由坐标的空间变换和灰度内插这两个基本操作组成，经过对内插效果和运算速度的综合考虑，采用双线性内插方法对旋转后的图像进行反向映射，以保证图像原有数据与形状不会因为旋转和插值而发生畸变。

旋转时以原图像的正中央点为中心，以逆时针方向为旋转正方向，仿射变换采用反向映射的方式，旋转所用的仿射变换公式为：

$\left[\begin{array}{c}v\\ w\\ 1\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\θ}& \sin \mathrm{\θ}& 0\\ -\sin & \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$

其中，(v,w)为原图像中像素的坐标，(x,y)为旋转变换后图像中像素的坐标。

（2）基于投影法和邻域搜索自适应调整的二级网格定位算法

本发明的目的之二是提供一种基于投影法和邻域搜索自适应调整的二级网格定位算法，综合考虑了背景噪声、信号微弱和样点重叠等因素，解决了现有技术中对样点自动定位和定位精度不高的问题。

为了达到上述目的二，经过比较常见的主流样点定位算法，包括手动设置、半自动定位和全自动定位的方法，本发明提出了基于投影法对样点进行初次定位，并生成优化后的网格；然后基于邻域搜索自适应地调整网格，对样点进行二次精确定位，该方法解决问题的技术方案由如下A、B、C、D和E五部分组成。

（A）芯片图像的投影计算

现对经过倾斜校正的准直芯片图像进行投影计算，得到在水平和竖直两个方向上的投影波形。在水平方向的投影表达式为：

$P_H\left(x\right)=\underset{y=0}{\overset{H-1}{\mathrm{\Σ}}}f(x,y)$

在竖直方向的投影表达式为：

$P_V\left(x\right)=\underset{x=0}{\overset{W-1}{\mathrm{\Σ}}}f(x,y)$

其中，芯片图像的宽度为W，高度为H，芯片图像的二维信号为f(x,y)。

考虑到样点周围存在背景噪声，或图像中存在微弱信号的样点和互相重叠的样点，造成投影波形的波峰、波谷位置常有异常的噪声，可能会影响网格定位结果。所以对投影波形P_H(x),P_V(y)进行开闭级联运算操作，可以消除一维信号中异常的凸起峰值和起伏噪声，平滑背景并保留了样点主体区域的信号。

（B）投影波形的自适应二值化

根据投影波形的特征，动态地确定阈值，将得到二值的方波。大于等于该阈值的投影信号认为是波峰，在相应位置上存在样点行（列），将信号值置为1；如果小于该阈值则认为是波谷，在相应位置上是行（列）间空隙，将信号值置为0。自适应二值化算法是根据投影信号的数学统计学信息，动态地选择最佳阈值，以水平投影信号P_H(x)为例，具体描述如下：

a)计算投影信号P_H(x)的平均值，选取该平均值为初始阈值T₀；

b)利用该阈值将投影信号分成两组，信号值小于阈值的划分为group1，大于等于阈值的划分为group2；

c)分别计算group1和group2的均值m₁和m₂，并计算新的阈值T，其中

d)根据第一次迭代得到的阈值T₁，获得判别参数φ=|T₀-T₁|；重复步骤2和3，迭代直到前后连续两次的阈值差|T_n-T_n-1|小于指定的参数φ*λ_T为止，即|T_n-T_n-1|≤φ*λ_T；其中λ_T为阈值参数，一般取值范围为λ_T∈[0.1,0.8]，实验证明λ=0.5较为合适；

e)最终得到在允许范围内的最佳分割阈值T_s，大于等于T_s的信号值置为1，小于T_s的置为0，将原始的投影波形转变为二值的方波。

二值的方波具有原始投影波形不可比拟的优点，使得波峰和波谷宽度的计算变得相当简单、标准化。投影方波波峰的中心就对应着每一行（列）样点的中心，在该位置划分网格线，就能确定样点在水平（竖直）方向的坐标，这坐标将会运用在后续的二次精确定位。波谷中心就对应着行（列）间空隙的中心位置，在该位置划分网格线，就能将每个样点分割到各自的子区域中。

（C）投影方波的分析和处理

步骤C1：精确地读取方波的各个波峰与波谷的宽度，以及波峰中心与波谷中心的位置。

步骤C2：运用数学统计学的原理计算最贴近真实的平均波峰宽度和平均波谷宽度。

将各个波峰宽度存放进一个数组，对该数组进行数学统计，得到平均的波峰宽度。将波峰宽度数组从小到大地排序，然后通过不断地剔除偏离整体较大的数值，得到剩余的符合条件的数值，最后对这些波峰宽度求平均值，即得平均波峰宽度，同样原理求得平均波谷宽度。具体的数值剔除原理如下：

将经过排序的数组第一个值（最小值）剔除，计算数组的方差std₁；将经过排序的数组最后一个值（最大值）剔除，计算数组的方差std₂。比较两种情况下的方差，方差小的情况说明剔除正确。例如，如果剔除最大值后数组的方差更小，即剩余数值的差异性更小，说明被剔除的最大值偏离整体水平比较大。

不断地剔除偏离较大的数值直到满足条件即可终止：

（a）两种剔除方式的方差之差低于一定值，即|std₁-std₂|≤std₀*λ_s；

其中，std₀为原始数组的标准方差，λ_s为方差系数，一般取值范围为λ_s∈[0.1,0.6]，实验证明λ_s=0.5较为合适。

（b）数值中剩余数值的数目到达下限N_min时，必须停止剔除，否则求得的平均波峰宽度不具代表性。（D）网格线优化

大量理论与实验显示投影法在定位时容易造成网格线的冗余和缺失，这一缺点应予以弥补，本发明结合了数学统计学的原理进行改进，对网格线的位置与数量进行调整，删除冗余和填补缺失。

步骤D1:剔除冗余的网格线

冗余网格线的判别需要两组数据：方波波峰宽度和中心距。将上面已经获取的方波波峰宽度数组记为D_p={d_p1,d_p2,...d_pn}，则平均波峰宽度记为

其物理意义是样点直径；分别计算每两个相邻方波波峰的中心距

波峰中心距数组记为D_c={d_c1,d_c2,...d_cn}。采用上述的数学统计学原理，计算出经过剔除的波峰中心距数组的平均值，将平均波峰中心距记为

其物理意义是投影方波的波长。平均波峰宽度

与平均中心距

将应用于冗余剔除，判别原则如下所述：

a)将当前两个相邻波峰的中心距d_ci与平均波峰中心距

进行比较，如果d_ci与平均水平

相差不大，即符合条件

则保留这两个波峰对应的网格线；否则视为冗余网格线并予以剔除。

b)将当前方波波峰宽度d_pi与平均波峰宽度

进行比较，如果d_pi与平均水平

相差不大，即符合条件

则保留该波峰对应的网格线；否则视为样点行（列）形状有误或受到噪声污染，应予以剔除。

其中，λ_p为波峰宽度系数，一般取值范围为λ_p∈[3,10]，实验证明λ_p=5较为合适。

步骤D2：添补缺失的网格线

经过冗余剔除后剩余的网格线存在缺失、且局部变得稀疏，应该进行合理的判别和添补。对剔除后剩余的网格线重新进行计算，相邻网格线距离的数组记为D′_c＝{d′_c1，d′_c2，...d′_cn}，则新的平均网格线距离记为

将当前相邻网格线距离d′_ci与平均网格线距离记为

进行比较，计算当前相邻的两条网格线之间存在的波长个数

（四舍五入地取整数），并添补N_i-1条网格线。

根据上述计算得到的平均波峰中心距，其物理意义为相邻样点平均中心距，可以作为网格初始直径，对每个样点都生成一个大网格，将每个样点分割到各自的子区域中。

（E）自适应调整网格和二次精确定位

基于邻域搜索的原理，自适应地调整网格的尺寸和位置，实现对样点的二次精确定位。水平和竖直两个方向网格线的每个交点对应存在一个样点，交点应该为样点的中心。即使存在或多或少的误差，交点也不会偏离样点太多，一般会落在样点亮斑内部。根据上述计算获得的平均波峰宽度

即得样点平均直径

对每个样点生成一个二次精确定位网格，该小网格的初始直径即为

对于圆形的样点图像，网格设计成圆形；而对于矩形的样点图像，网格设计成矩形。由于样点平均直径

略小于真实的样点直径，故这些网格都能定位样点内部。

只要通过精细地调整每个网格的位置和大小，即可将样点完全套住。网格自适应调整的过程结合了邻域搜索算法，其基本思想是只需要搜索网格中心初始位置的附近区域就能获得最佳匹配区域。自适应调整的最佳结果是将样点刚好套住，信号区域完全落在网格内部，同时套住的非背景区域要尽量少。这需要对每个网格进行精细的调整，包括中心位置的精确定位和尺寸大小的微小调节。精确定位与尺寸调节是交叉进行的，网格中心位置每移动一次，接着就要进行尺寸大小的调节，直到网格区域符合最佳匹配准则。

步骤E1：运用邻域搜索算法，对网格中心位置逐个进行精确调整：

（a）计算网格内部信号强度数据(I_sum,I_mean,I_std)，即圆（矩）形范围内的样点信号总值、均值和标准方差。将网格中心（圆心或者矩形中心）分别移到其八邻域上，在新位置上保持样点直径一致，分别计算样点在新位置上的信号强度数据。从而得到原网格和8个虚拟网格的三组信号强度信息数据，并应用于后续的判别。

（b）找出8个虚拟网格中信号总强度值I_sum的最大值，并跟原位置信号总强度值进行比较。由于原网格与虚拟网格的直径一致，故网格内的像素数目是一致的，只要通过比较信号总强度值I_sum，即可判断出信号局部最强的位置，这也是网格中心位置调整的方向。

（c）由于信号标准方差能反映网格内部的差异信息，故通过比较新位置的信号方差与原位置的方差，可以判断网格中心位置是否调整合理正确。如果符合方差较小的准则，即可用新的中心位置代替原位置。

步骤E2：网格尺寸的调整

在完成一次中心位置的调整后，应该接着调整一次网格的尺寸大小。每次调整网格直径都要确定调整方向，即增大直径还是缩小直径，然后按照直径调整准则进行调节。

直径调整准则是合理直观的，也是相当重要的：

（a）网格的最终直径应当大于初始直径（即样点平均直径

）的80%，并且小于

（b）确保网格不与其他相邻的网格重合；

根据直径调整准则，尺寸大小调整的步骤如下所述：

（a）将样点平均直径设置为网格的初始直径，且现位置网格的内部信号强度数据记为(I⁰ _sum,I⁰ _mean,I⁰ _std)；

（b）将网格直径增大两个像素长度（即半径增大一个像素长度），计算新的虚拟网格的三个信号强度数据(I¹ _sum,I¹ _mean,I¹ _std)；将网格直径缩小两个像素长度（即半径缩小一个像素长度），计算新的虚拟网格的三个信号强度数据(I^-1 _sum,I^-1 _mean,I^-1 _std)；

（c）将I⁰ _mean与(256*256-1)/3进行比较，如果满足条件I⁰ _mean>=(256*256-1)/3，则认为网格内样点信号较强、形状丰满；否则认为网格内样点信号较弱、边缘有缺口。

（d）对于信号较强、形状完好丰满的样点，对应的网格只需要往信号平均值更大、方差更小的方向调整。信号标准方差是一个相当重要的判别依据，网格中含有的背景像素越多，则信号方差越大；当方差最小时，网格内含有的背景像素最少，即调节到最佳尺寸。将增大尺寸的情况与原尺寸的网格进行比较，如果满足条件I¹ _mean>I⁰ _mean,I¹ _std<I⁰ _std，则确定将网格直径增大两个像素长度。将缩小尺寸的情况与原尺寸的网格进行比较，如果满足条件I^-1 _mean=I⁰ _mean,I^-1 _std<I⁰ _std，即尺寸更小的网格内信号平均强度保持不变，而方差小于原网格，则确定将网格直径缩小两个像素长度。其它情况皆保持网格直径不变。

（e）对于信号较弱、边缘有缺口的样点，对应的网格只需要往信号总强度值更大、方差更小的方向调整。将增大尺寸的情况与原尺寸的网格进行比较，如果满足条件I¹ _sum>I⁰ _sum,I¹ _std<I⁰ _std，则确定将网格直径增大两个像素长度。将缩小尺寸的情况与原尺寸的网格进行比较，如果满足条件I^-1 _sum=I⁰ _sum,I^-1 _median>I⁰ _median，即尺寸更小的网格内信号总强度值保持不变，而中值小于原网格，则确定将网格直径缩小两个像素长度。因为弱信号样点对应的网格内含有的背景像素越多，则信号中值越小；缩小尺寸能减少网格中背景像素数目，也使信号中值增大。其它情况皆保持网格直径不变。

步骤E3：snake初始轮廓的生成

本发明采用了snake模型对样点进行边缘提取，以实现样点分割，而snake模型需要将初始轮廓定位到目标样点边缘的附近。因此，本发明巧妙地利用了上述的初始网格和自适应调整后的网格，提出使用双snake分割内部有空洞的样点。初始网格是根据投影网格线的交点产生的，尺寸为样点平均直径，故一般能落在样点亮斑内部，用作内部snake的初始轮廓，非常适合用来追踪样点内部空洞的边缘。自适应调整后的网格将样点刚好套住，信号区域完全落在网格内部，用作外部snake的初始轮廓，适合用来追踪样点外边缘。本发明利用这两种网格成功地将有空洞样点的环形信号区域分割出来，内外snake初始轮廓起着不可替代的作用。

（3）基于Multiple snakes和CNN的样点分割算法

本发明的目的之三是提供一种基于主动轮廓模型（Multiple snakes）和细胞神经网络（CNN）的生物芯片样点分割算法，能够很好地解决形状不规则样点及弱信号样点分割困难的问题。

为了达到上述发明目的之三，本发明提出了利用CNN将样点按信号强弱进行分类，根据贪婪算法优化snake模型，并将Multiple snakes与CNN结合；CNN先学习强信号样点snake的收敛行为，再指导弱信号样点snake的收敛，得到弱信号样点的snake轮廓，最终实现样点的合理分割。

该方法针对复杂背景下生物芯片宏状多目标的特点，基于snake能量模型的思想，构造可变参数的snake神经元能量子模型实现单个样点轮廓跟踪；在此基础上，利用神经网络的智能互联特点、大规模并行判定能力以及分布式信息存储能力，将多个snake神经元网络互联，充分利用芯片样点既具有阵列特性，又是独立的个体的特点，实现样点识别的高速并行处理和分布信息存储的能力。这种新方法采用“先易后难”“边学习边判定”的策略，通过强信号点的形状等先验知识对模型进行训练和学习，获得相同制备条件下形成的生物芯片强信号样点的轮廓特征，然后启发后续弱信号样点轮廓的搜索。这种启发式的引导避免传统主动轮廓跟踪算法在生物芯片样点分割中的不足，极有效的提高分割准确度。

多snake神经元并行网络样点分割算法的主要思想是：

（a）强信号样点指导弱信号样点，利用贪婪算法优化snake模型，结合CNN提高处理速度与准确性。

（b）设定一样点分类模板，利用该模板的CNN对所有样点进行分类，将样点的灰度方差作为初始输入，将强信号样点提取出来，并组成一个新的CNN网络。

（c）根据“先易后难”的思想，CNN的参数模板设定为变量，随着snake的收敛和CNN的学习而更新，强信号样点的snake曲线不断迭代收敛，snake参数也不断改变并达到能量最小化。根据给出的神经元的期望输出，网络跟随着snake神经元进行全局迭代学习。

（d）当边缘最复杂的样点的snake曲线完成收敛时，CNN的学习训练终告结束，得到性能良好的参数模板，并生成供弱信号样点使用的网络。

（e）重新将所有样点当作神经元，利用上面得到的参数模板生成新的CNN，然后将强信号样点对应的神经元看作是记忆神经元，将弱信号样点对应的神经元看作是活动神经元，CNN通过将网络能量最小化指导弱信号样点的snake收敛。最终结合李雅普诺夫方法来判别网络是否稳定，从而结束弱信号样点的snake收敛。

多主动轮廓模型在生物芯片样点分割上的应用：将强信号样点与弱信号样点分开处理，使用不同的snake模型对样点进行分割。强信号样点适合用一般的snake模型；而由于弱信号样点分割困难，需要采用Multiple snake模型来处理，同时借助邻域强信号样点的snake信息。

（A）单独样点的snake模型

由于强信号样点的边缘清晰，图像信息丰富，非常适合使用snake轮廓进行分割，仅需用一般的snake模型即可得到良好的分割效果。故snake曲线的能量E_snake仅考虑单独样点的图像信息，可由个体能量来代表。

（A1）个体能量E_individual由snake曲线的内部能量E_int和外部能量E_ext构成，反映每个单独样点的图像信息，例如梯度信息和光滑度。样点k的个体能量的表达式为

$E_{\mathrm{snake}}^k=E_{\mathrm{in}\mathrm{div}\mathrm{idual}}^k={\mathrm{\&alpha;}}_k{E}_{\mathrm{ext}}^k+{\mathrm{\&beta;}}_k{E}_{\mathrm{int}}^k$

（A2）内部能量E_int是用来保持样点snake曲线的光滑度的，能很好地控制snake曲线上相邻三个连续轮廓点之间的角度。设v_i代表snake曲线上(x_i,y_i)处的轮廓点，该位置的平滑度为S，则S的表达式为

$S\left(v_i\right)=1-\cos {\mathrm{\&theta;}}_i=1-\frac{\stackrel{\&OverBar;}{v_{i-1}{v}_i}\&CenterDot;\stackrel{\&OverBar;}{v_i{v}_{i+1}}}{\left|\stackrel{\&OverBar;}{v_{i-1}{v}_i}\right|\&CenterDot;\left|\stackrel{\&OverBar;}{v_i{v}_{i+1}}\right|}$

其中，v_i-1,v_i和v_i+1为snake曲线上连续的三个轮廓点，

表示v_i-1,v_i之间的向量，

表示v_i,v_i+1之间的向量。样点k的snake曲线平滑度能量E_smooth为

$E_{\mathrm{int}}=E_{\mathrm{smooth}}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}S(x_i,y_i)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}S\left(v_i\right)$

（A3）外部能量E_ext由图像的梯度能量来表示，并使用sobel算子获得合理的图像梯度。图像中位置(x,y)的梯度的灰度强度为I(x,y)，梯度G(x,y)为

$G(x,y)={-|\&dtri;I(x,y)|}^2=-(G_x^2+G_y^2)$

Sobel算子

$G_x=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right],G_y=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& 2& 1\end{array}\right]$

而样点snake的梯度能量需要计算曲线上每一个轮廓点(x_ki,y_ki)的梯度；因此样点的总梯度能量E_gradient或外部能量E_ext为

$E_{\mathrm{ext}}=E_{\mathrm{grad}\mathrm{ient}}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}G(x_i,y_i)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}G\left(v_i\right)$

（A3）权重系数的确定

根据内部能量和外部能量的表达式，可以得到单独的强信号样点的snake能量表达式为：

$E_{\mathrm{snake}}^k={\mathrm{\&alpha;}}_k{E}_{\mathrm{ext}}^k+{\mathrm{\&beta;}}_k{E}_{\mathrm{int}}^k={\mathrm{\&alpha;}}_k\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{G}_k\left(v_i\right)+{\mathrm{\&beta;}}_k\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{S}_k\left(v_i\right)$

权重系数α_k,β_k能够根据样点的snake内部区域的像素灰度方差计算出来。设样点k起始轮廓内部区域D_k的像素灰度方差为σ_k，则有

${\mathrm{\&alpha;}}_k=\frac{{\mathrm{\&sigma;}}_k-{\mathrm{\&sigma;}}_{\min }}{{\mathrm{\&sigma;}}_{\max }-{\mathrm{\&sigma;}}_{\min }},{\mathrm{\&beta;}}_k=1-{\mathrm{\&alpha;}}_k$

其中， ${\mathrm{\&sigma;}}_{\min }=\underset{1<k<M}{\min }{\mathrm{\&sigma;}}_k,{\mathrm{\&sigma;}}_{\max }=\underset{1<k<M}{\max }{\mathrm{\&sigma;}}_k.$

强信号样点拥有清晰的边缘，故权重系数α_k,β_k的计算不需要先验知识的支持；而弱信号样点的边缘不清晰，梯度能量信息不可靠，权重系数α_k,β_k需要根据强信号样点的模板来初始化，否则snake曲线的收敛效果不佳。

（B）多个样点的snake模型

针对分割弱信号样点困难的问题，Multiple snake模型提出了增加集群能量和约束能量的方案，利用强信号样点轮廓的集群能量指导弱信号样点轮廓的更新，并通过周围样点的约束能量来避免收敛错误。

（B1）集群能量E_group通过寻找强信号样点之间的共同特征，来决定弱信号样点轮廓点的更新，从而得到正确的snake形状。由弱信号样点邻域内的强信号样点共同提供一个snake曲线模板，然后使用贝塞尔曲线模型比较模板与目标轮廓的相似性，指导弱信号样点snake轮廓改变形状。

根据弱信号样点邻域内的强信号样点的snake曲线，计算出模板曲线上各点的极坐标，并确定模板曲线形状。设弱信号样点k的3×3邻域内有|S_k|个强信号样点，其中要满足|S_k|≥2，否则扩大邻域尺寸。设

为强信号样点集S_k内snake曲线上轮廓点个数的平均值，并根据S_k的共同特征生成一条有

个的轮廓点的模板曲线。现根据n阶贝塞尔曲线的原理，计算模板曲线上每个轮廓点的极坐标。

对于强信号样点s∈S_k，设s的snake曲线上一轮廓点的极坐标为(r_si,θ_i)，其中

该极坐标的角度为θ_i，则极径为r_si=Q(θ_i)，其中Q(θ)为贝塞尔曲线的伯恩斯坦多项式表达式。

接着可以根据S_k的snake曲线极坐标集

计算出模板曲线t上各个轮廓点的极坐标 $(r_{\mathrm{ti}},{\mathrm{\&theta;}}_i),1\&le;i\&le;\stackrel{\&OverBar;}{N}:$

$r_{\mathrm{ti}}=\underset{s\&Element;S_k}{\mathrm{\&Sigma;}}\frac{r_{\mathrm{si}}}{\left|S_k\right|}$

其中，r_ti表示模板曲线t上轮廓点i的极径，也是角度为θ_i的所有轮廓点的平均极径。

设弱信号样点k的snake曲线上一轮廓点的极坐标为现有由邻域内强信号样点生成的模板曲线t，为强信号样点集S_k内snake曲线上轮廓点个数的平均值，则样点k的集群能量表达式为:

$E_{\mathrm{group}}^k=\underset{i=1}{\overset{\stackrel{\&OverBar;}{N}}{\mathrm{\&Sigma;}}}|r_{\mathrm{ki}}-r_{\mathrm{ti}}|$

（B2）约束能量E_constraint是通过统计相邻样点的空间信息，来提高样点间的约束力，防止snake轮廓互相重叠或者收敛错误。弱信号样点的约束能量也可以作为snake收敛稳定的参考，snake不断收敛直至约束能量的变化量要小于一个阈值。根据北南东西四个方向上与弱信号样点相邻的样点的空间距离C_ki，制定适当的约束能量。

当前弱信号样点k的约束能量为

$E_{\mathrm{constaint}}^k=\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\max {[0,C_{\mathrm{ki}}]}^2$

（B3）权重系数的确定

根据集群能量和约束能量的表达式，可以得到弱信号样点的snake能量表达式为：

$E_{\mathrm{snake}}^k={\mathrm{\&rho;}}_k{E}_{\mathrm{in}\mathrm{div}\mathrm{idual}}^k+{\mathrm{\&gamma;}}_k{E}_{\mathrm{group}}^k+E_{\mathrm{constraint}}^k$

权重系数ρ_k,γ_k可以根据所有弱信号样点的snake曲线的集群能量产生，设当前弱信号样点k的snake曲线的集群能量为则有

${\mathrm{\&rho;}}_k=\frac{E_{\mathrm{group}}^k-E_{\mathrm{group}}^{\min }}{E_{\mathrm{group}}^{\max }-E_{\mathrm{group}}^{\min }},{\mathrm{\&gamma;}}_k=1-{\mathrm{\&rho;}}_k,1\&le;k\&le;W$

其中，芯片图像中样点的总数目为M，强信号样点的数目为S，则弱信号样点的数目是W=M-S。

（C）Snake模型的优化

由于Snake模型对初始轮廓的位置比较敏感，需要依赖其它机制来将曲线约束在样点边缘附近。初始化位置不能离目标样点边缘太远，否则将错误地收敛到附近样点的边缘上，从而不能正确提取目标样点。为了减少对初始轮廓位置的依赖，提高snake的收敛速度和减少迭代的计算量，实现算法的准确性和快速性，现选择结合贪婪算法的原理来优化snake的收敛过程。

贪婪算法的基本思想是不追求最优解，只希望能得到比较满意的解，且一般可以快速地得到较优的解。因为该算法没有为最优解而穷尽所有可能的解，从而大大节省了计算时间。通常贪婪算法都是以当前的参数状况为基础，求出最优解而不考虑可能的整体最优解。

现将贪婪算法应用于优化snake的收敛过程中，在样点k的snake轮廓搜索能量极小的形状时，需要进行一次迭代运算，对snake上的每一个轮廓点的位置进行调整。根据snake能量公式，得到轮廓点i的能量为：

$E_i^k={\mathrm{\&alpha;}}_k{G}_k\left(v_i\right)+{\mathrm{\&beta;}}_k{S}_k\left(v_i\right)$

其中，G_k(v_i)为轮廓点i的梯度，而S_k(v_i)为曲线在轮廓点i处的平滑度。

为了优化snake模型，提出snake的收敛分两个阶段进行，第一阶段使snake快速地收敛到目标样点的附近捕获位置，而第二阶段使snake曲线精细化并尽可能地贴近目标样点边缘。尺度空间由粗到细，降低了计算复杂度与提高了局部位置的精度。

（C1）第一阶段，主要是使snake快速地收敛到目标样点边缘的附近位置，这过程无需考虑snake的光滑度和形状，轮廓点的移动距离可能较大，从而减少对初始轮廓的依赖。此时外部能量占主导地位，快速推动snake贴近样点边缘，其它能量相对较弱，可以不予以考虑。但要保证snake的轮廓点没有因为快速收敛而进入目标样点的内部，因此要对轮廓点的平均距离和像素灰度方差作出限制，这种对轮廓点的限制同样适用于第二阶段，制定合理性判断准则如下：

1）当前轮廓点i与两个相邻轮廓点的距离d_i-1,i,d_i,i+1要与平均距离

相近，即要符合条件

$\left|\frac{d_{i-1,i}+d_{i,i+1}-2\stackrel{\&OverBar;}{d}}{\stackrel{\&OverBar;}{d}}\right|\&le;0.5;$

2）轮廓点i处于当前位置v_n时，snake曲线内部区域的像素灰度方差记为轮廓点i处于下一位置v_n+1时，snake曲线内部区域的像素灰度方差记为

将

与

进行比较，判断轮廓点i的新位置v_n+1是否合理，是否偏离感兴趣区域ROI。如果不符合条件

说明像素灰度方差发生突变，说明该轮廓点移动过快，有可能进入目标样点的内部，应该暂停本次不合理的移动，保持在原位置v_n。其中，λ_σ为方差系数，一般取值范围为λ_σ∈[5,15]，实验证明λ_σ=10较为合适。

（C2）第二阶段，主要是使snake尽可能地贴近目标样点，这过程要细化并平滑snake曲线，轮廓点的移动缓慢。此时，snake已经接近目标样点的边缘，外部能量已经减弱，轮廓曲线的内部能量占主导地位。为了得到能量极小、贴近样点边缘的snake曲线，需要自适应地增删轮廓点，以及精细调整轮廓点的位置。

由于样点的边缘存在局部不规则，snake曲线应该适当细化，在曲率较大的位置自适应地增加轮廓点，在变化较缓、轮廓点距离较均匀的位置自适应地减少轮廓点，从而提高曲线的光滑度和局部位置的合理性。在本阶段中，自适应地增删轮廓点应该与合理性判别准则一起运用，才能保证snake收敛的效果。

自适应的轮廓点增删准则如下：

1）如果局部曲率大于某个阈值，即相邻连续的三个轮廓点组成的角小于某个阈值时，一般取值为75°，则在这三个点之间添加一个轮廓点，使得局部曲率降低。

2）如果相邻连续的三个轮廓点距离较近且近乎在一条直线上时，可将中间的轮廓点删除。

轮廓点的精细调整结合了邻域搜索算法与贪婪算法，其基本思想是通过搜索当前轮廓点的附近区域即可获得能量极小的位置。假设其它各轮廓点都已经处于最佳位置，当前轮廓点i的位置优劣与其它各点都不相关联。计算轮廓点i的八邻域内每个点的能量，将这9个点的能量进行比较，选择具有最小能量的位置作为轮廓点i的下一个位置。轮廓点每次移动一个像素长度，使得snake曲线缓慢收敛到目标样点的边缘。遍历每个轮廓点，完成一次迭代过程。不停地重复迭代，直到该样点的snake曲线趋于稳定为止。

设v_i-1,v_i,v_i+1为snake曲线上三个相邻连续的轮廓点，在计算点v_i的能量并调整位置时，假定相邻的v_i-1和v_i+1都处于最佳位置，这样就将v_i-1,v_i+1对v_i的影响忽略，不必每次迭代都往全局最优解发展，从而简化了计算。

假定v_i点的左右相邻轮廓点v_i-1,v_i+1均已处于各自的最佳位置，可以只考虑v_i八邻域的点，达到局部最优化。同样的原理，遍历所有轮廓点，找到每个点的最佳位置，完成一次迭代。在下一次迭代的时候，以上一次所决定的轮廓点位置作为初始位置，继而进行轮廓点位置的更新和调整，当绝大部分轮廓点的位置基本不再移动时，可以认为迭代运算完成。

因为在运用贪婪算法单独调整每个轮廓点的过程中，假定了相邻的轮廓点已经处于最佳位置，而实际上并不一定最佳，所以只能说相邻轮廓点处于比较好的位置。从而可以看出，尽管贪婪算法使得计算量显著降低，但所求得的解不一定是全局最优解，只能说是比较好的解。

多snake神经元并行网络样点分割算法的实现步骤：

（A）强信号样点的检测

由于要先学习强信号样点的轮廓特征，因此要将强、弱信号样点进行分类，而这种分类问题非常适合运用细胞神经网络来解决，原理如下所述：

（A1）由于强信号样点明显有别于背景区域，轮廓内的背景像素和信号像素之间有很大差别。可以选取像素灰度方差作为判别标准，将方差值较大的样点记为强信号样点，从而区分强和弱信号样点。

微阵列中的样点t起始轮廓的内部区域D_t的像素灰度方差为

${\mathrm{\&sigma;}}_t=\underset{(i,j)\&Element;D_t}{\mathrm{\&Sigma;}}\frac{{(I(i,j)-\stackrel{\&OverBar;}{I})}^2}{N_t}$

其中，内部区域D_t的像素数目为N_t，I(i,j)为D_t内像素点(i,j)的灰度值，

为D_t内所有像素灰度值的平均值。

将有M×N个样点的微阵列的所有灰度方差进行归一化，第i行第j列样点必须满足|σ_ij|≤1，记归一化的方差矩阵为

$\mathrm{\&sigma;}=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_{11}& {\mathrm{\&sigma;}}_{12}& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&sigma;}}_{1N}\\ {\mathrm{\&sigma;}}_{21}& {\mathrm{\&sigma;}}_{22}& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&sigma;}}_{2N}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ {\mathrm{\&sigma;}}_{M1}& {\mathrm{\&sigma;}}_{M1}& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&sigma;}}_{\mathrm{MN}}\end{array}\right]$

（A2）使用代数构造法确定CNN的模板A,B,I的取值，生成专用于强弱样点分类的细胞神经网络，其中A为反馈模板，B为控制模板，I为阈值模板。设用于分类的模板形式为

$A=\left[\begin{array}{ccc}0& 0& 0\\ 0& a& 0\\ 0& 0& 0\end{array}\right],B\left[\begin{array}{ccc}-c& -c& -c\\ -c& b& -c\\ -c& -c& -c\end{array}\right],I=-i$

经过推导与讨论，得到模板参数a,b,i的取值范围

1)a>1/R_x

2)|b-8c|<i

3)i<b-6c

通常取R_x=1，并a要求为大于1的整数，故取值a=2。阈值i取值为2时，b可以取值8，c则取值1。最终得到模板

$A=\left[\begin{array}{ccc}0& 0.& 0\\ 0& 2& 0\\ 0& 0& 0\end{array}\right],B\left[\begin{array}{ccc}-1& -1& -1\\ -1& 8& -1\\ -1& -1& -1\end{array}\right],I=-1$

（A3）将微阵列样点的灰度方差值当作神经元的输入，放入设计好模板的细胞神经网络中，即第i行第j列样点归一化的灰度方差值σ_ij对应着CNN中神经元C(i,j)的输入值u_ij(0)，将神经元C(i,j)初始状态设置为x_ij(0)=0，且满足|u_ij(0)|≤1，|x_ij(n)|≤1。该网络将不断迭代、更新输出并收敛，最终趋于稳定。

参数a=2>1/R_x起着正反馈的作用，一方面使与强信号样点对应的神经元状态值x_ij(n)越来越大，最后x_ij(n)趋于一个大于1的常数值，输出y_ij(n)稳定为1；另外一方面使与弱信号样点对应的神经元状态值x_ij(n)越来越小，最后x_ij(n)趋于小于-1的一常数值，输出y_ij(n)稳定为-1。

（A4）当CNN收敛稳定时，强、弱信号样点的分类完成，可以根据神经元C(i,j)的最终输出y_ij(n)来判断信号的强弱。将强信号样点提取出来，并组成一个新的细胞神经网络。

1）如果y_ij(n)=1，则神经元C(i,j)对应的样点是强信号的。

2）如果y_ij(n)=-1，则神经元C(i,j)对应的样点是弱信号的。

CNN跟一般的人工神经网络一样，在学习分类方面有很高的性能。当需要处理的样点数量越大时，CNN的速度和准确度的优势就更展露无遗。

（B）CNN学习强信号样点snake轮廓的收敛行为

将强、弱信号样点分开的目的是让CNN自组织地学习强信号样点的snake曲线的收敛行为，得到轮廓的模板和CNN的参数A,B,I，然后指导弱信号样点的snake快速收敛。

将强信号样点当作神经元，放入一个新的参数未知的细胞神经网络中，然后给出神经元的期望输出，CNN将跟随着强信号样点snake的收敛进行全局迭代学习。当边缘最复杂的样点的snake完成收敛时，CNN的学习训练终告结束，得到性能良好的模板A,B,I，并生成供弱信号样点使用的网络。

（B1）设强信号样点组成的CNN大小为w×h，设x_ij(n)和y_ij(n)分别为第i行第j列的神经元C(i,j)在第n次迭代的状态和输出，表示CNN的目标输出，本次学习中将

设置为全一矩阵，将A,B,I初始化为零矩阵。通过学习训练决定网络的A,B参数矩阵和I阈值矩阵，以便将CNN网络的输出y_ij(n)转变为期望输出

并在整个过程中学习迭代进行。

（B2）在网络的第n次迭代时，设u_ij(n)表示神经元C(i,j)的当前输入；y_kl(n)表示C(i,j)邻域内的神经元C(k,l)的当前输出，y_kl(n-1)表示上一状态的输出，η为学习率。此时神经元C(i,j)对应强信号样点的snake能量为

将snake能量作为C(i,j)的当前输入，即

根据状态方程和输出方程，可以求得C(i,j)的当前输出。

令 $q=\mathrm{\&eta;}[y_{\mathrm{ij}}\left(n\right)-{\hat{y}}_{\mathrm{ij}}\left(n\right)][1+y_{\mathrm{ij}}\left(n\right)][1-y_{\mathrm{ij}}\left(n\right)],n=\mathrm{1,2,3}\&CenterDot;\&CenterDot;\&CenterDot;,$ 经过推导得到模板的表达式

$A_{\mathrm{kl}}(n+1)=A_{\mathrm{kl}}\left(n\right)-q[y_{\mathrm{kl}}\left(n\right)+A_{\mathrm{kl}}\left(n\right)\frac{y_{\mathrm{kl}}\left(n\right)-y_{\mathrm{kl}}(n-1)}{A_{\mathrm{kl}}\left(n\right)-A_{\mathrm{kl}}\left(n\right)}]$

B_kl(n+1)=B_kl(n)-qu_kl

I_kl(n+1)=I_kl(n)-q

（B3）对于有w×h个神经元的网络，每次迭代都要进行w×h次学习，不断地更新A,B,I，使得输出y_ij(n)不断地逼近期望输出

从而学习如何获得强信号样点轮廓。强信号样点的snake曲线不断地调整收敛，曲线的参数也不断改变并使曲线能量达到最小化。当所有强信号样点的snake轮廓都收敛稳定，那么CNN的学习训练也结束。得到的模板A,B,I能够反映强信号样点的snake收敛特性，非常适合用来指导弱信号样点的snake快速收敛。

（C）经过学习的CNN指导弱信号样点snake轮廓的收敛

利用上述经过学习产生的模板A,B,I，生成用来指导弱信号样点snake收敛的细胞神经网络。将微阵列中样点当作神经元，不分信号强弱，放入模板为A,B,I的细胞神经网络中；然后将强信号样点对应的神经元看作是记忆神经元，它们不参与网络的迭代收敛，故状态、输出保持不变，用来指导弱信号样点的snake收敛；将弱信号样点对应的神经元看作是活动神经元，每次迭代都会改变状态、输入和输出；对应的各个snake神经元不断地调整参数，使得曲线轮廓点向正确的方向移动。同时运用李雅普诺夫方法建立网络能量监控单元，通过监控网络能量判断网络是否稳定；当网络趋于稳定时，得到弱信号样点的snake轮廓以及参数。

（C1）弱信号样点snake参数的初始值由3×3（或5×5）邻域内的强信号样点生成的snake模板α,β,ρ,γ决定。例如，弱信号样点对应的神经元C(i,j)的3×3邻域内有4个强信号样点，则由这4个样点的snake共同生成一个局部适用的模板，即α,β,ρ,γ。利用这组参数进行轮廓收敛迭代，然后不断更新snake参数α,β,ρ,γ，使得snake曲线上的轮廓点向正确的方向移动，snake完美地贴近弱信号样点的边缘。

（C2）李雅普诺夫方法判断网络稳定：

对于大小规模为M×N的CNN，其中神经元C(i,j)(1≤i≤M,1≤j≤N)有来自邻域神经元C(k,l)的外部输入为u_kl，输出变量为y_ij(t)，阈值为I；邻域神经元C(k,l)的输出为y_kl(t)，R_x表示网络常量系数。A(i,j;k,l)表示邻域神经元C(k,l)的输出y_kl(t)与神经元C(i,j)之间的连接权；B(i,j;k,l)表示邻域神经元C(k,l)的输入u_kl(t)与神经元C(i,j)之间的连接权。

根据状态方程（连续时间）定义一个Lyapunov函数，网络能量函数E(t):

$E\left(t\right)=-\frac{1}{2}\underset{(i,j)}{\mathrm{\&Sigma;}}\underset{(k,l)}{\mathrm{\&Sigma;}}\left|A(i,j;k,l)\right|y_{\mathrm{ij}}\left(t\right)y_{\mathrm{kl}}\left(t\right)-\underset{(i,j)}{\mathrm{\&Sigma;}}\underset{(k,l)}{\mathrm{\&Sigma;}}\left|B(i,j;k,l)\right|y_{\mathrm{ij}}\left(t\right)u_{\mathrm{kl}}$

$-\underset{(i,j)}{\mathrm{\&Sigma;}}I{y}_{\mathrm{ij}}\left(t\right)+\frac{1}{2R_x}\underset{(i,j)}{\mathrm{\&Sigma;}}{y}_{\mathrm{ij}}{\left(t\right)}^2$

CNN的能量函数E(t)为单调递减，即

当网络中各神经元的状态基本不变时，网络收敛稳定，能量函数也不再变化。

即当 $\frac{{\mathrm{du}}_{\mathrm{ij}}\left(t\right)}{\mathrm{dt}}=0$ 时，有 $\frac{\mathrm{dE}\left(t\right)}{\mathrm{dt}}=0,\frac{{\mathrm{dy}}_{\mathrm{ij}}\left(t\right)}{\mathrm{dt}}=\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N.$

能量函数的离散形式为：

$E\left(n\right)=-\frac{1}{2}\underset{(i,j)}{\mathrm{\&Sigma;}}\underset{(k,l)}{\mathrm{\&Sigma;}}\left|A(i,j;k,l)\right|y_{\mathrm{ij}}\left(n\right)y_{\mathrm{kl}}\left(n\right)-\underset{(i,j)}{\mathrm{\&Sigma;}}\underset{(k,l)}{\mathrm{\&Sigma;}}\left|B(i,j;k,l)\right|y_{\mathrm{ij}}\left(n\right)u_{\mathrm{kl}}$

$-\underset{(i,j)}{\mathrm{\&Sigma;}}I{y}_{\mathrm{ij}}\left(n\right)+\frac{1}{2R_x}\underset{(i,j)}{\mathrm{\&Sigma;}}{y}_{\mathrm{ij}}{\left(n\right)}^2$

且有E(n)-E(n-1)≤0。根据能量有界定理，能量的最大值为

$E_{\max }=\frac{1}{2}\underset{(i,j)}{\mathrm{\&Sigma;}}\underset{(k,k)}{\mathrm{\&Sigma;}}\left|A(i,j;k,l)\right|+\underset{(i,j)}{\mathrm{\&Sigma;}}\underset{(k,k)}{\mathrm{\&Sigma;}}\left|B(i,j;k,l)\right|+\mathrm{MN}(\frac{1}{2R_x}+|I\left|\right)$

$\underset{t}{\max }\left|E\left(n\right)\right|\&le;E_{\max }$

当E(n)-E(n-1)=0，E(n)稳定为一常数时，说明网络收敛稳定，可以结束迭代。这时弱信号样点的轮廓也收敛得非常完美。由于样点的snake能量在网络中不断传播，经过网络中神经元的前馈控制和反馈修正，使得网络和snake参数不断完善，出色地完成任务。

与现有技术相比，本发明所具有的有益效果为：本发明能自动、快速地识别和分割每个样点，准确地在样点区域内提取灰度强度作为后续处理的数据。基于Hough变换和Radon变换的倾斜快速校正算法创新地使用不同方法分别处理矩形和圆形样点，成功地解决了倾斜检测精度和速度不佳的问题，将倾角检测精度提高到0.01°。基于投影法和邻域搜索自适应调整的二级网格定位算法减少了背景噪声、信号微弱和样点重叠等因素的干扰，解决了对样点自动定位和定位精度不高的问题，为后续的snake分割奠定良好基础。基于Multiple snakes和CNN的样点分割算法很好地结合了snake和CNN的优点，利用优化的snake模型来提取期望的轮廓，联合CNN的自组织学习和并行处理，解决了重叠样点、形状不规则样点以及弱信号样点分割困难的问题。这一套处理流程提高了样点识别的速度、精确度、自动化程度，完全解决了样点识别中人为干预过多的问题，能满足高密度、高通量应用的需要。本发明解决了生物芯片图像倾斜校正、形状不规则样点及弱信号样点分割困难等问题，实现了生物芯片样点的自动识别，适合大规模的生物芯片样点快速分析。

附图说明

图1为本发明一实施例方法流程图；

图2为倾角为2.15°的矩形样点芯片图像的截图；

图3为二值化后矩形样点芯片图像的截图；

图4为经过旋转校正的矩形样点图像的截图；

图5为附图3对应的Hough参数空间沿ρ轴的截图；

图6为倾角为1.55°的圆形样点芯片图像的截图；

图7为对圆形样点芯片图像进行二值化以及提取样点边缘后的截图；

图8为经过旋转校正的圆形样点图像的截图；

图9为对附图7进行第一级Radon变换后对应的参数空间沿ρ轴的截图；

图10为对附图7进行第二级Radon变换后对应的参数空间沿ρ轴的截图；

图11为采用投影法对芯片图像进行初次定位的示意图；

图11（a）为一幅经过倾斜校正的圆形样点芯片图像；

图11（b）为芯片图像在水平方向上的投影波形；

图11（c）为芯片图像在竖直方向上的投影波形；

图11（d）为水平方向上经过自适应二值化的投影方波；

图11（e）为竖直方向上经过自适应二值化的投影方波；

图12为网格线确定行列中心位置的效果图；

图13为网格线将样点分割到子区域的效果图；

图14为对圆形样点进行初次网格定位的效果图；

图15为对圆形样点进行自适应二次定位的效果图；

图16为CNN对样点按信号强弱进行分类的示意图；

图17展示了CNN自组织学习强信号样点轮廓的收敛行为；

图18展示了运用贪婪算法优化的snake模型的收敛过程；

图18（a）为芯片图像中一个边缘有缺陷的样点；

图18（b）到图18（c）展示了snake曲线快速收敛到样点边缘附近；

图18（c）到图18（d）展示了snake曲线被样点边缘上的缺口所吸引；

图18（d）到图18（e）展示了snake曲线精细化并尽可能地贴近目标样点边缘；

图19比较了未经优化的snake模型与经本发明优化的snake模型在处理内部有空洞的样点的效果；

图19（a）为图像中一个内部有空洞的样点；

图19（b）为未经优化的snake模型收敛效果图；

图19（c）为经本发明优化的snake模型收敛效果图；

图20为CNN中强信号样点指导弱信号样点snake收敛的示意图；

图21为本发明对圆形样点的分割效果图；

图22为本发明对矩形样点的分割效果图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚明白，以下将结合具体实施例子，并参照附图，对本发明进一步详细说明。

附图1为本发明的技术解决方案采用C#2.0编程完成整个流程，同时参照附图，对本发明进行更详细说明。

（一）采用改进的Hough变换对矩形样点图像进行倾斜校正

现已知一幅矩形样点芯片图像的倾斜角度为2.15°（以逆时针为旋转正方向），参见附图2。对输入图像进行滤波降噪、灰度增强等图像预处理。采用Otsu算法对预处理后的图像进行二值化，然后进行一次形态学的开闭运算级联操作，得到二值图像如附图3所示。

观察附图3可见：二值化后矩形样点的边缘清晰可见，大部分样点的形状饱满，同一行样点的上（下）边缘呈同一直线分布，非常适合Hough直线检测。对二值图像进行Hough变换，仅计算矩形样点的上下边缘像素的Hough参数和累计值，得到对应的Hough参数空间的正弦图如附图5所示。提取占主导的直线组后，计算平均倾斜角度，得到检测结果为与实际倾角2.15°相差仅为0.04°，算法误差范围为±0.05°。可得最佳的倾斜校正角度为-2.19°，对倾斜图像顺时针旋转校正2.19°，运用双线性内插方法对旋转后的图像进行反向映射，得到校正后的图像如附图4。

由Hough变换的原理可知，倾角为θ_k=2.19°的一族直线占主导地位，对应ρθ参数空间中的累计值也是局部极大值。观察附图5可知，在θ=2.19°附近正弦曲线呈有规律的聚焦状，Hough正弦图中局部极亮点代表矩形样点上下边缘的所在位置。运用改进的Hough变换对矩形样点倾斜图像进行检测，结果证明本发明在运算速度、存储空间和精确度方面有良好的表现。

（二）采用改进的Radon变换对圆形样点图像进行倾斜校正

现已知一幅圆形样点芯片图像的倾斜角度为1.55°（以逆时针为旋转正方向），参见附图6。对输入图像进行滤波降噪、灰度增强等图像预处理。采用Otsu算法对预处理后的图像进行二值化，然后通过形态学腐蚀操作提取圆形样点的边缘，得到二值边缘图像如附图7所示。观察附图7可见：提取到的圆形样点的边缘清晰可见，大部分样点的形状饱满，行列间隙中基本不存在噪声，非常适合Radon投影。

采用二级不同尺度的搜索策略，对边缘图像进行Radon变换。在[-5°,5°]区间内θ以0.1°为步长，进行第一级Radon变换得到初次倾斜角度θ_m=1.6°，得到对应的Radon参数空间的正弦图如附图9所示。观察附图9可知，在θ=1.6°附近正弦曲线有规律地成束分布，Radon正弦图中正弦曲线束代表圆形样点的所在位置。

根据初级倾斜角度θ_m=1.6°，在[-1.50°,1.70°]区间内θ以0.01°为步长，进行第二级Radon变换，得到对应的Radon参数空间的正弦图如附图10所示。搜索得到精确度更高的最佳的倾斜校正角度θ_m=1.58°，与实际倾角1.55°相差仅为0.03°，算法误差范围为±0.04°。可得最佳的倾斜校正角度为-1.58°，对倾斜图像顺时针旋转校正1.58°，运用双线性内插方法对旋转后的图像进行反向映射，得到校正后的图像如附图8。

由Radon变换的原理可知，边缘图像在θ_k=1.58°的方向上的投影为g(ρ,θ_k)，其A_k累加值是局部极大值。运用改进的Radon变换对圆形样点倾斜图像进行检测，结果证明本发明在运算速度、存储空间和精确度方面有良好的表现。

（三）运用二级网格定位算法对生物芯片图像中多样点进行空间定位

现输入一幅经倾斜校正的圆形样点芯片图像，如附图11（a）所示，对生物芯片图像进行投影计算，得到在水平和竖直两个方向上的投影波形P_H(x),P_V(y)，并对波形进行开闭级联运算操作，如附图11（b）和附图11（c）所示。对投影波形进行自适应二值化，运用数学统计学的信息，动态地选择最佳阈值，将投影波形转变为方波，如附图11（d）和附图11（e）所示。

对投影方波进行分析，精确地读取方波的各个波峰与波谷的宽度，以及波峰中心与波谷中心的位置；运用数学统计学的原理计算最贴近真实的平均波峰和波谷宽度。可得水平方向的投影方波P_H(x)的平均波峰宽度为38，平均波谷宽度为37；而垂直方向的投影方波P_V(y)的平均波峰宽度为37，平均波谷宽度为38。

由上述统计得样点平均直径为

行（列）间距为37。对比真实的微阵列数据：样点直径为40，行（列）间距为35；可见，本发明所计算的样点直径略小于真实的样点直径，能够巧妙地将网格初步定位在样点信息区域内部。

根据投影方波P_H(x),P_V(y)波峰中心的位置，划分网格线初次定位到样点行（列）位置；根据投影方波P_H(x),P_V(y)波谷中心的位置，划分网格线将每个样点分割到各自的子区域中。但这两种网格线都存在冗余或缺失，需要进行优化。

结合了数学统计学的原理，对网格线的位置与数量进行调整，通过判断投影方波的波峰宽度和中心距，发现冗余网格线并予以剔除；对剩余的网格线重新进行计算，合理添补缺失的网格线。经过优化和调整的网格线效果图如附图12和附图13所示。附图12展示了根据投影方波的波峰中心划分网格线，确定样点行列中心的位置，水平和竖直方向两个方向网格线的每个交点对应存在着一个样点，且网格线交点应该定位于样点真正中心的邻域；附图13展示了根据方波波谷中心划分网格线，将样点分割到各自的子区域中，为后续的snake样点分割的实现奠定基础。

根据上述计算获得的样点平均直径

以行与列定位网格线的交点为网格中心初始位置，对每个圆形样点生成一个圆形小网格。由于样点平均直径

略小于真实的样点直径，故这些网格都能定位到样点内部。最终得到基于投影法和数学统计学的网格初次定位的效果图，如附图14所示。运用邻域搜索算法，对初次定位的网格中心位置逐个进行精确调整；通过比较网格内部样点信号数据以及样点边缘状况，自适应地确定尺寸调整方向，再按照调整准则进行调整。自适应二次定位的效果图如附图15所示。附图15中位于样点内部的网格是初次定位生成的网格，尺寸略小于样点真实直径，将会用作内部snake的初始轮廓，非常适合用来追踪样点内部空洞的边缘；而位于样点边缘外部的网格是经过精细调整的二次定位网格，能将样点刚好套住，信号区域完全落在网格内部，用作外部snake的初始轮廓，适合用来追踪样点的外轮廓。

（四）多snake神经元并行网络的构造以及在样点分割上的应用

设计一个样点分类专用的模板A,B,I，生成相应的CNN对芯片图像中所有样点进行分类，详细结构见附图16。每个样点对应着一个神经元，将样点归一化的灰度方差作为网络的初始输入，然后网络自动进行迭代收敛。当CNN收敛稳定时，根据网络输出将强信号样点提取出来，并组成一个新的CNN网络。

根据“先易后难”的思想，将分割难度低的强信号样点当作神经元，放入一个新的参数未知的细胞神经网络中。强信号样点组成的CNN大小为w×h，将A,B,I初始化为零矩阵，然后将神经元的期待输出设置为全一矩阵，强信号样点的snake能量作为神经元的输入。CNN将跟随强信号样点snake的收敛进行全局迭代学习，不断自主地调整网络参数A,B,I。当边缘最复杂的样点的snake完成收敛，CNN的学习终告结束，从而得到性能良好的参数模板，并生成供弱信号样点使用的网络。CNN自组织学习强信号样点轮廓的收敛行为，其示意图如附图17所示。

附图18展示了运用贪婪算法优化的snake模型的收敛过程，snake曲线先快速地收敛到样点边缘附近，然后被样点边缘的缺口吸引而凹陷，接着snake曲线精细化并尽可能地贴近目标样点边缘。

观察附图19可见，本发明对snake模型的优化是显而易见的，在运算效率与收敛效果上都有显著的提高。附图19（b）是未经优化的snake模型收敛效果图，未能完全贴近样点边缘，且没有将样点内部的空洞排除；附图19（c）是本发明优化的snake模型收敛效果图，既能很好地贴近样点边缘，又将样点内部的空洞排除在外，读取样点信号时只需计算类环形区域的强度数据。

重新将所有样点当作神经元，利用上面学习得到的参数模板生成新的CNN。将强信号样点对应的神经元看作是记忆神经元，它们不参与网络的迭代收敛，故状态、输出保持不变；将弱信号样点对应的神经元看作是活动神经元，CNN通过将网络能量最小化指导弱信号样点的snake收敛。每次迭代都会改变状态、输入和输出，对应的各个snake神经元不断地调整参数，使得曲线轮廓点向正确的方向移动。最终结合李雅普诺夫方法来判别网络是否稳定，从而结束弱信号样点的snake收敛。

弱信号样点snake参数的初始值由3×3（或5×5）邻域内的强信号样点生成的snake模板α,β,ρ,γ决定，如附图20所示。利用这组参数进行轮廓收敛迭代，然后不断更新snake参数，使得snake完美地贴近弱信号样点的边缘。

当网络能量稳定为一常数时，说明网络收敛稳定，可以结束迭代。这时弱信号样点的轮廓也收敛得非常完美。由于样点的snake能量在网络中不断传播，经过网络中神经元的前馈控制和反馈修正，使得网络和snake参数不断完善，出色地完成任务。

附图21展示了本发明对圆形样点的分割效果图，而附图22展示了本发明对矩形样点的分割效果图。可见，对于微弱信号样点、边缘有缺陷样点和内部有空洞的样点，本发明的分割效果良好。

Claims (13)

1.一种基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，该方法为：

1)若原始生物芯片图像未发生倾斜，即图像准直，则进入4）；若图像发生倾斜，即倾斜角度在[-5°,5°]范围内，则进入2）；若图像倾角超出[-5°,5°]范围，则认为芯片图像的质量不符合要求，不进行后续的分析；

2)将生物芯片图像二值化，计算二值的生物芯片图像的倾斜角度，得到其相反数为最佳的倾斜校正角度；

3)利用所述最佳的倾斜校正角度和图像几何空间变换的原理，对倾斜的原始生物芯片图像进行旋转校正并灰度化，得到准直的芯片灰度图像；

4)对所述准直的芯片灰度图像进行投影计算，得到准直图像在水平和竖直两个方向上的投影波形P_H(x),P_V(y)，将投影波形进行开闭级联运算操作和自适应二值化处理，得到二值投影方波；

5)二值投影方波P_H(x)的波峰中心对应着每一行样点的中心，在所述中心划分网格线，确定样点行在水平方向的坐标；二值投影方波P_H(x)的波谷中心对应着样点行间空隙的中心位置，在所述中心位置划分网格线，将每个样点分割到各自的子区域中；运用上述同样的方法处理二值投影方波P_V(y)；

6)读取上述二值投影方波的各个波峰和波谷的宽度，以及波峰中心与波谷中心的位置；运用数学统计学原理，计算出贴近真实的平均波峰宽度和平均波谷宽度，以及平均波峰中心距；

7)结合数学统计学原理，通过判断每条网格线对应的方波的波峰宽度和中心距，发现冗余网格线并予以剔除；对剩余的网格线重新进行计算，添补缺失的网格线；对网格线的位置与数量进行调整，生成优化后的网格线，完成对样点的初次定位；

8)将步骤6）获得的平均波峰宽度用作二次精确定位网格的初始直径，对每个样点生成一个二次精确定位初始网格；对于圆形样点的芯片图像，网格设计成圆形；对于矩形样点的芯片图像，网格设计成矩形；

9)基于邻域搜索和自适应调整的原理，对样点进行二次精确定位：运用邻域搜索算法，对上述产生的二次精确定位网格的中心位置逐个进行调整；通过比较网格内部样点信号数据，确定尺寸调整方向，再按照直径调整准则进行调整；使得网格最终能将样点套住，信号区域完全落在网格内部，同时套住的背景区域像素数目达到最小值；

10)根据贪婪算法优化主动轮廓模型，利用细胞神经网络CNN将芯片样点按信号强弱进行分类；将多主动轮廓模型Multiple snakes与CNN结合，CNN先学习强信号芯片样点snake的收敛行为，并得到强信号样点的snake轮廓；再指导弱信号芯片样点snake的收敛，得到弱信号芯片样点的snake轮廓，最终实现芯片样点信号区域的分割，提取芯片样点信号区域内灰度强度作为后续处理的数据；

11)计算每个样点信号区域内部的灰度总值、平均值和方差，统计得到芯片样点的信号值；计算样点snake轮廓外部环形区域的灰度强度信息，统计出每个样点的背景值，最后输出生物芯片样点分析数据。

2.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤2）中，计算生物芯片二值图像的倾斜角度时，对矩形样点图像采用改进的Hough变换检测倾斜角度，对于圆形样点图像则采用改进的Radon变换：

a)采用改进的Hough变换检测矩形样点的芯片图像的倾斜角度过程如下：

1)对生物芯片图像进行灰度化，然后进行一次形态学的开闭运算级联操作；接着运用Otsu算法将芯片图像二值化；

2)计算矩形样点的上下边缘像素的Hough参数和累计值；以二值芯片图像的中心为计算原点，以对角线长度为ρ轴最大值，步长为1；以5°为θ轴的最大值，步长为0.01°，计算Hough参数空间；统计ρθ参数空间中各曲线交点的信号累加值，对应着各个角度上特征点的累计数量；

3)选择图像宽度的1/10为最小直线长度A_min；遍历Hough参数空间，如果参数空间坐标(ρ_i,θ_j)的累加值为A(i,j)，代表直线xcosθ_j+ysinθ_j=ρ_i的累计长度为A(i,j)；如果线长度A(i,j)小于A_min，则表明该直线不占主导地位，舍去不记录；反之，如果线长度A(i,j)不小于A_min，则进一步判别它是否为局部极大值；根据局部峰值半径r，将Hough参数空间中A(i,j)的r×r邻域值与当前累加值A(i,j)比较，如果A(i,j)在局部区域内是极大值，则表明该直线起主导作用，保存到直线组并用于下一步计算；

4)将上述步骤3）中提取到的直线组按照累计长度从大到小排序，从而得到最大累计长度A_max，分别计算累计长度较长的前10条直线的相对长度

提取相对长度不小于0.5的直线，计算它们的平均倾斜角度

则最佳的倾斜校正角度为平均倾斜角度的相反数

$\stackrel{\&OverBar;}{\mathrm{\&theta;}}=\frac{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_m{R}_m}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{R}_m}$

其中，M为相对长度不小于0.5的直线的条数，R_m为相对长度不小于0.5的直线的相对长度，θ_m为相对长度不小于0.5的直线的倾斜角度；

b)采用改进的Radon变换检测圆形样点的芯片图像的倾斜角度过程如下：

1)选择Otsu算法将芯片图像二值化，得到二值的芯片图像A；通过形态学腐蚀操作提取圆形样点的边缘，使用最小的结构元B对二值图像A进行腐蚀，再用A减去腐蚀图像ε_B(A)，得到二值边缘图像β(A)=A-ε_B(A)，且二值边缘图像边缘宽度为1；

2)在倾斜角度[-5°,5°]区间内θ以0.1°为步长，对二值边缘图像进行第一级Radon变换，得到宽度为101的Radon参数空间g(ρ,θ)；固定θ值为θ_k，计算投影g(ρ,θ_k)的累计值

遍历所有累计值，最大累计值A_m所对应的θ_m为初次倾斜角度；记初次倾斜角度为

在

区间θ以0.01°为步长，对二值边缘图像进行第二级Radon变换，得到宽度为301的Radon参数空间g(ρ,θ)；通过计算投影的累计值并比较出最大值，从而得到对应的二次倾斜角度其相反数即为最佳的倾斜校正角度。

3.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤4）中，对投影信号P_H(x),P_V(y)分别进行自适应二值化处理的过程如下：

1)计算投影信号P_H(x)或P_V(y)的平均值，选取该平均值为初始阈值T₀；

2)利用所述初始阈值将投影信号分成两组，信号值小于阈值的划分为group1，不小于阈值的划分为group2；

3)分别计算group1和group2的均值m₁和m₂，并计算新的阈值T，其中

4)根据第一次迭代得到的阈值T₁，获得判别参数φ=|T₀-T₁|；重复步骤2和3，迭代直到前后连续两次的阈值差|T_n-T_n-1|小于指定的参数φ*λ_T为止，即|T_n-T_n-1|≤φ*λ_T；其中λ_T为阈值系数，取值范围为λ_T∈[0.1,0.8]；

5)得到最佳分割阈值T_s，大于等于T_s的信号值置为1，小于T_s的置为0，将原始的投影波形转变为二值的方波。

4.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤6）中，运用数学统计学方法分别计算平均波峰宽度和平均波谷宽度，二值投影方波的平均波峰宽度或平均波谷宽度的计算过程如下：

将所述步骤6）所获取的各个波峰宽度或波谷宽度存放在一个数组中，将波峰宽度数组或波谷宽度数组中的数值从小到大地排序，然后通过不断地剔除偏离整体较大的数值，得到剩余的符合条件的数值，最后对这些波峰宽度数值求平均值，得到平均波峰宽度或平均波谷宽度；其中剔除偏离整体较大的数值的过程为：将经过排序的波峰宽度数组或波谷宽度数组的第一个数值剔除，即剔除最小值，计算剔除最小值后数组的方差；将经过排序的波峰宽度数组或波谷宽度数组的最后一个数值剔除，即剔除最大值，计算剔除最大值后数组的方差；比较两种情况下的方差，方差较小情况下的剔除是正确的；对于经过排序的波峰宽度数值或波谷宽度数值，以此类推。

5.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤7）中，剔除冗余的网格线和添补缺失的网格线的过程为：

1)将当前两个相邻波峰的中心距d_ci与平均波峰中心距

进行比较，如果d_ci与平均波峰中心距相差不大，即符合条件

则保留这两个波峰对应的网格线；否则视为冗余网格线并予以剔除；其中，λ_c为波峰中心距系数，取值范围为λ_c∈[5,15]；

2)将当前方波波峰宽度d_pi与平均波峰宽度

进行比较，如果d_pi与平均水平

符合条件

则保留该波峰对应的网格线；否则视为样点行或列形状有误或受到噪声污染，应予以剔除；其中，λ_p为波峰宽度系数，取值范围为λ_p∈[3,10]；

3)对剔除后剩余的网格线重新进行计算，获取相邻网格线距离的数组D′_c＝{d′_c1，d′_c2，...d′_cn}，并计算新的平均网格线距离

将当前相邻网格线距离d′_ci与平均网格线距离

进行比较，计算当前相邻的两条网格线之间存在的波长整数个数

并添补N_i-1条网格线。

6.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤9）中，运用邻域搜索算法，对网格中心位置逐个进行调整的过程为：

1)计算网格内部信号强度数据(I_sum,I_mean,I_std)，I_sum,I_mean,I_std分别表示圆形或矩形网格范围内的样点信号总值、均值和标准方差；将网格中心分别移到网格中心的八邻域上，在新位置上保持样点直径一致，分别计算样点在新位置上的信号强度数据，得到原网格和8个虚拟网格的三组信号强度信息数据；

2)找出8个虚拟网格中样点信号总值I_sum的最大值，并跟原位置样点信号总值进行比较；判断出信号局部最强的位置，即确定网格中心位置调整的方向；

3)比较新位置的信号方差与原位置的信号方差，判断网格中心位置是否调整合理正确，如果新位置的信号方差较小，则用新的网格中心位置代替原位置。

7.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤9）中，自适应调整网格尺寸的步骤为：

1)根据步骤8）所生成的二次精确定位网格，计算所述初始网格的内部信号强度数据(I⁰ _sum,I⁰ _mean,I⁰ _std)；其中，I⁰ _sum,I⁰ _mean,I⁰ _std分别表示初始网格内的样点信号总值、均值和标准方差；

2)将二次精确定位网格的直径增大两个像素长度，计算增大直径后的虚拟网格的样点信号总值、均值和标准方差(I¹ _sum,I¹ _mean,I¹ _std)；将网格直径缩小两个像素长度，计算缩小直径后的虚拟网格的样点信号总值、均值和标准方差(I^-1 _sum,I^-1 _mean,I^-1 _std)；

3)将I⁰ _mean与(256*256-1)/3进行比较，如果满足条件I⁰ _mean≥(256*256-1)/3，则认为网格内的样点信号较强、形状丰满；否则认为网格内的样点信号较弱、边缘有缺口；

4)对于信号较强、形状丰满的样点，如果满足条件

则将网格的直径增大两个像素长度；如果满足条件

则将网格的直径缩小两个像素长度；其它情况皆保持网格直径不变；

5)对于信号较弱、边缘有缺口的样点，如果满足条件

则将网格的直径增大两个像素长度；如果满足条件

则将网格的直径缩小两个像素长度；其它情况皆保持网格直径不变；

6)每完成一次网格中心位置的调整，接着调整一次网格的尺寸大小。

8.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤10）中，运用贪婪算法优化主动轮廓模型的过程为：

1)通过搜索snake曲线上当前轮廓点的附近区域获得能量极小的位置；假定snake曲线上相邻的轮廓点都已经处于最佳位置，当前轮廓点位置的优劣与其它各轮廓点都不相关联；在当前轮廓点的八邻域内，通过搜索局部能量极小的位置来调整当前轮廓点，从而达到局部最优化；

2)将snake的收敛分成两个阶段进行：第一阶段贪婪地仅以外部能量为主导，不考虑snake曲线的光滑度和形状，通过将snake曲线能量最小化计算出轮廓点的移动距离，使snake收敛到目标样点的附近捕获位置；第二阶段仅以内部能量为主导，不考虑外部能量，自适应地增删轮廓点，调整各轮廓点的位置，使snake曲线贴近目标样点的边缘。

9.根据权利要求1所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，所述步骤10）中，CNN将芯片样点按信号强弱进行分类的过程为：

1)计算芯片灰度图像中样点t的snake起始轮廓的内部区域D_t的像素灰度方差：

其中，内部区域D_t的像素数目为N_t，I(i,j)为D_t内像素点(i,j)的灰度值，

为D_t内所有像素灰度值的平均值；将微阵列中M×N个样点的所有灰度方差进行归一化，第i行第j列样点必须满足|σ_ij|≤1，得到归一化的方差矩阵

$\mathrm{\&sigma;}=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_{11}& {\mathrm{\&sigma;}}_{12}& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&sigma;}}_{1N}\\ {\mathrm{\&sigma;}}_{21}& {\mathrm{\&sigma;}}_{22}& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&sigma;}}_{2N}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ {\mathrm{\&sigma;}}_{M1}& {\mathrm{\&sigma;}}_{M1}& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&sigma;}}_{\mathrm{MN}}\end{array}\right];$

2)使用代数构造法确定CNN的模板A,B,I的取值，生成用于强弱样点分类的细胞神经网络；

3)将微阵列样点的灰度方差值作为神经元的输入，即第i行第j列样点归一化的灰度方差值σ_ij对应着CNN中神经元C(i,j)的输入值u_ij(0)，将神经元C(i,j)初始状态设置为x_ij(0)=0，且满足|u_ij(0)|≤1，|x_ij(n)|≤1；所述网络将随着不断迭代而收敛，最终趋于稳定；

4)当CNN收敛稳定时，强、弱信号样点的分类完成，根据神经元C(i,j)的最终输出y_ij(n)来判断信号的强弱，将强信号样点提取出来，并组成一个新的细胞神经网络。

10.根据权利要求9所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，CNN学习强信号芯片样点snake收敛行为的过程为：

1)将所述的强信号样点当作神经元，放入一个新的参数未知的细胞神经网络中；设神经元C(i,j)对应强信号样点的snake曲线的能量为

将snake能量作为C(i,j)的当前输入，即

然后给出神经元C(i,j)的期望输出

通过学习训练决定网络的参数矩阵A,B和阈值矩阵I，将CNN网络的输出y_ij(n)转变为期望输出

并在整个过程中学习迭代进行；

2)对于有w×h个神经元的CNN，每次迭代都要进行w×h次学习，不断地更新A,B,I，使得输出y_ij(n)不断地逼近期望输出

从而学习获得强信号样点轮廓；当所有强信号样点的轮廓都收敛稳定，那么CNN的学习训练结束；得到的模板A,B,I用于指导弱信号样点的snake快速收敛。

11.根据权利要求10所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，经过学习的CNN指导弱信号芯片样点snake收敛的过程为：

1)根据所述经过学习产生的模板A,B,I，生成用来指导弱信号样点的snake收敛的细胞神经网络；将微阵列中样点当作神经元，不分信号强弱，放入经过学习的网络中；将强信号样点对应的神经元看作是记忆神经元，记忆神经元不参与网络的迭代收敛，故状态、输出保持不变，用来指导弱信号样点的snake收敛；将弱信号样点对应的神经元看作是活动神经元，且网络的每次迭代都会改变状态、输入和输出；

2)弱信号样点snake参数的初始值由3×3或5×5邻域内的强信号样点生成的snake模板α,β,ρ,γ决定；利用这组参数进行轮廓收敛迭代，然后不断更新snake参数α,β,ρ,γ，使得snake曲线上的轮廓点向向样点边缘的方向移动，贴近弱信号样点的边缘；同时运用李雅普诺夫方法建立网络能量监控单元，通过监控网络能量判断网络是否稳定；当网络趋于稳定时，得到弱信号样点的snake轮廓以及参数。

12.根据权利要求3所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，λ_T=0.5。

13.根据权利要求5所述的基于主动轮廓模型和细胞神经网络的生物芯片分析方法，其特征在于，λ_c=8，λ_p=5。

Images