CN103473776B - 一种比较图像清晰度的方法、系统及自动聚焦控制方法 - Google Patents

Abstract

本发明涉及高通量基因测序领域，提供了一种比较多个测序图像清晰度的方法和系统及自动聚焦控制方法。所述比较多个测序图像清晰度的方法通过对多个测序图像分别进行谐波分析处理，然后根据谐波分析处理图的纹理特征分别拟合出闭合曲线，并根据闭合曲线得到各闭合曲线对应的图形与其对应的测序图像之间的面积比值，进而通过比较面积比值的大小确定清晰度最高的测序图像。本发明的方法和系统相对于现有技术中具有量化精度更高的优点；此外，本发明还基于上述清晰度比较方法提供了一种自动聚焦控制方法，该方法能够精确定位自动聚焦位置。

Description

一种比较图像清晰度的方法、系统及自动聚焦控制方法

技术领域

本发明涉及高通量基因测序领域，更具体地说，涉及一种比较多个测序图像清晰度的方法和系统，以及一种自动聚焦控制方法。

背景技术

高通量基因测序图像承载着丰富的人类基因信息，极高的清晰度要求已经成为衡量实验成功的重要环节。高通量基因测序图像极具纹理密度多样性的特性，在检测图像质量过程中，属典型的无参考图像质量检测。无参考图像质量评价方法就是在没有一个可以进行参考和对比的原始图像的情况下，得出一个与人类视觉系统的视觉感知相一致的质量分数值的方法。目前，对数字图像清晰度的无参考图像质量评价方法主要有灰度平均梯度法、拉普拉斯算子和这两种评价方法。

灰度平均梯度法(GrayMeanGrads，GMG)是分别将图像长度和宽度方向上的相邻像素灰度值的差求平方和再求均方根，它能较好的反映图像的对比度变化特征，其值越大表示图像越清晰，高通量基因测序图像中的信号强度越突出，易被识别，说明图像质量越好。

令g为M×N的图像矩阵，则灰度平均梯度值表达式如下：

$GMG=\frac{1}{(M-1)(N-1)}\underset{i=1}{\overset{M-1}{\Σ}}\underset{j=1}{\overset{N-1}{\Σ}}\sqrt{\frac{{\[g(i,j+1)-g(i,j)\]}^2+{\[g(i+1,j)-g(i,j)\]}^2}{2}}.$

拉普拉斯算子和(LaplaceoperatorSum，LS)是将图像范围内每一个像素点的8邻域微分值求和，其值同样能较好的反映图像的对比度的变化特征，其值越大表示图像越清晰，高通量基因测序图像中的信号强度越突出，易被识别，说明图像质量越好。

令g为M×N的图像矩阵，则拉普拉斯算子和表达式如下：

$LS=\frac{1}{(M-2)(N-2)}\underset{i=2}{\overset{M-1}{\Σ}}\underset{j=2}{\overset{N-1}{\Σ}}|9g(i,j)-\underset{p=-1}{\overset{1}{\Σ}}\underset{q=-1}{\overset{1}{\Σ}}g(i+p,j+q)|.$

由以上模型可知，当图像越清晰，轮廓越鲜明，则每一个像素点附近的灰度值变化越大，即灰度平均梯度值与拉普拉斯算子和越大，可以说明高通量基因测序图像中的信号信息更易被识别。

但是，因为高通量基因测序图像的灰度级非常大，利用灰度评价梯度法与拉普拉斯算子和得到的结果的数量级数都很大，难以制定出合理的评价标量。此外，灰度评价梯度法与拉普拉斯算子和受图像中发光点密度的影响很大，对于清晰度不同的高通量基因测序图像，利用灰度评价梯度法或拉普拉斯算子和法得到的值却相同或相差不大。即，利用灰度评价梯度法或拉普拉斯算子和法，在对模糊程度不同的测序图像进行质量评价时存在着量化精度不高的问题。在现有技术中，基于灰度评价梯度法或拉普拉斯算子和的自动聚焦控制方法，所得的自动聚焦位置不精确。

因此，需要一种新的比较多个测序图像清晰度的方法和系统，上述方法和系统能够提高对测序图像进行质量评价时的量化精度，还需要一种新的能够精确定位自动聚焦位置的自动聚焦控制方法。

发明内容

本发明的目的在于提供一种新的比较多个测序图像清晰度的方法和系统，旨在解决现有技术对模糊程度不同的测序图像进行质量评价时存在着量化精度不高的问题。

为了实现发明目的，本发明提供了一种比较多个测序图像清晰度的方法，包括以下步骤：

A.对多个测序图像分别进行谐波分析处理，得多个与测序图像相对应的谐波分析处理图；

B.根据所述与测序图像相对应的谐波分析处理图的纹理特征分别拟合出多个与测序图像相对应的闭合曲线；

C.确定所述闭合曲线所对应的图形和相应的测序图像之间的面积比值；

D.根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像。

其中，所述谐波分析处理可为傅立叶变换、小波变换或沃尔什变换。

优选的，所述谐波分析处理为傅立叶变换；所述闭合曲线为椭圆曲线或圆形曲线；所述闭合曲线对应的图形为椭圆或圆。

其中，所述步骤A可包括以下步骤：

A1.根据式Ⅰ对多个测序图像分别进行预处理，得多个预处理后的图像；

A2.根据式Ⅱ对多个预处理后的图像分别进行傅立叶变换，得多个与测序图像相对应的傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}g(x,y)e^{-j2\mathrm{\π}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

其中，所述步骤B可包括以下步骤：

B1.分别对多个与测序图像相对应的傅立叶变换图进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

B2.将多个与测序图像相对应的傅立叶同态变换图中灰度级小于等于a的值重置为0；

B3.利用半径为r的圆形结构元对经步骤B2所得的结果图像进行闭运算；

B4.利用半径为R的圆形结构元对经步骤B3所得的结果图像进行开运算；

B5.利用边缘检测方法对经步骤B4所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线；

所述10＜a＜15，R＞r。

其中，步骤B所得多个与测序图像相对应的椭圆曲线的一般方程为：

Ax²+Bxy+Cy²+Dx+Ey+F＝0Ⅲ；

步骤B所得多个与测序图像相对应的圆形曲线的一般方程为：

x²+y²+Dx+Ey+F＝0(D²+E²-4F＞0)。

其中，所述步骤C包括以下步骤：

C1.分别根据步骤B所得椭圆曲线或圆形曲线的一般方程，计算出各椭圆曲线或圆形曲线所对应的椭圆或圆的面积；

C2.确定各椭圆曲线或圆形曲线所对应的椭圆或圆和相应的测序图像之间的面积比值。

其中，所述步骤B1为：根据式Ⅳ对多个与测序图像相对应的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。

其中，步骤B5所述边缘检测方法为Sobel边缘检测算法、canny边缘检测算法、robert边缘检测算法、prewite边缘检测算法或kirsch边缘检测算法。

其中，所述R等于1.5r、2r、2.5r或3r。

为了更好的实现本发明的目的，本发明还提供了一种比较多个测序图像清晰度的系统，包括：谐波分析处理单元、闭合曲线拟合单元、面积比值计算单元和清晰度判定单元；

所述谐波分析处理单元用于对测序图像进行谐波分析处理得谐波分析处理图；

所述闭合曲线拟合单元用于根据所述谐波分析处理图的纹理特征拟合出闭合曲线；

所述面积比值计算单元用于确定所述闭合曲线对应的图形和相应的测序图像之间的面积比值；

所述清晰度判定单元用于根据面积比值的大小确定所述测序图像的清晰度。

其中，所述谐波分析处理单元用于对测序图形进行谐波分析处理得谐波分析处理图可具体为：所述谐波分析处理单元用于对测序图像进行傅立叶变换得傅立叶变换图；或所述谐波分析处理单元用于对测序图像进行小波变换得小波变换图；或所述谐波分析处理单元用于对测序图像进行沃尔什变换得沃尔什变换图。

优选的，所述谐波分析处理单元用于对测序图像进行傅立叶变换得傅立叶变换图；所述闭合曲线拟合单元用于根据所述傅立叶变换图的纹理特征拟合出椭圆曲线或圆形曲线；所述面积比值计算单元用于确定所述椭圆曲线或圆形曲线所对应的椭圆或圆和相应的测序图像之间的面积比值。

其中，比较多个测序图像清晰度的系统还可包括预处理单元，所述预处理单元用于根据式Ⅰ对测序图像进行预处理，得预处理后的图像；所述谐波分析处理单元用于根据式Ⅱ对预处理后的图像进行傅立叶变换，得傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}g(x,y)e^{-j2\mathrm{\π}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

其中，所述闭合曲线拟合单元包括灰度级压缩模块、重置模块、闭运算模块、开运算模块和边缘提取曲线拟合模块；

所述灰度级压缩模块用于对傅立叶变换图进行灰度级压缩，得傅立叶同态变换图；

所述重置模块用于将傅立叶同态变换图中灰度级小于等于a的值重置为0；

所述闭运算模块用于利用半径为r的圆形结构元对经重置模块处理后所得的结果图像进行闭运算；

所述开运算模块用于利用半径为R的圆形结构元对经闭运算模块处理后所得的结果图像进行开运算；

所述边缘提取曲线拟合模块用于利用边缘检测方法对经开运算模块处理后所得的结果图像进行边缘提取，进而进行曲线拟合，得与测序图像相对应的椭圆曲线或圆形曲线；

所述10＜a＜15，R＞r。

其中，所述闭合曲线拟合单元拟合出的椭圆曲线的一般方程为：

Ax²+Bxy+Cy²+Dx+Ey+F＝0Ⅲ；

所述闭合曲线拟合单元拟合出的圆形曲线的一般方程为：

x²+y²+Dx+Ey+F＝0(D²+E²-4F＞0)。

其中，所述面积比值计算单元包括面积计算模块、比值计算模块；

所述面积计算模块用于计算闭合曲线所对应的图形的面积；

所述比值计算模块用于确定所述闭合曲线所对应的图形的面积和相应的测序图像的面积之间的比值。

其中，所述灰度级压缩模块根据式Ⅳ对傅立叶变换图进行灰度级压缩，得傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。

其中，所述边缘检测方法为Sobel边缘检测算法、canny边缘检测算法、robert边缘检测算法、prewite边缘检测算法或kirsch边缘检测算法。

其中，所述R等于1.5r、2r、2.5r或3r。

本发明的另一目的在于提供一种新的自动聚焦控制方法，旨在解决现有技术中的自动聚焦控制方法得到的自动聚焦位置不精确的问题。

一种自动聚焦控制方法，包括以下步骤：

A.调焦模块多次调节成像装置与采图位置之间的距离，并在每个距离均发出采集图像的指令给采图模块；

B.采图模块在接收到采集图像的指令后，利用成像装置采集测序图像，进而得到针对同一采图位置的多个测序图像；

C.图像处理模块在接收到经步骤B得到的多个测序图像后，根据上述的比较多个测序图像清晰度的方法，确定清晰度最高的测序图像，进而获得最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块；

D.调焦模块根据接收到的指令，控制成像装置与采图位置之间的距离为最佳聚焦距离。

其中，所述图像处理模块包括谐波分析处理单元、闭合曲线拟合单元、面积比值计算单元和最佳聚焦距离判定单元；

所述步骤C包括以下步骤：

C1.谐波分析处理单元对经步骤B所得的多个测序图像分别进行谐波分析处理，得多个与测序图像相对应的谐波分析处理图；

C2.闭合曲线拟合单元分别根据所述谐波分析处理图的纹理特征拟合出多个与测序图像相对应的闭合曲线；

C3.面积比值计算单元确定所述闭合曲线所对应的图形和其相应的测序图像之间的面积比值；

C4.所述最佳聚焦距离判定单元根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像，进而确定最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块。

其中，所述步骤C1中的谐波分析处理可为傅立叶变换、小波变换或沃尔什变换。

优选的，所述步骤C包括以下步骤：

C1.谐波分析处理单元对经步骤B所得的多个测序图像分别进行傅立叶变换，得多个与测序图像相对应的傅立叶变换图；

C2.闭合曲线拟合单元分别根据所述傅立叶变换图的纹理特征拟合出多个与测序图像相对应的椭圆曲线或圆形曲线；

C3.面积比值确定单元确定所述椭圆曲线或圆形曲线所对应的椭圆或圆和其相应的测序图像之间的面积比值；

C4.所述最佳聚焦距离判定单元根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像，进而确定最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块。

其中，所述图像处理模块还包括预处理单元；所述步骤C1包括以下步骤：

C11.所述预处理单元根据式Ⅰ对经步骤B所得的多个测序图像分别进行预处理，得多个预处理后的图像；

C12.所述谐波分析处理单元根据式Ⅱ对多个预处理后的图像分别进行傅立叶变换，得多个与测序图像相对应的傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}g(x,y)e^{-j2\mathrm{\π}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

其中，所述步骤C2包括以下步骤：

C21.对步骤C1所得的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

C22.将经步骤C21所得的傅立叶同态变换图中灰度级小于等于a的值重置为0；

C23.利用半径为r的圆形结构元对经步骤C22所得的结果图像进行闭运算；

C24.利用半径为R的圆形结构元对经步骤C23所得的结果图像进行开运算；

C25.利用边缘检测方法对经步骤C24所得的结果图像进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线；

所述10＜a＜15，R＞r。

其中，所述步骤C2所得多个与测序图像相对应的椭圆曲线的一般方程为：

Ax²+Bxy+Cy²+Dx+Ey+F＝0Ⅲ；

所述步骤C2所得多个与测序图像相对应的圆形曲线的一般方程为：

x²+y²+Dx+Ey+F＝0(D²+E²-4F＞0)。

其中，所述步骤C3包括以下步骤：

C31.根据步骤C2所得的闭合曲线的一般方程，计算出各闭合曲线所对应的椭圆的面积；

C32.确定各闭合曲线所对应的图形和相对应的测序图像之间的面积比值。

其中，所述步骤C21为：根据式Ⅳ对步骤C1所得的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。

其中，步骤C25所述边缘检测方法为Sobel边缘检测算法、canny边缘检测算法、robert边缘检测算法、prewite边缘检测算法或kirsch边缘检测算法。

其中，所述R等于1.5r、2r、2.5r或3r。

与现有技术相比，本发明的方法和系统创造性的在高通量测序领域通过对多个测序图像分别进行谐波分析处理，然后谐波分析处理图的纹理特征分别拟合出闭合曲线，并根据闭合曲线得到其所对应的图形与其对应的测序图像之间的面积比值，进而通过比较面积比值的大小确定清晰度最高的测序图像，本发明的方法和系统相对于现有技术中具有量化精度更高的优点；此外，本发明还基于上述清晰度比较方法提供了一种自动聚焦控制方法，该方法能够精确定位自动聚焦位置。

附图说明

图1a是本发明第二实施例中的测序图像。

图1b是本发明第二实施例中获得的前期傅立叶变换图。

图1c是本发明第二实施例中经StepA2’后获得的傅立叶变换图。

图1d是本发明第三实施例中经闭运算后得到的结果图像。

图1e是本发明第三实施例中经开运算后得到的结果图像。

图1f是以图1a作为测序图像，按照第一实施例和第三实施例的方法处理后，得到的结果图。

图2至图4是本发明第一具体实施例中的高通量基因测序图像。

图5是本发明第一具体实施例中根据3种对比不同测序图像清晰度方法所得结果经归一化处理后所做的对比图。

图6是本发明的比较多个测序图像清晰度的系统的简单结构示意图。

图7a至7d是本发明第二具体实施例中的高通量基因测序图像。

图8a至8e是本发明第三具体实施例中的高通量基因测序图像。

图9是图1a所示的测序图像的来源。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

本发明提出第一实施例，一种比较多个测序图像清晰度的方法，包括以下步骤：

StepA.对多个测序图像分别进行谐波分析处理，得多个与测序图像相对应的谐波分析处理图；

StepB.根据所述与测序图像相对应的谐波分析处理图的纹理特征分别拟合出多个与测序图像相对应的闭合曲线；

StepC.确定所述与测序图像相对应的闭合曲线所代表的图形与其相应的测序图像之间的面积比值；

StepD.根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像。

多个用于比较的测序图像经本实施例的方法处理后，通过StepC得到的面积比值确定其中哪个测序图像的清晰度最高，面积比值越大清晰度越高。本实施例的方法与灰度评价梯度法、拉普拉斯算子和等现有技术相比，得到的值不像灰度评价梯度法、拉普拉斯算子和等现有技术那样数量级很大，难以制定出合理的评价标量，且其受图像中发光点密度的影响更小，另外通过面积比值消除了因测序图像大小差异而引起的闭合曲线所对应的图形的面积变化，因此，本实施例的方法较现有技术中的清晰度比较方法具有量化精度高的优点。对此，本发明将在后续的实施例中进一步阐述。

需要说明的是：

所述测序图像为高通量基因测序图像或高通量基因测序图像中的一部分，不包括通过第一代Sanger测序法得到的测序图像或其中的一部分，尤指通过连接测序法或合成测序法获得的测序图像。高通量基因测序图像中存在大量的发光点，在聚焦准确测序图像清晰的情况下发光点的形状为圆形或似圆形，在聚焦不准测序图像不清晰的情况下发光点可能为圆形、似圆形或椭圆形。

所述高通量基因测序图像可以是黑白的灰度图也可以是彩色图像，当高通量基因测序图像为彩色图像时需先将高通量基因测序图像转换为黑白的灰度图。

当测序图像为高通量基因测序图像中的一部分时，可减少本方法的运算量，提高效率。

因为高通量基因测序图像中的物体都处于同一平面，且拍图组件中存在凸透镜或凸透镜组，所以，在拍摄测序图像时，将聚焦点设在该次采图范围的中间能够得到整体上清晰度最好的测序图像。

优选的，作为所述多个测序图像的来源，各高通量基因测序图像为通过同一测序系统在相同采图条件下获得的图像；更优选的，所述各高通量基因测序图像为通过同一测序系统在相同采图条件下针对相同或类似测序反应体系获得的图像。

为了使测序图像更有代表性，更能准确反应高通量基因测序图像的清晰度，本发明采用中心区域法或五区域法对高通量基因测序图像进行截取进而获得测序图像，优选中心区域法。例如，图1a即是利用中心区域法对图9进行截去后获得的测序图像；其中，图9的长、宽均为2048像素，图1a的长、宽均为512像素。

其中，所述谐波分析处理可为傅立叶变换、小波变换或沃尔什变换。

优选的，所述谐波分析处理为傅立叶变换。此时，所述闭合曲线为椭圆曲线或圆形曲线，所述闭合曲线对应的图形为椭圆或圆。

以下主要以谐波分析处理为傅立叶变换为例，对本发明进行进一步的阐述说明。

针对StepA，本发明提出第一实施例，包括以下步骤：

StepA1.根据式Ⅰ对多个测序图像分别进行预处理，得多个预处理后的图像；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；

StepA2.根据式Ⅱ对多个预处理后的图像分别进行傅立叶变换，使之由空域转换至频域，得多个傅立叶变换图；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}g(x,y)e^{-j2\mathrm{\π}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

本实施例通过上述预处理，使得测序图像在进行二维离散傅立叶变换后得到的傅立叶变换图更易拟合出椭圆曲线或圆形曲线。

针对StepA，还可通过第二实施例实现，包括以下步骤：

StepA1’.根据式Ⅱ分别对测序图像进行傅立叶变换，使之由空域转换至频域，得多个前期傅立叶变换图；

StepA2’.根据式Ⅴ分别对前期傅立叶变换图进行变换，从而最终得到后处理后的多个傅立叶变换图。

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}g(x,y)e^{-j2\mathrm{\π}(\frac{ux}{M}+\frac{vy}{N})}---II;$

$\left\{\begin{array}{c}0\&le;u\&le;\frac{M}{2},0\&le;v\&le;\frac{N}{2};G(u,v)=G(u+\frac{M}{2},v+\frac{N}{2})\\ \frac{M}{2}<u\&le;M,0\&le;v\&le;\frac{N}{2};G(u,v)=G(u-\frac{M}{2},v+\frac{N}{2})\\ 0\&le;u\&le;\frac{M}{2},\frac{N}{2}<v\&le;N;G(u,v)=G(u+\frac{M}{2},v-\frac{N}{2})\\ \frac{M}{2}<u\&le;M,\frac{N}{2}<v\&le;N;G(u,v)=G(u-\frac{M}{2},v-\frac{N}{2})\end{array}\right.---V;$

其中，所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

以图1a作为测序图像，本实施例获得的前期傅立叶变换图如图1b所示，但是其经StepA2’后获得的傅立叶变换图如图1c所示，与以图1a作为测序图像，经第一实施例获得的傅立叶变换图相同。

需要说明的是：第一实施例中的运算更少，效率更高。

针对StepB，以下将通过多个实施例来进行进一步阐述。

本发明提出第三实施例，包括以下步骤：

StepB1.利用半径为r的圆形结构元对经StepA得到的多个傅立叶变换图分别进行闭运算；

StepB2.利用半径为R的圆形结构元对经StepB1得到的结果图像分别进行开运算；

StepB3.利用边缘检测方法对经StepB2得到的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是：

在第三实施例中，R＞r；如图1d所示，通过对傅立叶变换图的闭运算可将亮点区域中的暗点细节去掉，使得椭圆或圆的轮廓得到明显的表现；因为闭运算后的结果图像中的椭圆或圆外存在较多的斑点，如果直接在此基础上进行边缘信息提取，干扰较大；因此，本实施例中通过开运算将小区域内的亮点细节去掉，如图1e所示，并同时达到平滑椭圆或圆边缘的效果，使得StepB3中的边缘提取受到的干扰更小，结果更准确。以图1a作为测序图像，按照第一实施例和第三实施例的方法处理后，得到的结果图如图1f所示。

为了进一步减少对后续边缘提取的干扰，本发明提出第四实施例，包括以下步骤：

StepB1.利用半径为r的圆形结构元对经StepA得到的多个傅立叶变换图分别进行闭运算；

StepB1’.将经StepB1得到的结果图像中的灰度级小于等于a的值重置为0；

StepB2.利用半径为R的圆形结构元对经StepB1’得到的结果图像分别进行开运算；

StepB3.利用边缘检测方法对经StepB2所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是，在第四实施例中，10＜a＜15，R＞r；与第三实施例相比，第四实施例通过将闭运算得到的结果图像中灰度级较低的值重置为0，然后再进行开运算的方式，使得StepB2得到的结果图像中的椭圆或圆更为明显，有利于StepB3中的边缘提取和曲线拟合。

为了降低图像处理过程中的运算复杂度，本发明提出第五实施例，包括以下步骤：

StepB0.分别对经StepA所得的多个傅立叶变换图进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

StepB1.利用半径为r的圆形结构元对经StepB0得到的多个傅立叶同态变换图分别进行闭运算；

StepB2.利用半径为R的圆形结构元对经StepB1得到的结果图像分别进行开运算；

StepB3.利用边缘检测方法对经StepB2所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是，在第五实施例中，R＞r。因为测序图像实为一种数字图像，其灰度级非常大，这使得计算出的频率谱的值级数也较大，若直接以计算出的频率谱的值进行后续运算，计算的复杂度较高。与第三、第四实施例相比，第五实施例在进行闭运算前，先对经StepA所得的多个傅立叶变换图分别进行了灰度级压缩，这降低了运算的复杂度，提高效率，并且有利于后续频率谱的平滑，使得边缘更易突出从而拟合出椭圆曲线或圆形曲线，且灰度级压缩对后续的计算精度微乎其微，可忽略不记。

本发明在第五实施例的基础上进一步提出了第六实施例，包括以下步骤：

StepB0.分别对经StepA所得的多个傅立叶变换图进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

StepB0’.将经StepB0所得的傅立叶同态变换图中灰度级小于等于a的值重置为0；

StepB1.利用半径为r的圆形结构元对经StepB0’所得的结果图像分别进行闭运算；

StepB2.利用半径为R的圆形结构元对经StepB1所得的结果图像分别进行开运算；

StepB3.利用边缘检测方法对经StepB2所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是，所述10＜a＜15，R＞r。第六实施例在灰度级压缩之后，闭运算之前，将傅立叶同态变换图中灰度级小于等于a的值重置为0，使StepB2得到的结果图像中的椭圆或圆更为明显，有利于StepB3中的边缘提取和曲线拟合。此外，该重置步骤还可以放在闭运算之后，开运算之前，此方案与第六实施例一样，能够使StepB2得到的结果图像中的椭圆或圆更为明显，有利于StepB3中的边缘提取和曲线拟合。

针对StepB0，本发明提出第七实施例，在本实施例中，所述StepB0为：根据式Ⅳ对多个与测序图像相对应的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。

需要说明的是，根据式Ⅳ对傅立叶变换图进行灰度级压缩后的数值为正整数，本步骤有利于后续频率谱的平滑，使边缘更易突出，更易拟合出准确的椭圆曲线或圆形曲线。

针对StepB3，所述边缘检测方法可采用Sobel边缘检测算法、canny边缘检测算法、robert边缘检测算法、prewite边缘检测算法或kirsch边缘检测算法。

优选的，所述边缘检测方法为Sobel边缘检测算法，所述Sobel边缘检测算法是提取水平和竖直两个方向的纹理梯度信息的算法，符合人眼对图像观测时的纹理提取过程，这使得通过本发明方法得到的结果与人眼的主观感觉一致。

优选的，所述R等于1.5r、2r、2.5r或3r；更优选的，所述R等于2r。

StepB所拟合出的椭圆曲线均可用式Ⅲ所示的椭圆曲线的一般方程表示，Ax²+Bxy+Cy²+Dx+Ey+F＝0Ⅲ。

StepB所拟合出的圆形曲线均可用下式所示的圆形曲线的一般方程表示：

x²+y²+Dx+Ey+F＝0(D²+E²-4F＞0)。

由上可知，实际上椭圆曲线和圆形曲线对应的一般方程均可用式Ⅲ表示，只不过所述曲线为圆形曲线时，式Ⅲ中的A和C等于1，B等于0。即，步骤B进行曲线拟合时，均可按照式Ⅲ进行曲线拟合。

针对StepC，本发明提出第八实施例，包括以下步骤：

StepC1.分别根据式Ⅲ计算出各椭圆曲线或圆形曲线所对应的椭圆或圆的面积；

StepC2.确定各椭圆曲线或圆形曲线所对应的椭圆或圆和相应的测序图像之间的面积比值。

其中，当所拟合出的的曲线为椭圆曲线时，StepC1可通过第九实施例进行，包括以下步骤：

StepC11.按照式Ⅵ计算经StepB得到的椭圆曲线所对应的椭圆的旋转角度α；

StepC12.将式Ⅶ带入经步骤B得到的椭圆曲线的一般方程式Ⅲ，得旋转前的椭圆的一般方程式Ⅷ；

StepC13.根据式Ⅸ计算得出式Ⅷ中A’、B’、C’、D’、E’、F’的值；

StepC14.根据式Ⅹ计算得出椭圆曲线所对应的椭圆的长半轴m的平方和短半轴n的平方；

StepC15.根据式Ⅺ计算得出椭圆曲线所代表的椭圆的面积S。

$\mathrm{\&alpha;}=\frac{1}{2}atan\&lsqb;B/(A-C)\&rsqb;---VI;$

$\left\{\begin{array}{c}x=x^{\&prime;}cos\mathrm{\&alpha;}-y^{\&prime;}sin\mathrm{\&alpha;}\\ y=y^{\&prime;}cos\mathrm{\&alpha;}+x^{\&prime;}sin\mathrm{\&alpha;}\end{array}\right.---VII;$

A'x'²+B'x'y'+C'y'²+D'x'+E'y'+F'＝0Ⅷ；

$\left\{\begin{array}{c}{A}^{\&prime;}=\frac{1}{2}\&lsqb;A+C+(A-C)cos2\mathrm{\&alpha;}+Bsin2\mathrm{\&alpha;}\&rsqb;\\ B^{\&prime;}=Bcos2\mathrm{\&alpha;}-(A-C)sin2\mathrm{\&alpha;}\\ C^{\&prime;}=\frac{1}{2}\&lsqb;A+C-(A-C)\cos 2\mathrm{\&alpha;}-Bsin2\mathrm{\&alpha;}\&rsqb;\\ D^{\&prime;}=D\cos \mathrm{\&alpha;}+E\sin \mathrm{\&alpha;}\\ E^{\&prime;}=Ecos\mathrm{\&alpha;}-Dsin\mathrm{\&alpha;}\\ F^{\&prime;}=F\end{array}\right.---IX;$

$\{\begin{array}{c}{a}^2=\{{\&lsqb;E^{\&prime;}-B^{\&prime;}{D}^{\&prime;}/\left(2A^{\&prime;}\right)\&rsqb;}^2/\&lsqb;4C^{\&prime;}-B^{\&prime;2}/A^{\&prime;}\&rsqb;-F^{\&prime;}+D^{\&prime;2}/\left(4A^{\&prime;}\right)\}/A^{\&prime;}\\ b^2=\{{\&lsqb;E^{\&prime;}-B^{\&prime;}{D}^{\&prime;}/\left(2A^{\&prime;}\right)\&rsqb;}^2/\&lsqb;4C^{\&prime;}-B^{\&prime;2}/A^{\&prime;}\&rsqb;-F^{\&prime;}+D^{\&prime;2}/\left(4A^{\&prime;}\right)\}/\&lsqb;C^{\&prime;}-B^{\&prime;2}/\left(4A^{\&prime;}\right)\&rsqb;\end{array}---X;$

S＝πmnⅪ。

按第九实施例的方法，可计算出图1f中的椭圆的长短半轴分别为139.8、134.7，椭圆面积与其对应的测序图像的面积比值为：22.56％。

以下将通过第一具体实施例来进一步对本发明进行详细说明。

以图2a-f，图3a-f、图4a-e所示的高通量基因测序图像为基础，这些高通量基因测序图像是在一个测序反应中，由同一高通量基因测序仪针对测序反应小室中的同一采图位置，在不同的聚焦位置，按照焦距逐步变大的方式采集的图像，这些图像的长和宽均为512像素。

以上述17个高通量基因测序图像为测序图像来进行后续的运算。

一、实验组。

1、根据式Ⅰ对上述选取的17个测序图像分别进行预处理，得17个预处理后的图像。

2、根据式Ⅱ对上述17个预处理后的图像分别进行二维离散傅立叶变换，使之由空域转换至频域，得17个傅立叶变换图。

3、根据式Ⅳ对17个傅立叶变换图分别进行灰度级压缩，得17个与测序图像相对应的傅立叶同态变换图。

4、将17个傅立叶同态变换图中灰度级小于等于13的值重置为0。

5、利用半径为5的圆形结构元对步骤4所得的结果图像分别进行闭运算。

6、利用半径为10的圆形结构元对步骤5所得的结果图像分别进行开运算。

7、利用Sobel边缘检测算法对步骤6所得的结果图像分别进行边缘提取，进而进行曲线拟合，得17个与测序图像相对应的椭圆曲线。

8、参照第九实施例所述的方法，计算上述17个与测序图像相对应的椭圆曲线所代表的椭圆的面积与测序图像的面积比值。

二、对照组。

1、根据下式分别计算图图2a-f，图3a-f、图4a-e的灰度平均梯度值

$GMG=\frac{1}{(M-1)(N-1)}\underset{i=1}{\overset{M-1}{\&Sigma;}}\underset{j=1}{\overset{N-1}{\&Sigma;}}\sqrt{\frac{{\&lsqb;g(i,j+1)-g(i,j)\&rsqb;}^2+{\&lsqb;g(i+1,j)-g(i,j)\&rsqb;}^2}{2}}.$

2、根据下式分别计算图7-23的的拉普拉斯算子和

$LS=\frac{1}{(M-2)(N-2)}\underset{i=2}{\overset{M-1}{\&Sigma;}}\underset{j=2}{\overset{N-1}{\&Sigma;}}|9g(i,j)-\underset{p=-1}{\overset{1}{\&Sigma;}}\underset{q=-1}{\overset{1}{\&Sigma;}}g(i+p,j+q)|.$

上述3种方法获得的值的结果如表1所示，对上述3种方法获得的值进行归一化处理，结果如表2所示，并根据表2作图得图5。

表1.第一具体实施例的计算结果

表2.对表1中的数据分别进行归一化处理后的结果

由表1、2，图5可知，上述3种方法得到的值中最大的值所对应的测序图像均为图3b，即，通过上述方法确定的清晰度最高的图像是相同的。但是本发明的方法对清晰度最佳位置所对应的焦距附近的图像的评分较对照组高，这与通过人眼对图像的清晰度评价是一致的，即清晰度最佳位置所对应的焦距是一个范围值。不过本发明的方法对距清晰度最佳位置所对应的焦距较远的图像的评分较对照组低，即差异度更大。综上，本发明的方法在对模糊程度不同的测序图像进行质量评价时量化精度更高。

另外，为了证明本发明的方法受图像中发光点密度的影响更小，本申请提出第二具体实施例进行说明。

以图7a至7d为基础，分别按照上述具体实施例的方法计算相应的值，其中，图7a至7d是由同一高通量基因测序仪，以相同的焦距，针对样品密度不同的测序反应小室采图得到的图像，图像的大小完全相同，图像的长宽均为512个像素，另外，图像中发光点的大小也相同，结果如表3所示。

表3.不同样品密度的测序图像的计算结果

由表3可知，随着样品密度的增加，即发光点密度的增加，灰度平均梯度值法和拉普拉斯算子和法得到的值的变动趋势明显较本发明的方法大很多，即，本发明的方法受图像中发光点密度的影响显著的小于灰度平均梯度值法和拉普拉斯算子和法。

另外，本申请提出了第三具体实施例对本发明进行进一步的说明。

以图8a至8e为基础，分别按照上述具体实施例的方法计算相应的值。其中，图8a至8e是由与第一实施例相同的高通量基因测序仪，在放大倍数为第一具体实施例中的2倍的情况下对样品密度不同的测序反应小室采图得到的图像，结果如表4所示。

表4.第三具体实施例的计算结果

需要说明的，第三具体实施例中的发光点的大小是第一具体实施例中的2倍大，且图8a至8e的清晰度可由肉眼直接判断是不断变好的。由表4的结果可知，虽然，现有技术中的灰度平均梯度值法、拉普拉斯算子和法与本发明的方法一样，均能确定图8e是清晰度最佳的图像，但是，本发明的方法对图8a至8e的区分度更好，且与实际的观感更为契合。例如，图8d和8e的发光点密度是一致的，但是现有技术中的灰度平均梯度值法、拉普拉斯算子和法与肉眼的观感有区别，而本发明的方法则是一致的；图8b和8c的清晰度差别明显，但是，因为现有技术中的灰度平均梯度值法、拉普拉斯算子和法受测序图像中的发光点密度的影响大，所以，图8b和8c经灰度平均梯度值法、拉普拉斯算子和法计算得出的值相差不大；而利用本发明的方法计算得到的面积比值则差距明显。

综上，本发明的方法与现有技术相比，受测序图像中发光点的密度影响更小，量化精度更高，结果更为准确。

需要说明的是，当本发明的方法中的谐波分析处理采用傅立叶变换时，得到的闭合曲线的形状主要与测序图像中的发光点形状相关。当采图范围中的发光点为基板上的发光簇或规则的圆形荧光磁珠时，所述谐波分析处理优选为傅立叶变换。

当然，针对与上述的相同或不同的发光点，本发明的方法中的谐波分析处理还可采用小波变换、沃尔什变换等谐波分析处理方法；得到的闭合曲线还可以为三角形、方形、菱形、多边形等闭合曲线，尤其是对称的闭合曲线。采用小波变换、沃尔什变换与采用傅立叶变换一样，能够提高本发明的方法的量化精度。

图6为本发明的比较多个测序图像清晰度的系统的简单结构示意图，包括：谐波分析处理单元1、闭合曲线拟合单元2、面积比值计算单元3和清晰度判定单元4；

所述谐波分析处理单元1用于对测序图像进行谐波分析处理得谐波分析处理图；

所述闭合曲线拟合单元2用于根据所述谐波分析处理图的纹理特征拟合出闭合曲线；

所述面积比值计算单元3用于确定所述闭合曲线所代表的图形和相应的测序图像之间的面积比值；

所述清晰度判定单元4用于根据面积比值的大小确定所述测序图像的清晰度。

现有技术中利用灰度评价梯度法、拉普拉斯算子和等方法的系统那样，得到的数量级很大，难以制定出合理的评价标量。本系统可用于评价测序图像的清晰度，面积比值越大清晰度越高，且其受图像中发光点密度的影响更小，另外通过面积比值消除了因输入图像大小差异而引起的闭合曲线所代表的图形的面积变化，因此，本系统与现有技术相比，量化精度更高。

其中，所述谐波分析处理单元用于对测序图形进行谐波分析处理得谐波分析处理图可具体为：所述谐波分析处理单元用于对测序图像进行傅立叶变换得傅立叶变换图；或所述谐波分析处理单元用于对测序图像进行小波变换得小波变换图；或所述谐波分析处理单元用于对测序图像进行沃尔什变换得沃尔什变换图。

优选的，所述谐波分析处理单元用于对测序图像进行傅立叶变换得傅立叶变换图；所述闭合曲线拟合单元用于根据所述傅立叶变换图的纹理特征拟合出椭圆曲线或圆形曲线；所述面积比值计算单元用于确定所述椭圆曲线或圆形曲线所对应的椭圆或圆和相应的测序图像之间的面积比值。

以下主要以谐波分析处理单元用于对测序图像进行傅立叶变换得傅立叶变换图为例，对本发明进行进一步的阐述说明。

进一步的，所述系统还可包括预处理单元；所述预处理单元用于根据式Ⅰ对测序图像进行预处理，得预处理后的图像；所述谐波分析处理单元用于根据式Ⅱ对预处理后的图像进行傅立叶变换，使之由空域转换至频域，得傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\&Sigma;}}\underset{y=0}{\overset{N-1}{\&Sigma;}}g(x,y)e^{-j2\mathrm{\&pi;}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

针对闭合曲线拟合单元，本发明提出一实施例，包括：闭运算模块、开运算模块和边缘提取曲线拟合模块；所述闭运算模块用于利用半径为r的圆形结构元对傅立叶变换图进行闭运算；所述开运算模块利用半径为R的圆形结构元对闭运算模块处理后所得的结果图像进行开运算；所述边缘提取曲线拟合模块用于利用边缘检测方法对开运算模块处理后所得的结果图像进行边缘提取，进而进行曲线拟合，得椭圆曲线或圆形曲线。

需要说明的是，本实施例中，R＞r；通过闭运算模块可将亮点区域中的暗点细节去掉，使得椭圆或圆的轮廓得到明显的表现；因为闭运算后的结果图像中的椭圆或圆外存在较多的斑点，如果边缘提取曲线拟合模块直接在此基础上进行边缘信息提取，干扰较大；因此，本实施例中通过开运算模块将小区域内的亮点细节去掉，并同时达到平滑椭圆或圆边缘的效果，使得边缘提取曲线拟合模块的边缘提取受到的干扰更小，结果更准确。

针对闭合曲线拟合单元，本发明提出另一实施例，包括：闭运算模块、重置模块、开运算模块和边缘提取曲线拟合模块；所述闭运算模块用于利用半径为r的圆形结构元对傅立叶变换图进行闭运算；所述重置模块用于将闭运算模块处理后所得的结果图像中的灰度级小于等于a的值重置为0；所述开运算模块利用半径为R的圆形结构元对重置模块处理后所得的结果图像进行开运算；所述边缘提取曲线拟合模块用于利用边缘检测方法对开运算模块处理后所得的结果图像进行边缘提取，进而进行曲线拟合，得椭圆曲线或圆形曲线。

需要说明的是，在本实施例中，10＜a＜15，R＞r；与上一实施例相比，本实施例通过重置模块将闭运算模块处理后得到的结果图像中灰度级较低的值重置为0，然后再通过开运算模块对重置模块处理后所得的结果图像进行开运算，这使得开运算模块处理后得到的结果图像中的椭圆或圆更为明显，即，减少了对后续边缘提取的干扰，有利于边缘提取曲线拟合模块的边缘提取和曲线拟合。

针对闭合曲线拟合单元，本发明还提出了另一实施例，包括灰度级压缩模块、闭运算模块、开运算模块和边缘提取曲线拟合模块；所述灰度级压缩模块用于对傅立叶变换图进行灰度级压缩，得傅立叶同态变换图；所述闭运算模块用于利用半径为r的圆形结构元对灰度级压缩模块处理后所得的结果图像进行闭运算；所述开运算模块用于利用半径为R的圆形结构元对闭运算模块处理后所得的结果图像进行开运算；所述边缘提取曲线拟合模块用于利用边缘检测方法对开运算模块处理后所得的结果图像进行边缘提取，进而进行曲线拟合，得椭圆曲线或圆形曲线。

需要说明的是，在本实施例中，R＞r。因为测序图像实为一种数字图像，其灰度级非常大，这使得计算出的频率谱的值级数也较大，若直接以谐波分析处理单元计算出的频率谱的值进行后续运算，计算的复杂度较高。与上述两个实施例相比，本实施例在进行闭运算前，先利用灰度级压缩模块对傅立叶变换图进行灰度级压缩，这降低了运算的复杂度，提高效率，并且有利于后续频率谱的平滑，使得边缘更易突出从而拟合出椭圆曲线或圆形曲线，且灰度级压缩对后续的计算精度微乎其微，可忽略不记。

针对闭合曲线拟合单元，本发明还提出了另一实施例，包括灰度级压缩模块、重置模块、闭运算模块、开运算模块和边缘提取曲线拟合模块；所述灰度级压缩模块用于对傅立叶变换图进行灰度级压缩，得傅立叶同态变换图；所述重置模块用于将傅立叶同态变换图中灰度级小于等于a的值重置为0；所述闭运算模块用于利用半径为r的圆形结构元对重置模块处理后所得的结果图像进行闭运算；所述开运算模块用于利用半径为R的圆形结构元对闭运算模块处理后所得的结果图像进行开运算；所述边缘提取曲线拟合模块用于利用边缘检测方法对开运算模块处理后所得的结果图像进行边缘提取，进而进行曲线拟合，得椭圆曲线或圆形曲线。

需要说明的是，本实施例中，所述10＜a＜15，R＞r。本实施例通过重置模块在进一步利用闭运算模块处理前，对灰度级压缩模块处理后得到的结果图像中灰度级小于等于a的值重置为0，使得闭运算模块处理前的图像中的椭圆或圆更为明显，有利于边缘提取曲线拟合模块的边缘提取和曲线拟合。另外，若通过重置模块在进一步利用开运算模块处理前，对闭运算模块处理后得到的结果图像中灰度级小于等于a的值重置为0，同样能够使开运算模块得到的结果图像中的椭圆或圆更为明显，有利于边缘提取曲线拟合模块的边缘提取和曲线拟合。

针对上述实施例中的灰度级压缩模块，本发明提出另一实施例，所述灰度级压缩模块根据式Ⅳ对傅立叶变换图进行灰度级压缩，得傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。

需要说明的是，本实施例中灰度级压缩模块根据式Ⅳ对傅立叶变换图进行灰度级压缩后的数值为正整数，有利于后续频率谱的平滑，使边缘更易突出，更易拟合出准确的椭圆曲线或圆形曲线。

针对所述边缘提取曲线拟合模块，需要指出的是，该模块可采用Sobel边缘检测算法、canny边缘检测算法、robert边缘检测算法、prewite边缘检测算法或kirsch边缘检测算法对开运算模块处理后得到的结果图像进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

优选的，所述边缘提取曲线拟合模块采用Sobel边缘检测算法对开运算模块处理后得到的结果图像进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。所述Sobel边缘检测算法是提取水平和竖直两个方向的纹理梯度信息的算法，符合人眼对图像观测时的纹理提取过程，这使得通过本发明的系统得到的结果与人眼的主观感觉一致。

优选的，所述面积比值计算单元包括面积计算模块、比值计算模块；所述面积计算模块用于确定闭合曲线所对应的图形的面积；所述比值计算模块用于确定所述闭合曲线所对应的图形的面积和相应的测序图像的面积之间的比值。

优选的，所述R等于1.5r、2r、2.5r或3r；更优选的，所述R等于2r。

需要说明的是，本发明的系统中的谐波分析处理单元用于对测序图像进行小波变换得小波变换图，或用于对测序图像进行沃尔什变换得沃尔什变换图，同样能够实现本发明的目的，量化精度更高。

在上述的比较多个测序图像清晰度的方法的基础上，本发明又提出一种自动聚焦控制方法，包括以下步骤：

A.调焦模块多次调节成像装置与采图位置之间的距离，并在每个距离均发出采集图像的指令给采图模块；

B.采图模块在接收到采集图像的指令后，利用成像装置采集测序图像，进而得到针对同一采图位置的多个测序图像；

C.图像处理模块在接收到经步骤B得到的多个测序图像后，根据上述的比较多个测序图像清晰度的方法，确定清晰度最高的测序图像，进而获得最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块；

D.调焦模块根据接收到的指令，控制成像装置与采图位置之间的距离为最佳聚焦距离。

上述方法与现有技术相比，利用前述的比较多个测序图像清晰度的方法确定最佳聚焦距离，能够准确高效的定位自动聚焦位置，采图模块在此条件下采集到的测序图像更为清晰，图像质量更高。

针对所述步骤C，本发明提出一实施例，所述图像处理模块包括谐波分析处理单元、闭合曲线拟合单元、面积比值计算单元和最佳聚焦距离判定单元。此时，所述步骤C包括以下步骤：

C1.谐波分析处理单元对经步骤B所得的多个测序图像分别进行谐波分析处理，得多个与测序图像相对应的谐波分析处理图；

C2.闭合曲线拟合单元分别根据所述谐波分析处理图的纹理特征拟合出多个与测序图像相对应的闭合曲线；

C3.面积比值计算单元确定所述闭合曲线所对应的图形和其相应的测序图像之间的面积比值；

C4.所述最佳聚焦距离判定单元根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像，进而确定最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块。

其中，所述步骤C1中的谐波分析处理可为傅立叶变换、小波变换或沃尔什变换。

优选的，所述步骤C包括以下步骤：

C1.谐波分析处理单元对经步骤B所得的多个测序图像分别进行傅立叶变换，得多个与测序图像相对应的傅立叶变换图；

C2.闭合曲线拟合单元分别根据所述傅立叶变换图的纹理特征拟合出多个与测序图像相对应的椭圆曲线或圆形曲线；

C3.面积比值确定单元确定所述椭圆曲线或圆形曲线所对应的椭圆或圆和其相应的测序图像之间的面积比值；

C4.所述最佳聚焦距离判定单元根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像，进而确定最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块。

以下主要以步骤C1中的谐波分析处理为傅立叶变换为例，对本发明进行进一步的阐述说明。

针对所述步骤C1，本发明提出另一实施例，所述图像处理模块还包括预处理单元；所述步骤C1包括以下步骤：

C11.所述预处理单元根据式Ⅰ对经步骤B所得的多个测序图像分别进行预处理，得多个预处理后的图像；

C12.所述谐波分析处理单元根据式Ⅱ对多个预处理后的图像分别进行傅立叶变换，使之由空域转换至频域，得多个与测序图像相对应的傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\&Sigma;}}\underset{y=0}{\overset{N-1}{\&Sigma;}}g(x,y)e^{-j2\mathrm{\&pi;}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

针对所述步骤C1，本发明提出另一实施例，包括以下步骤：

StepC11’.谐波分析处理单元根据式Ⅱ分别对测序图像进行傅立叶变换，使之由空域转换至频域，得多个前期傅立叶变换图；

StepC12’.根据式Ⅴ分别对前期傅立叶变换图进行变换，从而最终得到后处理后的多个傅立叶变换图。

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\&Sigma;}}\underset{y=0}{\overset{N-1}{\&Sigma;}}g(x,y)e^{-j2\mathrm{\&pi;}(\frac{ux}{M}+\frac{vy}{N})}---II;$

$\left\{\begin{array}{c}0\&le;u\&le;\frac{M}{2},0\&le;v\&le;\frac{N}{2};G(u,v)=G(u+\frac{M}{2},v+\frac{N}{2})\\ \frac{M}{2}<u\&le;M,0\&le;v\&le;\frac{N}{2};G(u,v)=G(u-\frac{M}{2},v+\frac{N}{2})\\ 0\&le;u\&le;\frac{M}{2},\frac{N}{2}<v\&le;N;G(u,v)=G(u+\frac{M}{2},v-\frac{N}{2})\\ \frac{M}{2}<u\&le;M,\frac{N}{2}<v\&le;N;G(u,v)=G(u-\frac{M}{2},v-\frac{N}{2})\end{array}\right.---V;$

其中，所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

需要说明的是：本实施例和上一实施例能实现相同的技术效果，但是上一实施例中的运算更少，效率更高。

针对所述StepC2，以下将通过多个实施例来进行进一步阐述。

本发明提出一实施例，包括以下步骤：

StepC21.利用半径为r的圆形结构元对经StepC1得到的多个傅立叶变换图分别进行闭运算；

StepC22.利用半径为R的圆形结构元对经StepC21得到的结果图像分别进行开运算；

StepC23.利用边缘检测方法对经StepC22得到的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是：

在本实施例中，R＞r；通过对傅立叶变换图的闭运算可将亮点区域中的暗点细节去掉，使得椭圆或圆的轮廓得到明显的表现；此外，因为闭运算后的结果图像中的椭圆或圆外存在较多的斑点，如果直接在此基础上进行边缘信息提取，难度较大，因此，本实施例中通过开运算将小区域内的亮点细节去掉，并同时达到平滑椭圆或圆边缘的效果，使得StepC23中的边缘提取受到的干扰更小，结果更准确。

为了进一步减少对后续边缘提取的干扰，本发明提出另一实施例，包括以下步骤：

StepC21.利用半径为r的圆形结构元对经StepC1得到的多个傅立叶变换图分别进行闭运算；

StepC21’.将经StepC21得到的结果图像中的灰度级小于等于a的值重置为0；

StepC22.利用半径为R的圆形结构元对经StepC21’得到的结果图像分别进行开运算；

StepC23.利用边缘检测方法对经StepC22所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是，在本实施例中，10＜a＜15，R＞r；与上一实施例相比，本实施例通过将闭运算得到的结果图像中灰度级较低的值重置为0，然后再进行开运算的方式，使得StepC22得到的结果图像中的椭圆或圆更为明显，有利于StepC23中的边缘提取和曲线拟合。

为了降低图像处理过程中的运算复杂度，本发明提出另一实施例，包括以下步骤：

StepC20.对经StepC1所得的多个傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

StepC21.利用半径为r的圆形结构元对经StepC20得到的多个傅立叶同态变换图分别进行闭运算；

StepC22.利用半径为R的圆形结构元对经StepC21得到的结果图像分别进行开运算；

StepC23.利用边缘检测方法对经StepC22所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是，在本实施例中，R＞r。因为测序图像实为一种数字图像，其灰度级非常大，这使得计算出的频率谱的值级数也较大，若直接以计算出的频率谱的值进行后续运算，计算的复杂度较高。与上两个实施例相比，本实施例在进行闭运算前，先对StepC1所得的多个傅立叶变换图分别进行了灰度级压缩，这降低了运算的复杂度，提高效率，并且有利于后续频率谱的平滑，使得边缘更易突出从而拟合出椭圆曲线或圆形曲线，且灰度级压缩对后续的计算精度微乎其微，可忽略不记。

本发明在上一实施例的基础上进一步提出了另一实施例，包括以下步骤：

StepC20.对经StepC1所得的多个傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

StepC20’.将经StepC20所得的傅立叶同态变换图中灰度级小于等于a的值重置为0；

StepC21.利用半径为r的圆形结构元对经StepC20’所得的结果图像分别进行闭运算；

StepC22.利用半径为R的圆形结构元对经StepC21所得的结果图像分别进行开运算；

StepC23.利用边缘检测方法对经StepC22所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线。

需要说明的是，所述10＜a＜15，R＞r。本实施例在灰度级压缩之后，闭运算之前，将傅立叶同态变换图中灰度级小于等于a的值重置为0，使经StepC22得到的结果图像中的椭圆或圆更为明显，有利于StepC23中的边缘提取和曲线拟合。此外，该重置步骤还可以放在闭运算之后，开运算之前，此方案与本实施例一样，能够使经StepC22得到的结果图像中的椭圆或圆更为明显，有利于StepC23中的边缘提取和曲线拟合。

针对StepC20，本发明提出另一实施例，在本实施例中，所述StepC20为：根据式Ⅳ对多个与测序图像相对应的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。

需要说明的是，根据式Ⅳ对傅立叶变换图进行灰度级压缩后的数值为正整数，本步骤有利于后续频率谱的平滑，使边缘更易突出，更易拟合出准确的椭圆曲线或圆形曲线。

针对StepC23，所述边缘检测方法可采用Sobel边缘检测算法、canny边缘检测算法、robert边缘检测算法、prewite边缘检测算法或kirsch边缘检测算法。

优选的，所述边缘检测方法为Sobel边缘检测算法，所述Sobel边缘检测算法是提取水平和竖直两个方向的纹理梯度信息的算法，符合人眼对图像观测时的纹理提取过程，这使得通过本发明方法得到的结果与人眼的主观感觉一致。

优选的，所述R等于1.5r、2r、2.5r或3r；更优选的，所述R等于2r。

需要说明的是，当本发明的方法中的谐波分析处理采用傅立叶变换时，得到的闭合曲线的形状主要与测序图像中的发光点形状相关。当采图范围中的发光点为基板上的发光簇或规则的圆形荧光磁珠时，所述谐波分析处理优选为傅立叶变换。

当然，针对与上述的相同或不同的发光点，本发明的方法中的谐波分析处理还可采用小波变换、沃尔什变换等谐波分析处理方法；得到的闭合曲线还可以为三角形、方形、菱形、多边形等闭合曲线，尤其是对称的闭合曲线。采用小波变换、沃尔什变换与采用傅立叶变换一样，能够精确定位自动聚焦位置。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (17)

1.一种比较多个测序图像清晰度的方法，其特征在于，包括以下步骤：

A.对多个测序图像分别进行谐波分析处理，得多个与测序图像相对应的谐波分析处理图；

B.分别根据所述多个与测序图像相对应的谐波分析处理图的纹理特征拟合出多个与测序图像相对应的闭合曲线；

C.确定所述闭合曲线对应的图形和相应的测序图像之间的面积比值；

D.根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像；

所述测序图像为高通量基因测序图像或高通量基因测序图像中的一部分。

2.根据权利要求1所述的比较多个测序图像清晰度的方法，其特征在于：所述谐波分析处理为傅立叶变换；所述闭合曲线为椭圆曲线或圆形曲线；所述闭合曲线对应的图形为椭圆或圆。

3.根据权利要求2所述的比较多个测序图像清晰度的方法，其特征在于，所述步骤A包括以下步骤：

A1.根据式Ⅰ对多个测序图像分别进行预处理，得多个预处理后的图像；

A2.根据式Ⅱ对多个预处理后的图像分别进行傅立叶变换，得多个与测序图像相对应的傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\&Sigma;}}\underset{y=0}{\overset{N-1}{\&Sigma;}}g(x,y)e^{-j2\mathrm{\&pi;}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

4.根据权利要求2所述的比较多个测序图像清晰度的方法，其特征在于，所述步骤B包括以下步骤：

B1.分别对多个与测序图像相对应的傅立叶变换图进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

B2.将多个与测序图像相对应的傅立叶同态变换图中灰度级小于等于a的值重置为0；

B3.利用半径为r的圆形结构元对经步骤B2所得的结果图像进行闭运算；

B4.利用半径为R的圆形结构元对经步骤B3所得的结果图像进行开运算；

B5.利用边缘检测方法对经步骤B4所得的结果图像分别进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线；

所述10＜a＜15，R＞r。

5.根据权利要求4所述的比较多个测序图像清晰度的方法，其特征在于，所述步骤B1为：根据式Ⅳ对多个与测序图像相对应的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。

6.根据权利要求4所述的比较多个测序图像清晰度的方法，其特征在于，所述R等于1.5r、2r、2.5r或3r。

7.一种比较多个测序图像清晰度的系统，其特征在于，包括：谐波分析处理单元、闭合曲线拟合单元、面积比值计算单元和清晰度判定单元；

所述谐波分析处理单元用于对测序图像进行谐波分析处理得谐波分析处理图；

所述闭合曲线拟合单元用于根据所述谐波分析处理图的纹理特征拟合出闭合曲线；

所述面积比值计算单元用于确定所述闭合曲线对应的图形和相应的测序图像之间的面积比值；

所述清晰度判定单元用于根据面积比值的大小确定所述测序图像的清晰度；

所述测序图像为高通量基因测序图像或高通量基因测序图像中的一部分。

8.根据权利要求7所述的比较多个测序图像清晰度的系统，其特征在于，所述谐波分析处理单元用于对测序图像进行傅立叶变换得傅立叶变换图；所述闭合曲线拟合单元用于根据所述傅立叶变换图的纹理特征拟合出椭圆曲线或圆形曲线；所述面积比值计算单元用于确定所述椭圆曲线或圆形曲线所对应的椭圆或圆和相应的测序图像之间的面积比值。

9.根据权利要求8所述的比较多个测序图像清晰度的系统，其特征在于，还包括预处理单元；所述预处理单元用于根据式Ⅰ对测序图像进行预处理，得预处理后的图像；所述谐波分析处理单元用于根据式Ⅱ对预处理后的图像进行傅立叶变换，得傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\&Sigma;}}\underset{y=0}{\overset{N-1}{\&Sigma;}}g(x,y)e^{-j2\mathrm{\&pi;}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

10.根据权利要求8所述的比较多个测序图像清晰度的系统，其特征在于，所述闭合曲线拟合单元包括灰度级压缩模块、重置模块、闭运算模块、开运算模块和边缘提取曲线拟合模块；

所述灰度级压缩模块用于对傅立叶变换图进行灰度级压缩，得傅立叶同态变换图；

所述重置模块用于将傅立叶同态变换图中灰度级小于等于a的值重置为0；

所述闭运算模块用于利用半径为r的圆形结构元对经重置模块处理后所得的结果图像进行闭运算；

所述开运算模块用于利用半径为R的圆形结构元对经闭运算模块处理后所得的结果图像进行开运算；

所述边缘提取曲线拟合模块用于利用边缘检测方法对经开运算模块处理后所得的结果图像进行边缘提取，进而进行曲线拟合，得与测序图像相对应的椭圆曲线或圆形曲线；所述10＜a＜15，R＞r。

11.根据权利要求10所述的比较多个测序图像清晰度的系统，其特征在于，所述灰度级压缩模块用于根据式Ⅳ对傅立叶变换图进行灰度级压缩，得傅立叶同态变换图；G'(u,v)＝log(|G(u,v)|)Ⅳ。

12.一种自动聚焦控制方法，其特征在于，包括以下步骤：

A.调焦模块多次调节成像装置与采图位置之间的距离，并在每个距离均发出采集图像的指令给采图模块；

B.采图模块在接收到采集图像的指令后，利用成像装置采集测序图像，进而得到针对同一采图位置的多个测序图像；

C.图像处理模块在接收到经步骤B得到的多个测序图像后，根据权利要求1所述的方法，确定清晰度最高的测序图像，进而获得最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块；

D.调焦模块根据接收到的指令，控制成像装置与采图位置之间的距离为最佳聚焦距离。

13.根据权利要求12所述的自动聚焦控制方法，其特征在于，所述图像处理模块包括谐波分析处理单元、闭合曲线拟合单元、面积比值计算单元和最佳聚焦距离判定单元；

所述步骤C包括以下步骤：

C1.谐波分析处理单元对经步骤B所得的多个测序图像分别进行谐波分析处理，得多个与测序图像相对应的谐波分析处理图；

C2.闭合曲线拟合单元分别根据所述谐波分析处理图的纹理特征拟合出多个与测序图像相对应的闭合曲线；

C3.面积比值计算单元确定所述闭合曲线所对应的图形和其相应的测序图像之间的面积比值；

C4.所述最佳聚焦距离判定单元根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像，进而确定最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块。

14.根据权利要求13所述的自动聚焦控制方法，其特征在于，所述步骤C包括以下步骤：

C1.谐波分析处理单元对经步骤B所得的多个测序图像分别进行傅立叶变换，得多个与测序图像相对应的傅立叶变换图；

C2.闭合曲线拟合单元分别根据所述傅立叶变换图的纹理特征拟合出多个与测序图像相对应的椭圆曲线或圆形曲线；

C3.面积比值确定单元确定所述椭圆曲线或圆形曲线所对应的椭圆或圆和其相应的测序图像之间的面积比值；

C4.所述最佳聚焦距离判定单元根据面积比值的大小确定所述多个测序图像中清晰度最高的测序图像，进而确定最佳聚焦距离，并将调焦至最佳聚焦距离的指令发给调焦模块。

15.根据权利要求14所述的自动聚焦控制方法，其特征在于，所述图像处理模块还包括预处理单元；所述步骤C1包括以下步骤：

C11.所述预处理单元根据式Ⅰ对经步骤B所得的多个测序图像分别进行预处理，得多个预处理后的图像；

C12.所述谐波分析处理单元根据式Ⅱ对多个预处理后的图像分别进行傅立叶变换，得多个与测序图像相对应的傅立叶变换图；

g(x,y)＝f(x,y)(-1)^x+yⅠ；

$G(u,v)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\&Sigma;}}\underset{y=0}{\overset{N-1}{\&Sigma;}}g(x,y)e^{-j2\mathrm{\&pi;}(\frac{ux}{M}+\frac{vy}{N})}---II;$

其中，所述f(x,y)为测序图像的灰度函数；x，y代表了图像的坐标位置；所述G(u,v)为频域函数，M、N分别代表测序图像的宽和高。

16.根据权利要求14所述的自动聚焦控制方法，其特征在于，所述步骤C2包括以下步骤：

C21.对步骤C1所得的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

C22.将经步骤C21所得的傅立叶同态变换图中灰度级小于等于a的值重置为0；

C23.利用半径为r的圆形结构元对经步骤C22所得的结果图像进行闭运算；

C24.利用半径为R的圆形结构元对经步骤C23所得的结果图像进行开运算；

C25.利用边缘检测方法对经步骤C24所得的结果图像进行边缘提取，进而进行曲线拟合，得多个与测序图像相对应的椭圆曲线或圆形曲线；

所述10＜a＜15，R＞r。

17.根据权利要求16所述的自动聚焦的控制方法，其特征在于，所述步骤C21为：根据式Ⅳ对步骤C1所得的傅立叶变换图分别进行灰度级压缩，得多个与测序图像相对应的傅立叶同态变换图；

G'(u,v)＝log(|G(u,v)|)Ⅳ。