CN106127767A - 一种针对傅立叶叠层显微成像技术的位置校正方法 - Google Patents

Abstract

本发明公开了一种针对傅立叶叠层显微成像技术的位置校正方法，首先摄一组低分辨率图像，对物体的高分辨率频谱进行初始化，并初始化迭代次数j＝1；计算第j次迭代的更新范围S_j；对更新范围S_j内的所有图像进行更新，更新前利用模拟退火法校正每幅图像对应的频谱孔径位置；第j次迭代完成后，利用非线性回归法更新LED阵列的位置参数，重新对物体的高分辨率频谱进行初始化；j＝j+1，若更新范围S_j不包含所有图像，回到迭代步骤，当更新范围S_j包含所有图像以后，执行下一步；继续迭代至少3次，每次迭代完不进行频谱初始化，最终获得物体的高分辨率光强和相位图。本发明避免了LED阵列定位误差对重构结果的影响，提高了傅立叶叠层显微成像技术重构的图像质量。

Description

一种针对傅立叶叠层显微成像技术的位置校正方法

技术领域

本发明属于显微成像技术的系统校正领域，特别是一种针对傅立叶叠层显微成像技术的位置校正方法。

背景技术

在显微成像领域，更高的分辨率一直是追求的目标，但是在提高分辨率的同时存在一个关键性问题，那就是并没有随分辨率一起提高的显微镜的空间带宽积，换言之即传统显微镜存在分辨率与视场大小难以同时兼顾的矛盾。因为，传统显微镜使用低倍物镜进行成像时视场大但是分辨率低，而使用高倍物镜进行成像时分辨率提高了但是相应的视场就会缩得很小。目前，为了突破分辨率与视场大小难以同时兼顾的矛盾，常见的方法是采用常规显微镜系统配合高精度机械扫描和后期空域图像拼接方法将多个小视场高分辨率图像拼接融合生成一幅大视场高分辨率图像([1]2013205777012，适用于结核杆菌抗酸染色图像拼接的装置)。但是由于引入了机械移动装置，所以系统成像时的稳定性和成像速度又成为一对难以调和的矛盾，提高扫描速度必将影响成像稳定性。所以，想要突破分辨率与视场大小难以同时兼顾的矛盾又不引入了机械移动装置，必须采用近年来提出的计算成像的方法，比如傅立叶叠层显微成像技术。

傅里叶叠层成像技术是近年来发展出的一种大视场高分辨率定量相位计算显微成像技术([2]谢宗良,马浩统,任戈,等.小孔扫描傅里叶叠层成像的关键参量研究[J].光学学报,2015,35(10):94-102)，该方法整合了相位恢复和合成孔径的概念。与其他相位恢复方法相似，傅里叶叠层成像技术的处理过程也是根据空域中记录的光强信息和频域中某种固定的映射关系来进行交替迭代的，特别的是该技术借用了合成孔径叠层成像的思想。在一个传统的傅里叶叠层成像的系统中，LED阵列中每个LED单元从不同角度照射样品并通过一个低数值孔径的物镜进行成像。由于二维的薄物体被来自不同角度的平面波照明，所以在物镜后焦面上物体的频谱被平移到对应的不同位置。因此，一些本来超出物镜数值孔径的频率成分被平移到物镜数值孔径以内从而能够传递到成像面进行成像。反过来看，不同角度的入射光可等效为在频谱上不同位置的交叠的光瞳函数(子孔径)，每次通过不同位置子孔径的频谱在频域上形成叠层。之后再利用相机拍摄到的一系列低分辨率图像在频域里迭代，依次更新对应的子孔径里的频谱信息，子孔径与子孔径交叠着扩展了频域带宽并恢复出超过物镜空间分辨率限制的高频信息(合成孔径)，最终同时重构出物体的大视场高分辨率光强和相位图像(相位恢复)。这样就实现了使用一个低数值孔径、低放大率的物镜同时获得大视场和高分辨率成像，最终重构的分辨率取决于频域中合成数值孔径的大小。

傅里叶叠层显微成像技术是一种新型的大视场高分辨显微成像技术，但目前其重构质量往往会受到LED阵列定位误差的影响，尤其是调节傅里叶叠层显微成像系统光路时LED单元的空间位置难以非常准确的进行定位，因此如何有效的校正LED阵列的定位误差，使其不影响最终的重构质量就成为了傅里叶叠层显微成像技术必须克服的一个技术难题。

发明内容

本发明的目的在于提供一种针对傅立叶叠层显微成像技术的位置校正方法，以提高傅里叶叠层显微成像重构的图像质量。

实现本发明目的的技术解决方案为：一种针对傅立叶叠层显微成像技术的位置校正方法，步骤如下：

步骤一，拍摄一组低分辨率图像，对物体的高分辨率频谱进行初始化，并初始化迭代次数j＝1；

步骤二，计算第j次迭代的更新范围S_j；

步骤三，对更新范围S_j内的所有图像进行更新，更新前利用模拟退火法校正每幅图像对应的频谱孔径位置；

步骤四，第j次迭代完成后，利用非线性回归法更新LED阵列的位置参数，重新对物体的高分辨率频谱进行初始化；

步骤五，j＝j+1，若更新范围S_j不包含所有图像，回到第二步进行迭代；当更新范围S_j包含所有图像以后，执行第六步；

步骤六，继续迭代至少3次，每次迭代完不进行频谱初始化，最终获得物体的高分辨率光强和相位图。

本发明与现有技术相比，其显著优点：(1)采用模拟退火法校正每个频谱孔径的位置，这样能够有效避免优化过程陷入局部最优，提高了LED阵列位置校正的精度。(2)在每轮迭代完成后都采用非线性回归法计算LED阵列的位置参数，配合更新范围逐步扩大，这样可以提高算法的鲁棒性，降低图像噪声对校正精度的影响。

下面结合附图对本发明作进一步详细描述。

附图说明

图1为本发明针对傅立叶叠层显微成像技术的位置校正方法的流程图。

图2为以USAF分辨率板作为待测样品的傅立叶叠层显微成像结果，图2(a)表示的没有进行位置误差校正的重构结果，图2(b)表示的是使用本方法进行位置误差校正的重构结果。图2(c)表示的是使用本方法校正的频谱孔径中心的位置，其中三角点表示的是校正前错误的频谱孔径中心的位置，圆形点表示的是实际真实的频谱孔径中心的位置，菱形点表示的是使用本方法校正得到的频谱孔径中心的位置。

具体实施方式

本发明首先计算每次迭代的更新范围，然后对更新范围内的所有图像进行更新，更新前利用模拟退火法校正每幅图像对应的频谱孔径位置；每次迭代完成后，利用非线性回归法更新LED阵列的位置参数，重新对物体的高分辨率频谱进行初始化，然后判断更新范围是否包含所有拍摄到图像，当更新范围包含所有图像以后，继续迭代至少3次，每次迭代完不进行频谱初始化，最终获得物体的高分辨率光强和相位图。

结合图1，本发明针对傅立叶叠层显微成像技术的位置校正方法的具体实现步骤如下：

步骤一，拍摄一组低分辨率图像，对物体的高分辨率频谱进行初始化，并初始化迭代次数j＝1。对物体的高分辨率频谱进行初始化的公式为：

其中，O₀为物体初始化的高分辨率频谱，表示傅立叶变换求频谱，表示对一幅图像进行双线性插值，I_0,0为第0行第0列的LED单元垂直照明是拍摄到的低分辨率图像，为第0行第0列的LED单元对应的频谱里的孔径函数，(u_0,0,v_0,0)表示第0行第0列的LED单元对应的频谱里的孔径中心的频域坐标。

步骤二，计算第j次迭代的更新范围S_j。第j次迭代的更新范围S_j的公式为：

S_j＝{(m,n)|m＝-(2j+3),,,(2j+3),n＝-(2j+3),,,(2j+3)}

其中，S_j为第j次迭代的更新范围，(m,n)分别为要更新的图像对应的LED单元的行数和列数。

步骤三，对更新范围S_j内的所有图像进行更新，更新前利用模拟退火法校正每幅图像对应的频谱孔径位置。利用模拟退火法校正每幅图像对应的频谱孔径位置的具体流程为：

①初始化模拟退火搜索次数r，r＝1。

②计算第r次模拟退火搜索得到的第m行第n列LED对应的频谱孔径中心的频域坐标(u_r,m,n,v_r,m,n)

(u_r,m,n,v_r,m,n)＝(u_m,n,v_m,n)+(Δu_r,Δv_r)

其中，(u_m,n,v_m,n)为第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标，(Δu_r,Δv_r)为第r次模拟退火搜索随机的频谱平移距离。

③计算第r次模拟退火选择的第m行第n列LED对应的频谱

$O_{r,j,m,n}^e(u,v)=O_j(u-u_{r,m,n},v-v_{r,m,n})P_j(u,v)$

其中，(u,v)表示频域坐标，P_j(u,v)为第j次迭代的频谱里的孔径函数。

④计算第r次模拟退火选择的第m行第n列LED对应的复振幅

其中，表示逆傅立叶变换。

⑤计算第r次模拟退火搜索的光强均方误差E(r)

$E\left(r\right)=\underset{x,y}{\Σ}{\[{\left|o_{r,j,m,n}^e(x,y)\right|}^2-I_{m,n}^c(x,y)\]}^2$

其中，为第m行第n列的LED单元对应的低分辨率图像，(x,y)表示空间坐标。

⑥判断r是否大于8。若r<8，则r＝r+1并回到子步骤②。若r>＝8，则进行子步骤⑦。这里为了平衡模拟退火搜索的精度和效率，选择模拟退火搜索总次数为8次，若要提高模拟退火搜索的精度，可以适当提高搜索总次数，如要增加模拟退火搜索的效率，可以适当减少搜索总次数。

⑦计算更新后的第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标

s＝argmin[E(r)]

$(u_{m,n}^u,v_{m,n}^u)=(u_{m,n},v_{m,n})+({\mathrm{\&Delta;u}}_s,{\mathrm{\&Delta;v}}_s)$

其中，s表示一共8次模拟退火搜索的8个光强均方误差值中的最小值所对应的搜索次数，argmin[E(r)]表示找到E(r)最小值所对应的r，(Δu_s,Δv_s)为第s次模拟退火搜索随机的频谱平移距离。

步骤四，第j次迭代完成后，利用非线性回归法更新LED阵列的位置参数，重新对物体的高分辨率频谱进行初始化。利用非线性回归法更新LED阵列的位置参数的公式为：

$Q(\mathrm{\&theta;},\mathrm{\&Delta;}x,\mathrm{\&Delta;}y,h)=\underset{m,n}{\&Sigma;}\&lsqb;{(u_{m,n}-u_{m,n}^u)}^2+{(v_{m,n}-v_{m,n}^v)}^2\&rsqb;$

$x_{m,n}^i=d_{LED}\&lsqb;cos\left(\mathrm{\&theta;}\right)m+sin\left(\mathrm{\&theta;}\right)n\&rsqb;+\mathrm{\&Delta;}x$

$y_{m,n}^i=d_{LED}\&lsqb;-sin\left(\mathrm{\&theta;}\right)m+cos\left(\mathrm{\&theta;}\right)n\&rsqb;+\mathrm{\&Delta;}y$

$u_{m,n}=\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{x^o-x_{m,n}^i}{\sqrt{{(x^o-x_{m,n}^i)}^2+{(y^o-y_{m,n}^i)}^2+h^2}}$

$v_{m,n}=\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{y^o-y_{m,n}^i}{\sqrt{{(x^o-x_{m,n}^i)}^2+{(y^o-y_{m,n}^i)}^2+h^2}}$

其中，Q(θ,Δx,Δy,h)为非线性回归法的目标函数，(θ^u,Δx^u,Δy^u,h^u)为更新后的LED阵列的四个位置参数，分别是旋转误差、x方向的平移误差、y方向的平移误差、高度误差，表示进行非线性回归运算，d_LED为LED阵列相邻两个单元之间间距，表示第m行第n列的LED单元的空间位置坐标，(x^o,y^o)为重构小区域中心的空间位置坐标，λ为照明光的波长，(u_m,n,v_m,n)为第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标，为步骤三中获得的更新后的第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标。

步骤五，j＝j+1，若更新范围S_j不包含所有图像，回到第二步进行迭代。当更新范围S_j包含所有图像以后，执行第六步。

步骤六，继续迭代至少3次，每次迭代完不进行频谱初始化，最终获得物体的高分辨率光强和相位图。

通过上述步骤可以看出，本发明采用模拟退火法校正每个频谱孔径的位置，这样能够有效避免优化过程陷入局部最优，提高了LED阵列位置校正的精度。此外，本方法在每轮迭代完成后都采用非线性回归法计算LED阵列的位置参数，配合更新范围逐步扩大，这样可以提高算法的鲁棒性，降低图像噪声对校正精度的影响。

为了测试一种针对傅立叶叠层显微成像技术的位置校正方法的校正效果，我们以USAF分辨率板作为待测样品进行了傅立叶叠层显微成像，图2(a)表示的没有进行位置误差校正的重构结果，图2(b)表示的是使用本方法进行位置误差校正的重构结果。图2(c)表示的是使用本方法校正的频谱孔径中心的位置，其中三角点表示的是校正前错误的频谱孔径中心的位置，圆形点表示的是实际真实的频谱孔径中心的位置，菱形点表示的是使用本方法校正得到的频谱孔径中心的位置。从图2(a)和2(b)中可以看出，经过本方法校正LED阵列位置以后，所有高频细节信息都得到了准确的重构。从图2(c)中可以看出，使用本方法校正得到的频谱孔径中心的位置与实际真实的频谱孔径中心的位置完全重合，说明本方法能够实现高精度的LED阵列定位误差校正，显著提高重构图像质量。

Claims (5)

1.一种针对傅立叶叠层显微成像技术的位置校正方法，其特征在于步骤如下：

步骤一，拍摄一组低分辨率图像，对物体的高分辨率频谱进行初始化，并初始化迭代次数j＝1；

步骤二，计算第j次迭代的更新范围S_j；

步骤三，对更新范围S_j内的所有图像进行更新，更新前利用模拟退火法校正每幅图像对应的频谱孔径位置；

步骤四，第j次迭代完成后，利用非线性回归法更新LED阵列的位置参数，重新对物体的高分辨率频谱进行初始化；

步骤五，j＝j+1，若更新范围S_j不包含所有图像，回到第二步进行迭代；当更新范围S_j包含所有图像以后，执行第六步；

步骤六，继续迭代至少3次，每次迭代完不进行频谱初始化，最终获得物体的高分辨率光强和相位图。

2.根据权利要求1所述的针对傅立叶叠层显微成像技术的位置校正方法，其特征在于步骤一中，对物体的高分辨率频谱进行初始化的公式为：

其中，O₀为物体初始化的高分辨率频谱，表示傅立叶变换求频谱，表示对一幅图像进行双线性插值，I_0,0为第0行第0列的LED单元垂直照明是拍摄到的低分辨率图像，为第0行第0列的LED单元对应的频谱里的孔径函数，(u_0,0,v_0,0)表示第0行第0列的LED单元对应的频谱里的孔径中心的频域坐标。

3.根据权利要求1所述的针对傅立叶叠层显微成像技术的位置校正方法，其特征在于步骤二中，第j次迭代的更新范围S_j的公式为：

S_j＝{(m,n)|m＝-(2j+3),,,(2j+3),n＝-(2j+3),,,(2j+3)}

其中，S_j为第j次迭代的更新范围，(m,n)分别为要更新的图像对应的LED单元的行数和列数。

4.根据权利要求1所述的针对傅立叶叠层显微成像技术的位置校正方法，其特征在于步骤三中，利用模拟退火法校正每幅图像对应的频谱孔径位置的过程为：

①初始化模拟退火搜索次数r，r＝1；

②计算第r次模拟退火搜索得到的第m行第n列LED对应的频谱孔径中心的频域坐标(u_r，m，_n，v_r，m，n)

(u_r，m，n，v_r，m，n)＝(u_m，n，v_m，n)+(Δu_r，Δv_r)

其中，(u_m，n，v_m，n)为第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标，(Δu_r，Δv_r)为第r次模拟退火搜索随机的频谱平移距离；

③计算第r次模拟退火选择的第m行第n列LED对应的频谱

$O_{r,j,m,n}^e(u,v)=O_j(u-u_{r,m,n},v-v_{r,m,n})P_j(u,v)$

其中，(u，v)表示频域坐标，P_j(u，v)为第j次迭代的频谱里的孔径函数；

④计算第r次模拟退火选择的第m行第n列LED对应的复振幅

其中，表示逆傅立叶变换；

⑤计算第r次模拟退火搜索的光强均方误差E(r)

$E\left(r\right)=\underset{x,y}{\&Sigma;}{\&lsqb;{\left|o_{r,j,m,n}^e(x,y)\right|}^2-I_{m,n}^c(x,y)\&rsqb;}^2$

其中，为第m行第n列的LED单元对应的低分辨率图像，(x，y)表示空间坐标；

⑥判断r是否大于8，若r<8，则r＝r+1并回到子步骤②；若r>＝8，则进行子步骤⑦；

⑦计算更新后的第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标

s＝argmin[E(r)]

$(u_{m,n}^u,v_{m,n}^u)=(u_{m,n},v_{m,n})+({\mathrm{\&Delta;u}}_s,{\mathrm{\&Delta;v}}_s)$

其中，s表示模拟退火搜索的光强均方误差值中的最小值所对应的搜索次数，argmin[E(r)]表示找到E(r)最小值所对应的r，(Δu_s，Δv_s)为第s次模拟退火搜索随机的频谱平移距离。

5.根据权利要求1所述的针对傅立叶叠层显微成像技术的位置校正方法，其特征在于步骤四中，利用非线性回归法更新LED阵列的位置参数的公式为：

$Q(\mathrm{\&theta;},\mathrm{\&Delta;}x,\mathrm{\&Delta;}y,h)=\underset{m,n}{\&Sigma;}\&lsqb;{(u_{m,n}-u_{m,n}^u)}^2+{(v_{m,n}-v_{m,n}^v)}^2\&rsqb;$

$x_{m,n}^i=d_{LED}\&lsqb;cos\left(\mathrm{\&theta;}\right)m+sin\left(\mathrm{\&theta;}\right)n\&rsqb;+\mathrm{\&Delta;}x$

$y_{m,n}^i=d_{LED}\&lsqb;-sin\left(\mathrm{\&theta;}\right)m+cos\left(\mathrm{\&theta;}\right)n\&rsqb;+\mathrm{\&Delta;}y$

$u_{m,n}=\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{x^o-x_{m,n}^i}{\sqrt{{(x^o-x_{m,n}^i)}^2+{(y^o-y_{m,n}^i)}^2+h^2}}$

$v_{m,n}=\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{y^o-y_{m,n}^i}{\sqrt{{(x^o-x_{m,n}^i)}^2+{(y^o-y_{m,n}^i)}^2+h^2}}$

其中，Q(θ,Δx,Δy,h)为非线性回归法的目标函数，(θ^u,Δx^u,Δy^u,h^u)为更新后的LED阵列的四个位置参数，分别是旋转误差、x方向的平移误差、y方向的平移误差、高度误差，表示进行非线性回归运算，d_LED为LED阵列相邻两个单元之间间距，表示第m行第n列的LED单元的空间位置坐标，(x^o,y^o)为重构小区域中心的空间位置坐标，λ为照明光的波长，(u_m,n,v_m,n)为第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标，为步骤三中获得的更新后的第m行第n列的LED单元对应的频谱里的孔径中心的频域坐标。