CN105243677A - 一种保证精度的实时单分子定位方法及系统 - Google Patents

Abstract

本发明公开了一种保证精度的实时单分子定位方法及系统。其中，该方法包括：先基于像素级定位参数对待处理图像进行区域提取；再对提取到的图像进行粗定位，得到提取到的图像的亚像素级定位参数；接着基于提取到的图像的亚像素级定位参数对提取到的图像进行精定位。由于经过粗定位得到的初始值(亚像素级定位参数)更接近于真实值，再据此进行精定位计算，因而减少了在精定位时的迭代次数，从而加快了超分辨图像的定位速度。这样，本发明在保证定位精度的情况下，加快了超分辨图像的定位速度。

Description

一种保证精度的实时单分子定位方法及系统

技术领域

本发明涉及超分辨定位成像技术领域，尤其涉及一种保证精度的实时单分子定位方法及系统。

背景技术

在荧光显微镜领域，超分辨显微成像技术由于突破了衍射极限，达到了数十纳米的分辨率，因而为生物学家提供了一种研究微小结构的有效方法，并被广泛应用。光激活定位显微技术(PhotoActivatedLocalizationMicroscopy,PALM)和随机光学重建显微法(StochasticOpticalReconstructionMicroscopy,STORM)是超分辨显微成像技术中的一类，它们基于全内反射的光学成像系统，加之软件算法，实现超分辨成像的纳米分辨率，并被统称为超分辨定位成像技术。

超分辨定位成像技术利用时间换取空间分辨率，一般由上千乃至上万帧的原始图像得到一幅超分辨图像。这样，图像处理速度就成了限制该技术发展的重要因素。正因如此，众多科学家们开始寻找在不牺牲空间分辨率的前提下能够提高图像处理速度的方法。

同时，科研级互补金属氧化物半导体传感器(Scientific-gradeComplementaryMetalOxideSemiconductor,sCMOS)由于其并行读出的优势，弥补了常用的探测器——电子倍增耦合器件(ElectronMultiplyingChargeCoupledDevice,EMCCD)串行读出速度慢的缺点，并逐渐应用于超分辨定位成像技术中。但是，这同时也带来了另一项挑战——由于全速率sCMOS的探测速度约是EMCCD的探测速度的47倍，因而sCMOS一秒钟能够采集到的图像数据量更多，这便为提高超分辨图像的定位速度增加了难度。

发明内容

本发明实施例通过提供一种保证精度的实时单分子定位方法及系统，在保证定位精度的情况下，加快了超分辨图像的定位速度的技术效果。

本发明实施例提供了一种保证精度的实时单分子定位方法，包括：

基于像素级定位参数对待处理图像进行区域提取；

对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数；

基于所述提取到的图像的亚像素级定位参数对所述提取到的图像进行精定位。

进一步地，所述基于像素级定位参数对待处理图像进行区域提取，包括：

获取所述待处理图像中每个像素点的像素值；

判断所述获取到的像素值是否大于预设阈值；

若是，获取最大像素值所对应的像素点的坐标，所述最大像素值所对应的像素点的坐标为所述像素级定位参数；

提取出以所述获取到的最大像素值所对应的像素点的坐标为中心的一定区域的图像；

其中，所述一定区域由所述待处理图像而定。

进一步地，所述对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数，包括：

通过公式 $C_x={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(x_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 和 $C_y={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(y_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 得到所述提取到的图像中像素点的X坐标C_x和Y坐标C_y；其中，I_i,j表示像素点(i，j)处的灰度值，x_i表示像素点(i，j)在x轴方向的坐标，y_i表示像素点(i，j)在y轴方向的坐标，n表示所述提取到的图像的像素的个数；

获取所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；

其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

进一步地，所述对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数，包括：

将所述提取到的图像与X轴的梯度算子 $\left[\begin{array}{ccccc}-1& 1& 0& 1& 1\\ -1& -1& 0& 1& 1\\ -1& -1& 0& 1& 1\end{array}\right]$ 进行卷积运算，得到x轴的梯度分布Gx；

将所述提取到的图像与Y轴的梯度算子 $\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\\ -1& -1& -1\end{array}\right]$ 进行卷积运算，得到y轴的梯度分布Gy；

通过公式Θ(i,j)＝arctan2(Gy(i,j)/Gx(i,j))，计算得到像素梯度方向Θ(i,j)；

通过公式Cy＝K(i,j)*Cx+C(i,j)计算得到平行于所述像素梯度方向Θ(i,j)的直线，即辐射线；

其中，K(i,j)为所述辐射线的斜率，由所述像素梯度方向Θ(i,j)决定；C(i,j)为所述辐射线在y轴上的截距；

通过求解方程 $D_{sum}=\mathrm{\Σ}\mathrm{\Σ}{\left(\frac{|K(i,j)*Cx-Cy+C(i,j)|}{\sqrt{1+K{(i,j)}^2}}\right)}^2$ 和 $\left\{\begin{array}{c}\frac{\∂D_{sum}}{\∂Cx}=0\\ \frac{\∂D_{sum}}{\∂Cy}=0\end{array}\right.$ 得到与所述辐射线之间的距离的平方和最小的像素点的X坐标C_x和Y坐标C_y；

获取所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；

其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

进一步地，所述基于所述提取到的图像的亚像素级定位参数对所述提取到的图像进行精定位，包括：

建立计算模型： ${\stackrel{\‾}{I}}_{i,j}=A_0\exp \[-\frac{{(i-C_x)}^2+{(j-C_y)}^2}{{s_0}^2}\]+b_0;$

其中，为像素点(i,j)处的理论光强，A为所述信号亮度的估计值A₀，i为所述提取到的区域的x轴坐标，j为所述提取到的区域的y轴坐标，x₀为所述提取到的图像中像素点的X坐标C_x，y₀为所述提取到的图像中像素点的Y坐标C_y，s为所述荧光点扩散函数的估计标准差s₀，b为所述背景估计值b₀；

对所述计算模型进行估计，即以所述信号亮度的估计值A₀、所述提取到的图像中像素点的X坐标C_x、所述提取到的图像中像素点的Y坐标C_y、所述荧光点扩散函数的估计标准差s₀和所述背景估计值b₀作为初值代入所述计算模型进行迭代计算直至与观测值的差值在误差阈值之内，将此时的所述信号亮度的估计值A₀记作信号亮度的精确值A，所述提取到的图像中像素点的X坐标C_x记作像素点精度的X坐标x₀、所述提取到的图像中像素点的Y坐标C_y记作像素点精度的Y坐标y₀、所述荧光点扩散函数的估计标准差s₀记作荧光点扩散函数精确的标准差s，所述背景估计值b₀记作背景精确值b。

本发明实施例提供的保证精度的实时单分子定位系统，包括：

区域提取模块，用于基于像素级定位参数对待处理图像进行区域提取；

粗定位模块，用于对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数；

精定位模块，用于基于所述提取到的图像的亚像素级定位参数对所述提取到的图像进行精定位。

进一步地，所述区域提取模块，具体包括：

像素值获取单元，用于获取所述待处理图像中每个像素点的像素值；

判断单元，用于判断所述获取到的像素值是否大于预设阈值；

像素点获取单元，用于若所述判断单元的判断结果为是，获取最大像素值所对应的像素点的坐标，所述最大像素值所对应的像素点的坐标为所述像素级定位参数；

区域提取执行单元，用于提取出以所述获取到的最大像素值所对应的像素点的坐标为中心的一定区域的图像；其中，所述一定区域由所述待处理图像而定。

进一步地，所述粗定位模块，具体包括：

第一计算单元，用于通过公式和得到所述提取到的图像中像素点的X坐标C_x和Y坐标C_y；其中，I_i,j表示像素点(i，j)处的灰度值，x_i表示像素点(i，j)在x轴方向的坐标，y_i表示像素点(i，j)在y轴方向的坐标，n表示所述提取到的图像的像素的个数；

第二计算单元，用于计算所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

第三计算单元，用于获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

第四计算单元，用于通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

进一步地，所述粗定位模块，具体包括：

第一运算单元，用于将所述提取到的图像与X轴的梯度算子 $\left[\begin{array}{ccccc}-1& 1& 0& 1& 1\\ -1& -1& 0& 1& 1\\ -1& -1& 0& 1& 1\end{array}\right]$ 进行卷积运算，得到x轴的梯度分布Gx；

第二运算单元，用于将所述提取到的图像与Y轴的梯度算子 $\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\\ -1& -1& -1\end{array}\right]$ 进行卷积运算，得到y轴的梯度分布Gy；

第三运算单元，用于通过公式Θ(i,j)＝arctan2(Gy(i,j)/Gx(i,j))，计算得到像素梯度方向Θ(i,j)；

第四运算单元，用于通过公式Cy＝K(i,j)*Cx+C(i,j)计算得到平行于所述像素梯度方向Θ(i,j)的直线，即辐射线；其中，K(i,j)为所述辐射线的斜率，由所述像素梯度方向Θ(i,j)决定；C(i,j)为所述辐射线在y轴上的截距；

第五运算单元，用于通过求解方程 $D_{sum}=\mathrm{\Σ}\mathrm{\Σ}{\left(\frac{|K(i,j)*Cx-Cy+C(i,j)|}{\sqrt{1+K{(i,j)}^2}}\right)}^2$ 和 $\left\{\begin{array}{c}\frac{\∂D_{sum}}{\∂Cx}=0\\ \frac{\∂D_{sum}}{\∂Cy}=0\end{array}\right.$ 得到与所述辐射线之间的距离的平方和最小的像素点的X坐标C_x和Y坐标C_y；

第六运算单元，用于获取所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

第七运算单元，用于获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

第八运算单元，用于通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

进一步地，所述精定位模块，具体包括：

计算模型建立单元，用于建立计算模型： ${\stackrel{\‾}{I}}_{i,j}=A_0\exp \[-\frac{{(i-C_x)}^2+{(j-C_y)}^2}{{s_0}^2}\]+b_0;$ 其中，为像素点(i,j)处的理论光强，A为所述信号亮度的估计值A₀，i为所述提取到的区域的x轴坐标，j为所述提取到的区域的y轴坐标，x₀为所述提取到的图像中像素点的X坐标C_x，y₀为所述提取到的图像中像素点的Y坐标C_y，s为所述荧光点扩散函数的估计标准差s₀，b为所述背景估计值b₀；

精定位执行单元，用于对所述计算模型进行估计，即以所述信号亮度的估计值A₀、所述提取到的图像中像素点的X坐标C_x、所述提取到的图像中像素点的Y坐标C_y、所述荧光点扩散函数的估计标准差s₀和所述背景估计值b₀作为初值代入所述计算模型进行迭代计算直至与观测值的差值在误差阈值之内，将此时的所述信号亮度的估计值A₀记作信号亮度的精确值A，所述提取到的图像中像素点的X坐标C_x记作像素点精度的X坐标x₀、所述提取到的图像中像素点的Y坐标C_y记作像素点精度的Y坐标y₀、所述荧光点扩散函数的估计标准差s₀记作荧光点扩散函数精确的标准差s，所述背景估计值b₀记作背景精确值b。

本发明实施例中提供的一个或多个技术方案，至少具有如下技术效果或优点：

1、先基于像素级定位参数对待处理图像进行区域提取；再对提取到的图像进行粗定位，得到提取到的图像的亚像素级定位参数；接着基于提取到的图像的亚像素级定位参数对提取到的图像进行精定位。由于经过粗定位得到的初始值(亚像素级定位参数)更接近于真实值，再据此进行精定位计算，因而减少了在精定位时的迭代次数，从而加快了超分辨图像的定位速度。这样，本发明实施例在保证定位精度的情况下，加快了超分辨图像的定位速度。

2、若将区域提取、粗定位和精定位的步骤在具有并行处理能力和流处理能力的硬件(如FPGA(ProgrammableGateArray，现场可编程门阵列)或ASIC(ApplicationSpecificIntegratedCircuit，特定用途集成电路)等)中进行，由于在具有并行处理能力的硬件中处理数据，待处理的数据之间的操作属于并行操作，因而进一步加快了超分辨图像的定位速度。同时，每个硬件时钟内，数据都会向前传输，这样，来一个数据就可以处理一个数据，数据处理就处于流水线处理(流处理)的状态，这样，就可以做超分辨图像的实时处理，更进一步地加快了超分辨图像的定位速度。

3、若将粗定位模块和精定位模块配置在具有并行处理能力和流处理能力的硬件上，粗定位模块和精定位模块就可以配置在不同的硬件单元上，使粗定位模块和精定位模块能够相互独立地工作。这样，粗定位模块就不占用精定位模块的资源，从而缩短了超分辨图像定位的时间。

附图说明

图1为本发明实施例一提供的保证精度的实时单分子定位方法的流程图；

图2为本发明实施例二提供的保证精度的实时单分子定位系统的模块图；

图3为通过本发明实施例二提供的系统对超分辨图像进行定位的流程图。

具体实施方式

本发明实施例通过提供一种保证精度的实时单分子定位方法及系统，在保证定位精度的情况下，加快了超分辨图像的定位速度的技术效果。

本发明实施例中的技术方案为解决上述技术问题，总体思路如下：

先基于像素级定位参数对待处理图像进行区域提取；再对提取到的图像进行粗定位，得到提取到的图像的亚像素级定位参数；接着基于提取到的图像的亚像素级定位参数对提取到的图像进行精定位。由于经过粗定位得到的初始值(亚像素级定位参数)更接近于真实值，再据此进行精定位计算，因而减少了在精定位时的迭代次数，从而加快了超分辨图像的定位速度。这样，本发明实施例在保证定位精度的情况下，加快了超分辨图像的定位速度。

为了更好地理解上述技术方案，下面将结合说明书附图以及具体的实施方式对上述技术方案进行详细的说明。

实施例一

参见图1，本发明实施例提供的保证精度的实时单分子定位方法，包括：

步骤S110：基于像素级定位参数对待处理图像进行区域提取；

对本步骤进行说明，步骤S110具体包括：

获取待处理图像中每个像素点的像素值；

判断获取到的像素值是否大于预设阈值；

若是，获取其中的最大像素值B所对应的像素点的坐标，最大像素值B所对应的像素点的坐标为像素级定位参数；

提取出以获取到的最大像素值B所对应的像素点的坐标为中心的一定区域的图像；

其中，一定区域由待处理图像而定。

为了提高本发明实施例的超分辨图像的定位精度，在步骤S110之前，还包括：

对待处理图像进行去噪处理；

在本实施例中，对待处理图像进行去噪处理的步骤具体包括：

将待处理图像与均值模板进行卷积运算，得到去噪后的待处理图像，即对待处理图像进行去噪处理。

步骤S120：对提取到的图像进行粗定位，得到提取到的图像的亚像素级定位参数；

本发明实施例至少提出了两种对提取到的图像进行粗定位的方法，具体地，第一种粗定位的方法包括：

通过公式 $C_x={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(x_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 和 $C_y={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(y_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 得到提取到的图像中像素点的X坐标C_x和Y坐标C_y；其中，I_i,j表示像素点(i，j)处的灰度值，x_i表示像素点(i，j)在x轴方向的坐标，y_i表示像素点(i，j)在y轴方向的坐标，n表示提取到的图像的像素的个数。

获取提取到的图像的外围像素点的像素值的平均值，将平均值作为背景估计值b₀；

获取最大像素值B，将最大像素值B减去背景估计值b₀，得到信号亮度的估计值A₀；

通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；

其中，λ为全内反射的光学成像系统的光的波长，由全内反射的光学成像系统中使用的光源与滤光片决定；NA为数值孔径，由全内反射的光学成像系统中使用的物镜决定；pixelsize为像素尺寸的大小，由全内反射的光学成像系统的放大倍数和探测器的像素大小决定。

第二种粗定位的方法包括：

将提取到的图像与X轴的梯度算子 $\left[\begin{array}{ccccc}-1& 1& 0& 1& 1\\ -1& -1& 0& 1& 1\\ -1& -1& 0& 1& 1\end{array}\right]$ 进行卷积运算，得到x轴的梯度分布Gx；

将提取到的图像与Y轴的梯度算子 $\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\\ -1& -1& -1\end{array}\right]$ 进行卷积运算，得到y轴的梯度分布Gy；

通过公式Θ(i,j)＝arctan2(Gy(i,j)/Gx(i,j))，计算得到像素梯度方向Θ(i,j)；

通过公式Cy＝K(i,j)*Cx+C(i,j)计算得到平行于像素梯度方向Θ(i,j)的直线，即辐射线；

其中，K(i,j)为辐射线的斜率，由像素梯度方向Θ(i,j)决定，即K(i，j)＝tanΘ(i,j)；C(i,j)为辐射线在y轴上的截距；

通过求解方程 $D_{sum}=\mathrm{\Σ}\mathrm{\Σ}{\left(\frac{|K(i,j)*Cx-Cy+C(i,j)|}{\sqrt{1+K{(i,j)}^2}}\right)}^2$ 和 $\left\{\begin{array}{c}\frac{\∂D_{sum}}{\∂Cx}=0\\ \frac{\∂D_{sum}}{\∂Cy}=0\end{array}\right.$ 得到与辐射线之间的距离的平方和最小的像素点的X坐标C_x和Y坐标C_y；

在本实施例中，求解得到的该像素点为最大辐射对称中心点。

获取提取到的图像的外围像素点的像素值的平均值，将平均值作为背景估计值b₀；

获取最大像素值B，将最大像素值B减去背景估计值b₀，得到信号亮度的估计值A₀；

通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；

其中，λ为全内反射的光学成像系统的光的波长，由全内反射的光学成像系统中使用的光源与滤光片决定；NA为数值孔径，由全内反射的光学成像系统中使用的物镜决定；pixelsize为像素尺寸的大小，由全内反射的光学成像系统的放大倍数和探测器的像素大小决定。

步骤S130：基于提取到的图像的亚像素级定位参数对提取到的图像进行精定位。

对本步骤进行说明，步骤S130具体包括：

建立计算模型： ${\stackrel{\‾}{I}}_{i,j}=A_0\exp \[-\frac{{(i-C_x)}^2+{(j-C_y)}^2}{{s_0}^2}\]+b_0;$

其中，为像素点(i,j)处的理论光强，A为信号亮度的估计值A₀，i为提取到的区域的x轴坐标，j为提取到的区域的y轴坐标，x₀为提取到的图像中像素点的X坐标C_x，y₀为提取到的图像中像素点的Y坐标C_y，s为荧光点扩散函数的估计标准差s₀，b为背景估计值b₀；

对计算模型进行估计，即以信号亮度的估计值A₀、提取到的图像中像素点的X坐标C_x、提取到的图像中像素点的Y坐标C_y、荧光点扩散函数的估计标准差s₀和背景估计值b₀作为初值代入计算模型进行迭代计算直至与观测值的差值在误差阈值之内，将此时的信号亮度的估计值A₀记作信号亮度的精确值A，提取到的图像中像素点的X坐标C_x记作像素点精度的X坐标x₀、提取到的图像中像素点的Y坐标C_y记作像素点精度的Y坐标y₀、荧光点扩散函数的估计标准差s₀记作荧光点扩散函数精确的标准差s，背景估计值b₀记作背景精确值b。

具体地，在本实施例中，采用极大似然法对上述计算模型进行估计，即调整参数A,x_o,y_o,s,b使理论光强(即理论灰度值)逼近观测值I(x_i,j)(即实际灰度值)，得到的估计方程如下：

${(A,x_o,y_o,s,b)}^{*}=\arg \max (-\underset{i,j}{\Σ}{\stackrel{\‾}{I}}_{i,j}+\underset{i,j}{\Σ}I\left(x_{i,j}\right)\log \left({\stackrel{\‾}{I}}_{i,j}\right)).$

为了得到精定位之后的图像，在步骤S130之后，还包括：

基于得到的精确的像素点的X坐标x₀和Y坐标y₀、背景估计值A、信号亮度的估计值b和荧光点PSF标准差估计值s₀，进行图像渲染，得到最终的图像。

这里需要说明的是，对待处理图像进行去噪处理的步骤、步骤S110和步骤S120均在FPGA中完成，步骤S130在电脑的GPU(GraphicsProcessingUnit，图形处理器)中完成，得到最终的图像的步骤在电脑的CPU(CentralProcessingUnit，中央处理器)中完成。

实施例二

参见图2，本发明实施例提供的保证精度的实时单分子定位系统，包括：

区域提取模块100，用于基于像素级定位参数对待处理图像进行区域提取；

在本实施例中，区域提取模块100，具体包括：

像素值获取单元，用于获取待处理图像中每个像素点的像素值；

判断单元，用于判断获取到的像素值是否大于预设阈值；

像素点获取单元，用于若判断单元的判断结果为是，获取其中的最大像素值B所对应的像素点的坐标，最大像素值B所对应的像素点的坐标为像素级定位参数；

区域提取执行单元，用于提取出以获取到的最大像素值B所对应的像素点的坐标为中心的一定区域的图像；其中，一定区域由待处理图像而定。

为了提高本发明实施例的超分辨图像的定位精度，还包括：

去噪模块，用于对待处理图像进行去噪处理；

在本实施例中，去噪模块，具体用于将待处理图像与均值模板进行卷积运算，得到去噪后的待处理图像，即对待处理图像进行去噪处理。

粗定位模块200，用于对提取到的图像进行粗定位，得到提取到的图像的亚像素级定位参数；

本发明实施例至少提出了两种粗定位模块200，具体地，第一种粗定位模块200，具体包括：

第一计算单元，用于通过公式和得到提取到的图像中像素点的X坐标C_x和Y坐标C_y；其中，I_i,j表示像素点(i，j)处的灰度值，x_i表示像素点(i，j)在x轴方向的坐标，y_i表示像素点(i，j)在y轴方向的坐标，n表示提取到的图像的像素的个数。

第二计算单元，用于计算提取到的图像的外围像素点的像素值的平均值，将平均值作为背景估计值b₀；

第三计算单元，用于获取最大像素值B，将最大像素值B减去背景估计值b₀，得到信号亮度的估计值A₀；

第四计算单元，用于通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

在本实施例中，λ由全内反射的光学成像系统中使用的光源与滤光片决定，NA由全内反射的光学成像系统中使用的物镜决定，pixelsize由全内反射的光学成像系统的放大倍数和探测器的像素大小决定。

第二种粗定位模块200，具体包括：

第一运算单元，用于将提取到的图像与X轴的梯度算子 $\left[\begin{array}{ccccc}-1& 1& 0& 1& 1\\ -1& -1& 0& 1& 1\\ -1& -1& 0& 1& 1\end{array}\right]$ 进行卷积运算，得到x轴的梯度分布Gx；

第二运算单元，用于将提取到的图像与Y轴的梯度算子 $\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\\ -1& -1& -1\end{array}\right]$ 进行卷积运算，得到y轴的梯度分布Gy；

第三运算单元，用于通过公式Θ(i,j)＝arctan2(Gy(i,j)/Gx(i,j))，计算得到像素梯度方向Θ(i,j)；

第四运算单元，用于通过公式Cy＝K(i,j)*Cx+C(i,j)计算得到平行于像素梯度方向Θ(i,j)的直线，即辐射线；其中，K(i,j)为辐射线的斜率，由像素梯度方向Θ(i,j)决定，即K(i，j)＝tanΘ(i,j)；C(i,j)为辐射线在y轴上的截距；

第五运算单元，用于通过求解方程 $D_{sum}=\mathrm{\Σ}\mathrm{\Σ}{\left(\frac{|K(i,j)*Cx-Cy+C(i,j)|}{\sqrt{1+K{(i,j)}^2}}\right)}^2$ 和 $\left\{\begin{array}{c}\frac{\∂D_{sum}}{\∂Cx}=0\\ \frac{\∂D_{sum}}{\∂Cy}=0\end{array}\right.$ 得到与辐射线之间的距离的平方和最小的像素点的X坐标C_x和Y坐标C_y；

第六运算单元，用于计算提取到的图像的外围像素点的像素值的平均值，将平均值作为背景估计值b₀；

第七运算单元，用于获取最大像素值B，将最大像素值B减去背景估计值b₀，得到信号亮度的估计值A₀；

第八运算单元，用于通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

在本实施例中，λ由全内反射的光学成像系统中使用的光源与滤光片决定，NA由全内反射的光学成像系统中使用的物镜决定，pixelsize由全内反射的光学成像系统的放大倍数和探测器的像素大小决定。

精定位模块300，用于基于提取到的图像的亚像素级定位参数对提取到的图像进行精定位。

在本实施例中，精定位模块300，具体包括：

计算模型建立单元，用于建立计算模型： ${\stackrel{\‾}{I}}_{i,j}=A_0\exp \[-\frac{{(i-C_x)}^2+{(j-C_y)}^2}{{s_0}^2}\]+b_0;$ 其中，为像素点(i,j)处的理论光强，A为信号亮度的估计值A₀，i为提取到的区域的x轴坐标，j为提取到的区域的y轴坐标，x₀为提取到的图像中像素点的X坐标C_x，y₀为提取到的图像中像素点的Y坐标C_y，s为荧光点扩散函数的估计标准差s₀，b为背景估计值b₀；

精定位执行单元，用于对计算模型进行估计，即以信号亮度的估计值A₀、提取到的图像中像素点的X坐标C_x、提取到的图像中像素点的Y坐标C_y、荧光点扩散函数的估计标准差s₀和背景估计值b₀作为初值代入计算模型进行迭代计算直至与观测值的差值在误差阈值之内，将此时的信号亮度的估计值A₀记作信号亮度的精确值A，提取到的图像中像素点的X坐标C_x记作像素点精度的X坐标x₀、提取到的图像中像素点的Y坐标C_y记作像素点精度的Y坐标y₀、荧光点扩散函数的估计标准差s₀记作荧光点扩散函数精确的标准差s，背景估计值b₀记作背景精确值b。

具体地，在本实施例中，采用极大似然法对上述计算模型进行估计，即调整参数A,x_o,y_o,s,b使理论光强(即理论灰度值)逼近观测值I(x_i,j)(即实际灰度值)，得到的估计方程如下：

${(A,x_o,y_o,s,b)}^{*}=\arg \max (-\underset{i,j}{\Σ}{\stackrel{\‾}{I}}_{i,j}+\underset{i,j}{\Σ}I\left(x_{i,j}\right)\log \left({\stackrel{\‾}{I}}_{i,j}\right)).$

为了得到精定位之后的图像，还包括：

显示模块，用于基于得到的精确的像素点的X坐标x₀和Y坐标y₀、背景估计值A、信号亮度的估计值b和荧光点PSF标准差估计值s₀，进行图像渲染，得到最终的图像。

这里需要说明的是，去噪模块、区域提取模块100和粗定位模块200均在FPGA中，精定位模块300在电脑的GPU中，显示模块在电脑的CPU中。

实施例三

通过上述实施例二提供的系统对超分辨图像进行定位的具体步骤包括：

参见图3，先通过探测器对原始图像数据进行探测，并将获得的原始图像数据通过一分二的CameraLink子卡，一路传输到电脑端，进行传统的控制与数据传输，另一路将图像数据传输到FPGA板卡上，并在FPGA板卡上的DRAM缓存区中进行缓存。再将DRAM缓存区中的图像数据读出到行缓存区中，得到3x3的待处理图像像素。去噪模块使用3x3的均值模板与得到的3x3待处理图像像素进行卷积，得到滤波后的图像。去噪模块将滤波后的图像传入到区域提取模块100，进行荧光点的判断。其判断方法为先求得一列滤波后图像的灰度值的标准差，并以该标准差的5倍作为预设阈值。再判断滤波后的图像的像素点的像素值是否大于该预设阈值，同时判断其是否为其3x3像素临域中的最大值。若是，则将由原始图像的行缓存区扫描到的像素以该像素坐标为中心的7x7区域提取出来，得到区域提取后的图像数据；若不是，则不对该区域做任何处理。接着，求得区域提取的图像数据的最外围一周的24个像素的平均值作为背景估计值b₀，以中心像素值减去背景估计值b₀作为信号亮度的估计值A₀，荧光点扩散函数(PSF)标准差估计值使用固定的常数s₀进行估计，通过公式 $C_x={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(x_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 和 $C_y={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(y_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 计算出提取到的图像中像素点的X坐标C_x和Y坐标C_y，完成粗定位。再将粗定位得到的背景估计值b₀、信号亮度的估计值A₀、荧光点PSF标准差估计值s₀、XY坐标估计值C_x、C_y和区域提取得到的7x7区域通过通用接口(PCIE或者以太网，USB)传输到电脑端，并通过电脑内存传输给GPU进行精定位。最后，将精定位得到的精确的像素点的X坐标x₀和Y坐标y₀、背景估计值A、信号亮度的估计值b和荧光点扩散函数PSF传入到CPU中，调用系统渲染函数进行图像渲染，得到最终的图像。

【技术效果】

1、先基于像素级别的分子定位坐标对待处理图像进行区域提取；再对提取到的图像进行粗定位，得到提取到的图像的亚像素级别的分子定位坐标C_x和C_y、信号亮度的估计值A₀、荧光点扩散函数的估计标准差s₀和背景估计值b₀；接着基于提取到的图像的亚像素级别的分子定位坐标C_x和C_y、信号亮度的估计值A₀、荧光点扩散函数的估计标准差s₀和背景估计值b₀对提取到的图像进行精定位，得到精确定位的像素点的X坐标x₀和Y坐标y₀、背景估计值A、信号亮度的估计值b和荧光点PSF标准差估计值s₀。由于经过粗定位得到的初始值(亚像素级定位参数)更接近于真实值，再据此进行精定位计算，因而减少了在精定位时的迭代次数，从而加快了超分辨图像的定位速度。这样，本发明实施例在保证定位精度的情况下，加快了超分辨图像的定位速度。

2、在对待处理图像进行区域提取之前，对待处理图像进行去噪处理，进一步地保证了本发明实施例的定位精度。

3、在本发明实施例中，去噪、区域提取、粗定位的步骤在FPGA中完成，精定位和显示的步骤在电脑上的GPU中完成，由于FPGA和GPU均具有数据并行处理的能力，而在FPGA中可以将同一功能配置多个，将数据拆分成多份，分别进行处理，因而可进一步加快超分辨图像的定位速度。同时，每个硬件时钟内，数据都会向前传输，这样，来一个数据就可以处理一个数据，数据处理就处于流水线处理(流处理)的状态，这样，就可以做超分辨图像的实时处理，更进一步地加快了超分辨图像的定位速度。

4、在本发明实施例中，若将粗定位模块和精定位模块配置在具有并行处理能力和流处理能力的硬件上，粗定位模块和精定位模块就可以配置在不同的硬件单元上，使粗定位模块和精定位模块能够相互独立地工作。这样，粗定位模块就不占用精定位模块的资源，从而缩短了超分辨图像定位的时间。

本发明实施例中的硬件的加速时钟可以达到150MHz，这相对于现有的探测器的速度来说，已经完全能满足超分辨图像实时定位的要求。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (10)

1.一种保证精度的实时单分子定位方法，其特征在于，包括：

基于像素级定位参数对待处理图像进行区域提取；

对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数；

基于所述提取到的图像的亚像素级定位参数对所述提取到的图像进行精定位。

2.如权利要求1所述的方法，其特征在于，所述基于像素级定位参数对待处理图像进行区域提取，包括：

获取所述待处理图像中每个像素点的像素值；

判断所述获取到的像素值是否大于预设阈值；

若是，获取最大像素值所对应的像素点的坐标，所述最大像素值所对应的像素点的坐标为所述像素级定位参数；

提取出以所述获取到的最大像素值所对应的像素点的坐标为中心的一定区域的图像；

其中，所述一定区域由所述待处理图像而定。

3.如权利要求2所述的方法，其特征在于，所述对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数，包括：

通过公式 $C_x={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(x_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 和 $C_y={\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n(y_i{I}_{i,j)}\backslash {\mathrm{\Σ}}_{i=1}^n{\mathrm{\Σ}}_{j=1}^n\left(I_{i,j}\right)$ 得到所述提取到的图像中像素点的X坐标C_x和Y坐标C_y；其中，I_i,j表示像素点(i，j)处的灰度值，x_i表示像素点(i，j)在x轴方向的坐标，y_i表示像素点(i，j)在y轴方向的坐标，n表示所述提取到的图像的像素的个数；

获取所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；

其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

4.如权利要求2所述的方法，其特征在于，所述对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数，包括：

将所述提取到的图像与X轴的梯度算子 $\left[\begin{array}{ccccc}-1& 1& 0& 1& 1\\ -1& -1& 0& 1& 1\\ -1& -1& 0& 1& 1\end{array}\right]$ 进行卷积运算，得到x轴的梯度分布Gx；

将所述提取到的图像与Y轴的梯度算子 $\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\\ -1& -1& -1\end{array}\right]$ 进行卷积运算，得到y轴的梯度分布Gy；

通过公式Θ(i,j)＝arctan2(Gy(i,j)/Gx(i,j))，计算得到像素梯度方向Θ(i,j)；

通过公式Cy＝K(i,j)*Cx+C(i,j)计算得到平行于所述像素梯度方向Θ(i,j)的直线，即辐射线；

其中，K(i,j)为所述辐射线的斜率，由所述像素梯度方向Θ(i,j)决定；C(i,j)为所述辐射线在y轴上的截距；

通过求解方程 $D_{sum}=\mathrm{\Σ}\mathrm{\Σ}{\left(\frac{|K(i,j)*Cx-Cy+C(i,j)|}{\sqrt{1+K{(i,j)}^2}}\right)}^2$ 和 $\left\{\begin{array}{c}\frac{\∂D_{sum}}{\∂Cx}=0\\ \frac{\∂D_{sum}}{\∂Cy}=0\end{array}\right.$ 得到与所述辐射线之间的距离的平方和最小的像素点的X坐标C_x和Y坐标C_y；

获取所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；

其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

5.如权利要求3或4所述的方法，其特征在于，所述基于所述提取到的图像的亚像素级定位参数对所述提取到的图像进行精定位，包括：

建立计算模型： ${\stackrel{\‾}{I}}_{i,j}=A_0\exp \[-\frac{{(i-C_x)}^2+{(j-C_y)}^2}{{s_0}^2}\]+b_0;$

其中，为像素点(i,j)处的理论光强，A为所述信号亮度的估计值A₀，i为所述提取到的区域的x轴坐标，j为所述提取到的区域的y轴坐标，x₀为所述提取到的图像中像素点的X坐标C_x，y₀为所述提取到的图像中像素点的Y坐标C_y，s为所述荧光点扩散函数的估计标准差s₀，b为所述背景估计值b₀；

对所述计算模型进行估计，即以所述信号亮度的估计值A₀、所述提取到的图像中像素点的X坐标C_x、所述提取到的图像中像素点的Y坐标C_y、所述荧光点扩散函数的估计标准差s₀和所述背景估计值b₀作为初值代入所述计算模型进行迭代计算直至与观测值的差值在误差阈值之内，将此时的所述信号亮度的估计值A₀记作信号亮度的精确值A，所述提取到的图像中像素点的X坐标C_x记作像素点精度的X坐标x₀、所述提取到的图像中像素点的Y坐标C_y记作像素点精度的Y坐标y₀、所述荧光点扩散函数的估计标准差s₀记作荧光点扩散函数精确的标准差s，所述背景估计值b₀记作背景精确值b。

6.一种保证精度的实时单分子定位系统，其特征在于，包括：

区域提取模块，用于基于像素级定位参数对待处理图像进行区域提取；

粗定位模块，用于对提取到的图像进行粗定位，得到所述提取到的图像的亚像素级定位参数；

精定位模块，用于基于所述提取到的图像的亚像素级定位参数对所述提取到的图像进行精定位。

7.如权利要求6所述的系统，其特征在于，所述区域提取模块，具体包括：

像素值获取单元，用于获取所述待处理图像中每个像素点的像素值；

判断单元，用于判断所述获取到的像素值是否大于预设阈值；

像素点获取单元，用于若所述判断单元的判断结果为是，获取最大像素值所对应的像素点的坐标，所述最大像素值所对应的像素点的坐标为所述像素级定位参数；

区域提取执行单元，用于提取出以所述获取到的最大像素值所对应的像素点的坐标为中心的一定区域的图像；其中，所述一定区域由所述待处理图像而定。

8.如权利要求7所述的系统，其特征在于，所述粗定位模块，具体包括：

第一计算单元，用于通过公式和得到所述提取到的图像中像素点的X坐标C_x和Y坐标C_y；其中，I_i,j表示像素点(i，j)处的灰度值，x_i表示像素点(i，j)在x轴方向的坐标，y_i表示像素点(i，j)在y轴方向的坐标，n表示所述提取到的图像的像素的个数；

第二计算单元，用于计算所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

第三计算单元，用于获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

第四计算单元，用于通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

9.如权利要求7所述的系统，其特征在于，所述粗定位模块，具体包括：

第一运算单元，用于将所述提取到的图像与X轴的梯度算子 $\left[\begin{array}{ccccc}-1& 1& 0& 1& 1\\ -1& -1& 0& 1& 1\\ -1& -1& 0& 1& 1\end{array}\right]$ 进行卷积运算，得到x轴的梯度分布Gx；

第二运算单元，用于将所述提取到的图像与Y轴的梯度算子 $\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\\ -1& -1& -1\end{array}\right]$ 进行卷积运算，得到y轴的梯度分布Gy；

第三运算单元，用于通过公式Θ(i,j)＝arctan2(Gy(i,j)/Gx(i,j))，计算得到像素梯度方向Θ(i,j)；

第四运算单元，用于通过公式Cy＝K(i,j)*Cx+C(i,j)计算得到平行于所述像素梯度方向Θ(i,j)的直线，即辐射线；其中，K(i,j)为所述辐射线的斜率，由所述像素梯度方向Θ(i,j)决定；C(i,j)为所述辐射线在y轴上的截距；

第五运算单元，用于通过求解方程 $D_{sum}=\mathrm{\Σ}\mathrm{\Σ}{\left(\frac{|K(i,j)*Cx-Cy+C(i,j)|}{\sqrt{1+K{(i,j)}^2}}\right)}^2$ 和 $\left\{\begin{array}{c}\frac{\∂D_{sum}}{\∂Cx}=0\\ \frac{\∂D_{sum}}{\∂Cy}=0\end{array}\right.$ 得到与所述辐射线之间的距离的平方和最小的像素点的X坐标C_x和Y坐标C_y；

第六运算单元，用于获取所述提取到的图像的外围像素点的像素值的平均值，将所述平均值作为背景估计值b₀；

第七运算单元，用于获取所述最大像素值，将所述最大像素值减去所述背景估计值b₀，得到信号亮度的估计值A₀；

第八运算单元，用于通过公式s₀＝0.21*λ/(NA*pixelsize)计算得到荧光点扩散函数的估计标准差s₀；其中，λ为全内反射的光学成像系统的光的波长，NA为数值孔径，pixelsize为像素尺寸的大小。

10.如权利要求8或9所述的系统，其特征在于，所述精定位模块，具体包括：

计算模型建立单元，用于建立计算模型： ${\stackrel{\‾}{I}}_{i,j}=A_0\exp \[-\frac{{(i-C_x)}^2+{(j-C_y)}^2}{{s_0}^2}\]+b_0;$ 其中，为像素点(i,j)处的理论光强，A为所述信号亮度的估计值A₀，i为所述提取到的区域的x轴坐标，j为所述提取到的区域的y轴坐标，x₀为所述提取到的图像中像素点的X坐标C_x，y₀为所述提取到的图像中像素点的Y坐标C_y，s为所述荧光点扩散函数的估计标准差s₀，b为所述背景估计值b₀；

精定位执行单元，用于对所述计算模型进行估计，即以所述信号亮度的估计值A₀、所述提取到的图像中像素点的X坐标C_x、所述提取到的图像中像素点的Y坐标C_y、所述荧光点扩散函数的估计标准差s₀和所述背景估计值b₀作为初值代入所述计算模型进行迭代计算直至与观测值的差值在误差阈值之内，将此时的所述信号亮度的估计值A₀记作信号亮度的精确值A，所述提取到的图像中像素点的X坐标C_x记作像素点精度的X坐标x₀、所述提取到的图像中像素点的Y坐标C_y记作像素点精度的Y坐标y₀、所述荧光点扩散函数的估计标准差s₀记作荧光点扩散函数精确的标准差s，所述背景估计值b₀记作背景精确值b。