CN113671682B - 一种基于傅里叶叠层显微成像的频域光源位置精确校正方法 - Google Patents

Abstract

本发明基于傅里叶叠层显微成像的频域光源位置精确校正方法中，(1)针对校正位置不精确的问题，对LED阵列中每个LED引入独立的频域位置误差模型进行精确校正；(2)针对校正过程需要参与图像重建过程而带来的时间消耗过长的问题，引入一幅高倍物镜采集的高分辨率强度图像当做参考图像，模拟傅里叶叠层显微成像技术的图像采集过程，生成每个LED对应的虚拟低分辨强度图像。结合虚拟与实际采集的图像，建立损失函数，以二维粒子群算法搜索每个LED对应的最佳频域位置参数。因此，本发明利用引入的高倍物镜下高分辨强度图像与频域位置误差模型实现了基于傅里叶叠层成像原理的LED光源频域位置精确校正，具有实现简易、鲁棒性高、位置精确等优点。

Description

一种基于傅里叶叠层显微成像的频域光源位置精确校正方法

技术领域

本发明涉及定量位相计算显微成像技术领域，具体而言，涉及一种基于傅里叶叠层显微成像的频域光源位置精确校正方法。

背景技术

传统光学显微成像技术以记录光强度信息方式形成图像对比度，为人类打开了显微世界的大门。然而，对于一些透明生物样本，仅凭光强度测量值成像的传统光学显微技术却显得无能为力。1930年，Frits Zernike发明了一种可对位相物体进行高对比成像技术。至此，定量位相成像逐渐成为一种对透明生物样本进行无标记成像的主要方法。不同于传统光学显微成像，定量位相成像根据由样本引入的光程延迟形成高对比度定量图像，并提供对目标特征的客观位相测量值。由于位相与样本的折射率和厚度的乘积为正比，这些测量值可进一步揭示样本的深层结构特征，为后续判断、决策提供相应条件。当前，定量位相成像技术已被广泛研究并应用于生物观察、细胞检测、药物筛选、微量元素研究、精密材料等领域。

为打破传统光学显微成像技术中的一些固有限制条件，如系统分辨率与物镜视场之间的相互制约关系，系统衍射极限条件等，定量位相计算显微成像技术应运而生。作为一种新兴的定量位相计算显微成像技术，傅里叶叠层显微成像技术(Fourier PtychographicMicroscopy，FPM)通过有效结合合成孔径、位相检索、最优化理论等计算成像相关技术，可同时实现宽视场、高分辨、多模态、数字重聚焦、像差校正、定量位相成像等功能。相比于传统数字光学显微成像系统，该技术仅需利用一个可编程的LED阵列替换传统的反射式或主动式光源即可完成系统硬件部分改装，廉价且易于实现。

傅里叶叠层显微成像技术的实现分为两个过程：图像采集过程与图像重建过程。对于图像采集过程，其在低倍物镜下逐个点亮阵列中的LED作为照明光源。当点亮阵列中不同位置的单个LED斜照明一个二维薄样本时，其等效于样本的频谱在频域发生了相对平移效应。这种由不同斜照明引入的位置平移变换使得物镜可收集到原截止频率外的样本频谱信息并在相机面捕获一系列低分辨强度图像。之后在图像重建过程中，这些强度图像将在频域中进行迭代更新样本的频谱信息并最终合成为一张新频谱。经反傅里叶变换后，重建频谱信息即可转化为高分辨的样本复振幅信息。类似于合成孔径理念，傅里叶叠层显微成像技术提升了系统的等效截止频率，从而同时实现了低倍镜下宽视场、高分辨率及定量位相成像。

LED阵列作为傅里叶叠层显微镜中的光源部分，其每个LED的位置参数决定每个子频谱在频域的更新位置，其准确程度直接影响最终图像重建的质量。因此，在实际系统中，对每个LED位置实现精确校正，是实现傅里叶叠层显微成像技术的重要一环。2016年，Sun等人首次建立空域全局误差模型描述LED阵列位置误差，并以模拟退火算法与图像重建过程相结合，完成LED的位置校正。相关实验结果表明该方法有效自动校正了LED位置误差，降低了实验系统搭建中对光源位置精度需求。2017年，Pan等人在此基础上结合嵌入式瞳孔函数恢复算法(Embedded Pupil Function Recovery,EPRY)，进一步提升了校正算法的鲁棒性。然而这些算法都基于空域统一校正模型，以一致性参数表征每个LED与实际位置的误差。不仅如此，使用的模拟退火算法极易陷入局部最优解而非全局最优。综合以上两点分析，现有校正算法依旧存在校正不准确的风险。另一方面，基于模拟退火算法的校正过程必须参与图像重建过程，导致整个校正过程显得十分冗长。

发明内容

本发明旨在提供一种基于傅里叶叠层显微成像的频域光源位置精确校正方法，以解决现有基于傅里叶叠层显微成像技术中的校正算法存在校正位置不精确以及校正过程时间消耗过长的问题。

本发明提供的一种基于傅里叶叠层显微成像的频域光源位置精确校正方法，包括如下步骤：

步骤S100，低倍镜下带位置误差的FPM采集与重建：低倍物镜下逐个点亮可编程LED阵列中的LED，在相机面上捕获关于样本的一系列处于不同LED对应位置的低分辨强度图像；然后利用重建算法对捕获的低分辨强度图像在频域中完成重建，获得重建的带位置误差的高分辨强度与位相复振幅；

步骤S200，高倍物镜下高分辨率强度图像采集与区域配准：使用高倍物镜替换低倍物镜并完成调焦操作后，采集高倍物镜下样本的高分辨强度图像；然后对低分辨强度图像以及高分辨强度图像使用图像配准算法匹配并框选二者图像中表征相同信息区域，对高分辨强度图像中表征相同信息区域经过插值后获得配准后的高分辨强度图像；

步骤S300，LED对应位置的虚拟低分辨强度图像生成：使用带位置误差的高分辨强度与位相复振幅以及配准后的高分辨强度图像组成样本的新频谱图；根据可编程LED阵列中每个LED理想频域位置，从新频谱图中获取对应每个LED理想频域位置的子频谱信息，并对子频谱信息通过傅里叶反变换获得一系列虚拟低分辨强度图像；

步骤S400，阵列中独立LED位置参数的搜索：根据空域中LED对应位置的平移效应等效于频域中孔径位置的平移效应，对每个LED引入频域位置误差模型；由此，将当前LED理想频域位置作为搜索中心，规定搜索区域，并以二维粒子群算法进行最佳位置参数搜索；当所产生的虚拟低分辨强度图像与当前实际采集的低分辨强度图像之间差异最小时，即认为此位置为该LED对应的实际频域位置，完成该LED的位置校正。

进一步地，步骤S100中低倍镜下带位置误差的FPM采集的方法包括如下子步骤：

步骤S101，使用规格为m×n的可编程LED阵列作为光源；假定该可编程LED阵列中相邻LED之间的间隔为d，那么处于第i行、第j列(1≤i≤m,1≤j≤n)的LED理想位置P_i,j表示为：

P_i,j＝(x_i,j,y_i,j)＝(i·d,j·d)

式中，x_i,j与y_i,j表示LED在空域所处水平面的二维坐标分量值；

步骤S102，假定光源位置距样本足够远，那么从LED出射的照明光波可视为单色平面波。因此，使用某个LED照明样本时，对应的波矢量k_i,j表示为：

式中，

与

分别表示波矢量k_i,j在x与y方向上的波矢分量值；λ表示LED辐射的中心波长值；h表示样本距离当前照明样本的LED的垂直距离；

步骤S103，将样本的入射光场的强度值归一化为1，那么从样本出射的频谱表示为：

F{e(r)}＝F{o(r)exp(ik_i,jr)}＝O(k-k_i,j)

式中，F{}表示二维傅里叶变换，e(r)表示从样本出射的复振幅，o(r)表示样本空域的复振幅调制函数，即样本的复振幅函数；o(k)表示样本在频域中的频谱；O(k-k_i,j)表示将样本出射的频谱等效于样本的频谱中心被平移至波矢量k_i,j处；

步骤S104，从样本的出射复振幅在频域被低倍物镜的光瞳函数P(k)低通滤波，该低通滤波过程表示为：

G_i,j(k)＝O(k-k_i,j)P(k)

式中G_i,j(k)表示通过低倍物镜后的样本频谱；

步骤S105，样本频谱传播至像面进行成像，并被相机采集转换为低分辨率强度图像

表示为：

式中，g_i,j(r)表示像面上样本的复振幅，F^-1{}表示二维傅里叶逆变换；

步骤S106，逐个点亮可编程LED阵列中的LED，并对每个点亮的LED执行步骤S102～S105，则捕获关于样本的一系列处于不同LED对应位置的低分辨强度图像

进一步地，步骤S100中的重建方法包括如下子步骤：

步骤S111，频谱初始化，即将可编程LED阵列的中心位置采集的低分辨强度图像插值后的结果当作重建样本频谱的强度部分初始值，并选择零位相作为相关位相的初始值，再结合低倍物镜的光瞳函数，求解重建样本频谱的初始值为：

式中，

表示重建样本频谱的初始值；u与v表示二维频域中的坐标分量；R()表示相关插值过程，I_center(r)表示LED阵列中心位置的LED照明样本时采集的低分辨强度图像，k_center表示LED阵列中心位置的LED照明样本时对应的波矢量；

表示初始位相；

步骤S112，获取子孔径低分辨复振幅，即利用平移效应原理从当前样本频谱中截取相应子孔径样本频谱，然后利用二维傅里叶逆变换，获得对应LED下的子孔径低分辨复振幅估计值：

式中，n表示当前迭代次数，1≤n≤n_max，n_max表示设定的最大迭代次数，

表示更新位置序列(i,j)对应子孔径的样本频谱；(u_i,j,v_i,j)表示位置序列(i,j)在频域中对样本的平移效应，P(u+u_i,j,v+v_i,j)表示平移效应结果；

表示第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；

步骤S113，更新子频谱复振幅估计值：

以位置序列(i,j)实际采集的低分辨强度图像替换子频谱复振幅估计值的振幅部分，并使位相部分保持不变，由此子频谱复振幅估计值的更新公式为：

式中，

表示更新后的第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；

步骤S114，对更新后的第n次迭代对应位置序列(i,j)的子频谱复振幅估计值利用二维傅里叶逆变换即获得对应子孔径频谱信息：

式中，

表示对应

的子孔径频谱信息，F{}表示二维傅里叶变换；

之后，根据位置序列(i,j)对应的子孔径位置更新整体的样本频谱：

式中，

表示更新后的样本频谱；

步骤S115，以步骤S112～S114更新所有位置序列{(i,j)1≤i≤m,1≤j≤n}对应的子孔径频谱信息，得到第n次迭代更新的样本频谱作为重建结果Oⁿ(u,v)；

步骤S116，以步骤S112～S115迭代循环直至重建结果收敛或循环结束后，得到最终重建结果，即重建的带位置误差的样本频谱，表示为

其中，每一次迭代循环的初始值为上一次循环更新后的样本频谱结果；

步骤S117，将重建的带位置误差的样本频谱经过二维傅里叶逆变换，得到重建的带位置误差的样本高分辨复振幅

并将样本高分辨复振幅

分为强度部分与位相部分：

式中，

表示重建的带位置误差的样本高分辨复振幅，

表示重建的带位置误差的样本高分辨复振幅中的强度部分，

表示重建的带位置误差的样本高分辨复振幅中的位相部分。

进一步地，步骤S116中判断重建结果收敛的方法为：

对于第n次迭代循环对应的损失评判函数为：

式中，Eⁿ取得最小值或(Eⁿ-E^n-1)/E^n-1小于设定阈值，则认为重建结果已经收敛。

进一步地，步骤S200包括如下子步骤：

步骤S201，通过相关匹配操作，进行低分辨强度图像与高分辨强度图像显示区域的粗匹配，获得二者图像的粗匹配区域；

步骤S202，对粗匹配区域进行二值化操作，提取用于描述区域的边界中上、下、左、右的四点坐标，以此四点框选区域即认定为表征相同信息区域；

步骤S203，以FPM重建图像分辨率为基准，对高分辨强度图像中表征相同信息区域经过插值后获得配准后的高分辨强度图像I^hr。

进一步地，步骤S300包括如下子步骤：

步骤S301，使用配准后的高分辨采集图像I^hr替换带位置误差的高分辨强度与位相复振幅中的强度部分

位相部分保持不变，经二维傅里叶变换后获得虚拟频谱O^vir：

步骤S302，仿照步骤S100中FPM采集过程，根据位置序列(i,j)对应的理想频域位置(u+u_i,j,v+v_i,j)，进行子孔径采样处理：

式中，

表示位置序列(i,j)虚拟的采样子孔径频谱；

步骤S303，将虚拟的采样子孔径频谱经二维傅里叶变换，获得对应的虚拟低分辨强度图

步骤S304，重复步骤S302～S303直至所有位置序列LED采样完毕，生成一系列虚拟低分辨强度图像

进一步地，步骤S400包括如下子步骤：

步骤S401，对每个位置序列(i,j)引入独立的频域位置误差模型(Δu_i,j,Δv_i,j)，那么，引入频域位置误差模型后对应的虚拟低分辨强度图像表示为：

式中，Δu_i,j与Δv_i,j分别表示二维频域中位置序列(i,j)的位置误差分量；

步骤S402，使用二维相关函数Corr2评价位置序列(i,j)对应的虚拟低分辨强度图像与步骤S100实际采集的低分辨强度图像之间的差异性：

式中，

表示位置序列(i,j)对应的虚拟低分辨强度图像的均值，

表示位置序列(i,j)对应的步骤S100实际采集的低分辨强度图像的均值；

步骤S403，使用二维粒子群算法(PSO)搜索位置序列(i,j)在频域中的最佳频域位置参数：

(Δu_i,j,Δv_i,j)＝argmin(fitness)

fitness＝1-Corr2for(i,j)

式中，fitness表示描述粒子群的适应函数，即损失函数；当全体粒子的适应函数值收敛于最小值时，即表示所产生的虚拟低分辨强度图像与当前实际采集的低分辨强度图像之间差异最小，此时的(Δu_i,j,Δv_i,j)为位置序列(i,j)的频域校正参数；然后使用频域校正参数对位置序列(i,j)的频域位置参数进行校正，得到校正后的最佳频域位置参数；

进一步地，步骤S404还包括：

步骤S404，使用校正后的最佳频域位置参数再次按照步骤S100的方法进行FPM采集与重建，即获得校正后的重建的带位置误差的高分辨强度与位相复振幅。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

本发明基于傅里叶叠层显微成像的频域光源位置精确校正方法中，(1)针对校正位置不精确的问题，对LED阵列中每个LED引入独立的频域位置误差模型进行精确校正；(2)针对校正过程需要参与图像重建过程而带来的时间消耗过长的问题，引入一幅高倍物镜采集的高分辨率强度图像当做参考图像，模拟傅里叶叠层显微成像技术的图像采集过程，生成每个LED对应的虚拟低分辨强度图像。结合虚拟与实际采集的图像，建立损失函数，以二维粒子群算法搜索每个LED对应的最佳频域位置参数。因此，本发明利用引入的高倍物镜下高分辨强度图像与频域位置误差模型实现了基于傅里叶叠层成像原理的LED光源频域位置精确校正，具有实现简易、鲁棒性高、位置精确等优点。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明实施例的基于傅里叶叠层显微成像的频域光源位置精确校正方法的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例

如图1所示，本实施例提出一种基于傅里叶叠层显微成像的频域光源位置精确校正方法，包括如下步骤：

步骤S100，低倍镜下带位置误差的FPM采集与重建：低倍物镜下逐个点亮可编程LED阵列中的LED，在相机面上捕获关于样本的一系列处于不同LED对应位置的低分辨强度图像；然后利用重建算法对捕获的低分辨强度图像在频域中完成重建，获得样本带位置误差的高分辨强度与位相图像；

具体地：

(1)低倍镜下带位置误差的FPM采集的方法包括如下子步骤：

步骤S101，使用规格为m×n的可编程LED阵列作为光源；假定该可编程LED阵列中相邻LED之间的间隔为d，那么处于第i行、第j列(1≤i≤m,1≤j≤n)的LED理想位置P_i,j表示为：

P_i,j＝(x_i,j,y_i,j)＝(i·d,j·d)

式中，x_i,j与y_i,j表示LED在空域所处水平面的二维坐标分量值；

步骤S102，假定光源位置距样本足够远，那么从LED出射的照明光波可视为单色平面波，因此，使用某个LED照明样本时，对应的波矢量k_i,j表示为：

式中，

与

分别表示波矢量k_i,j在x与y方向上的波矢分量值；λ表示LED辐射的中心波长值；h表示样本距离当前照明样本的LED的垂直距离；

步骤S103，将样本的入射光场的强度值归一化为1，那么从样本出射的频谱表示为：

F{e(r)}＝F{o(r)exp(ik_i,jr)}＝O(k-k_i,j)

式中，F{}表示二维傅里叶变换，e(r)表示从样本出射的复振幅，o(r)表示样本空域的复振幅调制函数，即样本的复振幅函数；o(k)表示样本在频域中的频谱；O(k-k_i,j)表示将样本出射的频谱等效于样本的频谱中心被平移至波矢量k_i,j处；

步骤S104，从样本的出射复振幅在频域被低倍物镜的光瞳函数P(k)低通滤波，该低通滤波过程表示为：

G_i,j(k)＝O(k-k_i,j)P(k)

式中G_i,j(k)表示通过低倍物镜后的样本频谱；

步骤S105，样本频谱传播至像面进行成像，并被相机采集转换为低分辨率强度图像

表示为：

式中，g_i,j(r)表示像面上样本的复振幅，F^-1{}表示二维傅里叶逆变换；

步骤S106，逐个点亮可编程LED阵列中的LED，并对每个点亮的LED执行步骤S102～S106，则捕获关于样本的一系列处于不同LED对应位置的低分辨强度图像

(2)重建方法包括如下子步骤：

步骤S111，频谱初始化：

为了使迭代求解算法有着更快的收敛速度，需要先给定求解重建频谱的初始值。一般选择正入射，将可编程LED阵列的中心位置采集的低分辨强度图像插值后的结果当作重建样本频谱的强度部分初始值，并选择零位相作为相关位相的初始值，再结合低倍物镜的光瞳函数，求解重建样本频谱的初始值为：

式中，

表示重建样本频谱的初始值；u与v表示二维频域中的坐标分量；R()表示相关插值过程，I_center(r)表示LED阵列中心位置的LED照明样本时采集的低分辨强度图像，k_center表示LED阵列中心位置的LED照明样本时对应的波矢量；

表示初始位相，一般取0；

步骤S112，获取子孔径低分辨复振幅估计值：

对于由波矢量k_i,j决定的入射角度，利用平移效应原理从当前样本频谱中截取相应子孔径频谱信息，然后利用二维傅里叶逆变换，获得对应LED下的子孔径低分辨复振幅估计值：

式中，n表示当前迭代次数(1≤n≤n_max)，n_max表示设定的最大迭代次数，

表示更新位置序列(i,j)对应子孔径的样本频谱；(u_i,j,v_i,j)表示位置序列(i,j)在频域中对样本的平移效应，P(u+u_i,j,v+v_i,j)表示平移效应结果；

表示第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；

步骤S113，更新子频谱复振幅估计值：

以位置序列(i,j)实际采集的低分辨强度图像替换子频谱复振幅估计值的振幅部分，并使位相部分保持不变，子频谱复振幅估计值的更新公式为：

步骤S114，更新频谱内的子孔径区域：

对更新后的第n次迭代对应位置序列(i,j)的子频谱复振幅估计值利用二维傅里叶逆变换即获得对应子孔径频谱信息：

之后，根据位置序列(i,j)对应的子孔径位置更新整体的样本频谱：

式中，

表示更新后的样本频谱；该式表示，在进行一次子频谱复振幅后，整体样本频谱只对位置序列(i,j)子孔径部分更新频谱信息，其余部分则保持不变。

步骤S115，以步骤S112～S114更新所有位置序列{(i,j)|1≤i≤m,1≤j≤n}对应的子孔径频谱信息，得到第n次迭代更新的样本频谱作为重建结果Oⁿ(u,v)；

步骤S116，以步骤S112～S115迭代循环直至重建结果收敛或循环结束后，得到最终重建结果，即重建的带位置误差的样本频谱，表示为

其中，每一次迭代循环的初始值为上一次循环更新后的样本频谱结果。

其中判断重建结果收敛的方法为：

对于第n次迭代循环对应的损失评判函数为：

式中，Eⁿ取得最小值或(Eⁿ-E^n-1)/E^n-1小于设定阈值，则认为重建结果已经收敛。

步骤S117，将重建的带位置误差的样本频谱经过二维傅里叶逆变换，得到重建的带位置误差的样本高分辨复振幅

并将样本高分辨复振幅

分为强度部分与位相部分：

式中，～表示带位置误差的相关结果，即

表示重建的带位置误差的样本高分辨复振幅，

表示重建的带位置误差的样本高分辨复振幅中的强度部分，

表示重建的带位置误差的样本高分辨复振幅中的位相部分。

步骤S200，高倍物镜下高分辨率强度图像采集与区域配准：使用高倍物镜替换低倍物镜并完成调焦操作后，采集高倍物镜下样本的高分辨强度图像；然后对低分辨强度图像以及高分辨强度图像使用图像配准算法匹配并框选二者图像中表征相同信息区域，对高分辨强度图像中表征相同信息区域经过插值后获得配准后的高分辨强度图像；

由于视场差异，高倍物镜下采集的高分辨率强度图像的显示区域为低倍物镜下的捕获的低分辨强度图像的显示区域的一个子区域，为匹配二者图像的视场并表征相同信息区域，需要对高分辨强度图像进行配准，具体包括如下子步骤：

步骤S201，通过相关匹配操作，进行低分辨强度图像与高分辨强度图像显示区域的粗匹配，获得二者图像的粗匹配区域，也即二者图像的大致区域范围；

步骤S202，对粗匹配区域进行二值化操作，提取用于描述区域的边界中上、下、左、右的四点坐标，以此四点框选区域即认定为表征相同信息区域；

步骤S203，以FPM重建图像分辨率为基准，对高分辨强度图像中表征相同信息区域经过插值后获得配准后的高分辨强度图像I^hr。可见，由于本发明是对表征相同信息区域进行插值，在之后的位置校正过程中，也仅对高分辨强度图像中表征相同信息区域内信息作相关处理。

步骤S300，LED对应位置的虚拟低分辨强度图像生成：使用带位置误差的高分辨强度与位相复振幅以及配准后的高分辨强度图像组成样本的新频谱图；根据可编程LED阵列中每个LED理想频域位置，从新频谱图中获取对应每个LED理想频域位置的子频谱信息，并对子频谱信息通过傅里叶反变换获得一系列虚拟低分辨强度图像；具体包括如下子步骤：

步骤S301，使用配准后的高分辨采集图像I^hr替换带位置误差的高分辨强度与位相复振幅中的强度部分

位相部分保持不变，经二维傅里叶变换后获得虚拟频谱O^vir：

步骤S302，仿照步骤S100中FPM采集过程，根据位置序列(i,j)对应的理想频域位置(u+u_i,j,v+v_i,j)，进行子孔径采样处理：

式中，

表示位置序列(i,j)虚拟的采样子孔径频谱；

步骤S303，将虚拟的采样子孔径频谱经二维傅里叶变换，获得对应的虚拟低分辨强度图

步骤S304，重复步骤S302～S303直至所有位置序列LED采样完毕，生成一系列虚拟低分辨强度图像

步骤S400，阵列中独立LED位置参数的搜索：根据空域中LED对应位置的平移效应等效于频域中孔径位置的平移效应，对每个LED引入频域位置误差模型；由此，将当前LED理想频域位置作为搜索中心，规定搜索区域，并以二维粒子群算法进行最佳位置参数搜索；当所产生的虚拟低分辨强度图像与当前实际采集的低分辨强度图像之间差异最小时，即认为此位置为该LED对应的实际频域位置，完成该LED的位置校正；具体包括如下子步骤：

步骤S401，为了对每个LED进行位置误差精确校正，本发明对每个位置序列(i,j)引入独立的频域位置误差模型(Δu_i,j,Δv_i,j)，那么，引入频域位置误差模型后对应的虚拟低分辨强度图像表示为：

式中，Δu_i,j与Δv_i,j分别表示二维频域中位置序列(i,j)的位置误差分量；

步骤S402，使用二维相关函数Corr2评价位置序列(i,j)对应的虚拟低分辨强度图像与步骤S100实际采集的低分辨强度图像之间的差异性：

式中，

表示位置序列(i,j)对应的虚拟低分辨强度图像的均值，

表示位置序列(i,j)对应的步骤S100实际采集的低分辨强度图像的均值。

步骤S403，使用二维粒子群算法(PSO)搜索位置序列(i,j)在频域中的最佳频域位置参数：

(Δu_i,j,Δv_i,j)＝argmin(fitness)

fitness＝1-Corr2 for(i,j)

式中，fitness表示描述粒子群的适应函数，即损失函数；当全体粒子的适应函数值收敛于最小值时，即表示所产生的虚拟低分辨强度图像与当前实际采集的低分辨强度图像之间差异最小，此时的(Δu_i,j,Δv_i,j)为位置序列(i,j)的频域校正参数；然后使用频域校正参数对位置序列(i,j)的频域位置参数(二维频域中的坐标分量)进行校正，得到校正后的最佳频域位置参数；

步骤S404，使用校正后的最佳频域位置参数再次按照步骤S100的方法进行FPM采集与重建，即获得校正后的重建的带位置误差的高分辨强度与位相复振幅。

至此，本发明基于傅里叶叠层显微成像的频域光源位置精确校正方法中，(1)针对校正位置不精确的问题，对LED阵列中每个LED引入独立的频域位置误差模型进行精确校正；(2)针对校正过程需要参与图像重建过程而带来的时间消耗过长的问题，引入一幅高倍物镜采集的高分辨率强度图像当做参考图像，模拟傅里叶叠层显微成像技术的图像采集过程，生成每个LED对应的虚拟低分辨强度图像。结合虚拟与实际采集的图像，建立损失函数，以二维粒子群算法搜索每个LED对应的最佳频域位置参数。因此，本发明利用引入的高倍物镜下高分辨强度图像与频域位置误差模型实现了基于傅里叶叠层成像原理的LED光源频域位置精确校正，具有实现简易、鲁棒性高、位置精确等优点。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于傅里叶叠层显微成像的频域光源位置精确校正方法，其特征在于，包括如下步骤：

步骤S100，低倍镜下带位置误差的FPM采集与重建：低倍物镜下逐个点亮可编程LED阵列中的LED，在相机面上捕获关于样本的一系列处于不同LED对应位置的低分辨强度图像；然后利用重建算法对捕获的低分辨强度图像在频域中完成重建，获得重建的带位置误差的高分辨强度与位相复振幅；

步骤S200，高倍物镜下高分辨率强度图像采集与区域配准：使用高倍物镜替换低倍物镜并完成调焦操作后，采集高倍物镜下样本的高分辨强度图像；然后对低分辨强度图像以及高分辨强度图像使用图像配准算法匹配并框选二者图像中表征相同信息区域，对高分辨强度图像中表征相同信息区域经过插值后获得配准后的高分辨强度图像；

步骤S300，LED对应位置的虚拟低分辨强度图像生成：使用带位置误差的高分辨强度与位相复振幅以及配准后的高分辨强度图像组成样本的新频谱图；根据可编程LED阵列中每个LED理想频域位置，从新频谱图中获取对应每个LED理想频域位置的子频谱信息，并对子频谱信息通过傅里叶反变换获得一系列虚拟低分辨强度图像；

步骤S400，阵列中独立LED位置参数的搜索：根据空域中LED对应位置的平移效应等效于频域中孔径位置的平移效应，对每个LED引入频域位置误差模型；由此，将当前LED理想频域位置作为搜索中心，规定搜索区域，并以二维粒子群算法进行最佳位置参数搜索；当所产生的虚拟低分辨强度图像与当前实际采集的低分辨强度图像之间差异最小时，即认为此位置为该LED对应的实际频域位置，完成该LED的位置校正；

步骤S100中低倍镜下带位置误差的FPM采集的方法包括如下子步骤：

步骤S101，使用规格为m×n的可编程LED阵列作为光源；假定该可编程LED阵列中相邻LED之间的间隔为d，那么处于第i行、第j列(1≤i≤m,1≤j≤n)的LED理想位置P_i,j表示为：

P_i,j＝(x_i,j,y_i,j)＝(i·d,j·d)

式中，x_i,j与y_i,j表示LED在空域所处水平面的二维坐标分量值；

步骤S102，假定光源位置距样本足够远，那么从LED出射的照明光波可视为单色平面波，因此，使用某个LED照明样本时，对应的波矢量k_i,j表示为：

式中，

与

分别表示波矢量k_i,j在x与y方向上的波矢分量值；λ表示LED辐射的中心波长值；h表示样本距离当前照明样本的LED的垂直距离；

步骤S103，将样本的入射光场的强度值归一化为1，那么从样本出射的频谱表示为：

F{e(r)}＝F{o(r)exp(ik_i,jr)}＝O(k-k_i,j)

式中，F{}表示二维傅里叶变换，e(r)表示从样本出射的复振幅，o(r)表示样本空域的复振幅调制函数，即样本的复振幅函数；o(k)表示样本在频域中的频谱；O(k-k_i,j)表示将样本出射的频谱等效于样本的频谱中心被平移至波矢量k_i,j处；

步骤S104，从样本的出射复振幅在频域被低倍物镜的光瞳函数P(k)低通滤波，该低通滤波过程表示为：

G_i,j(k)＝O(k-k_i,j)P(k)

式中G_i,j(k)表示通过低倍物镜后的样本频谱；

步骤S105，样本频谱传播至像面进行成像，并被相机采集转换为低分辨率强度图像

表示为：

式中，g_i,j(r)表示像面上样本的复振幅，F^-1{}表示二维傅里叶逆变换；

步骤S106，逐个点亮可编程LED阵列中的LED，并对每个点亮的LED执行步骤S102～S105，则捕获关于样本的一系列处于不同LED对应位置的低分辨强度图像

步骤S100中的重建方法包括如下子步骤：

步骤S111，频谱初始化，即将可编程LED阵列的中心位置采集的低分辨强度图像插值后的结果当作重建样本频谱的强度部分初始值，并选择零位相作为相关位相的初始值，再结合低倍物镜的光瞳函数，求解重建样本频谱的初始值为：

式中，

表示重建样本频谱的初始值；u与v表示二维频域中的坐标分量；R()表示相关插值过程，I_center(r)表示LED阵列中心位置的LED照明样本时采集的低分辨强度图像，k_center表示LED阵列中心位置的LED照明样本时对应的波矢量；

表示初始位相；

步骤S112，获取子孔径低分辨复振幅，即利用平移效应原理从当前样本频谱中截取相应子孔径样本频谱，然后利用二维傅里叶逆变换，获得对应LED下的子孔径低分辨复振幅估计值：

式中，n表示当前迭代次数，1≤n≤n_max，n_max表示设定的最大迭代次数，

表示更新位置序列(i,j)对应子孔径的样本频谱；(u_i,j,v_i,j)表示位置序列(i,j)在频域中对样本的平移效应，P(u+u_i,j,v+v_i,j)表示平移效应结果；

表示第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；

步骤S113，更新子频谱复振幅估计值：

以位置序列(i,j)实际采集的低分辨强度图像替换子频谱复振幅估计值的振幅部分，并使位相部分保持不变，由此子频谱复振幅估计值的更新公式为：

式中，

表示更新后的第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；

步骤S114，对更新后的第n次迭代对应位置序列(i,j)的子频谱复振幅估计值利用二维傅里叶逆变换即获得对应子孔径频谱信息：

式中，

表示对应

的子孔径频谱信息，F{}表示二维傅里叶变换；

之后，根据位置序列(i,j)对应的子孔径位置更新整体的样本频谱：

式中，

表示更新后的样本频谱；

步骤S115，以步骤S112～S114更新所有位置序列{(i,j)|1≤i≤m,1≤j≤n}对应的子孔径频谱信息，得到第n次迭代更新的样本频谱作为重建结果Oⁿ(u,v)；

步骤S116，以步骤S112～S115迭代循环直至重建结果收敛或循环结束后，得到最终重建结果，即重建的带位置误差的样本频谱，表示为

其中，每一次迭代循环的初始值为上一次循环更新后的样本频谱结果；

步骤S117，将重建的带位置误差的样本频谱经过二维傅里叶逆变换，得到重建的带位置误差的样本高分辨复振幅

并将样本高分辨复振幅

分为强度部分与位相部分：

式中，

表示重建的带位置误差的样本高分辨复振幅，

表示重建的带位置误差的样本高分辨复振幅中的强度部分，

表示重建的带位置误差的样本高分辨复振幅中的位相部分；

步骤S116中判断重建结果收敛的方法为：

对于第n次迭代循环对应的损失评判函数为：

式中，Eⁿ取得最小值或(Eⁿ-E^n-1)/E^n-1小于设定阈值，则认为重建结果已经收敛；

步骤S200包括如下子步骤：

步骤S201，通过相关匹配操作，进行低分辨强度图像与高分辨强度图像显示区域的粗匹配，获得二者图像的粗匹配区域；

步骤S202，对粗匹配区域进行二值化操作，提取用于描述区域的边界中上、下、左、右的四点坐标，以此四点框选区域即认定为表征相同信息区域；

步骤S203，以FPM重建图像分辨率为基准，对高分辨强度图像中表征相同信息区域经过插值后获得配准后的高分辨强度图像I^hr；

步骤S300包括如下子步骤：

步骤S301，使用配准后的高分辨采集图像I^hr替换带位置误差的高分辨强度与位相复振幅中的强度部分

位相部分保持不变，经二维傅里叶变换后获得虚拟频谱O^vir：

步骤S302，仿照步骤S100中FPM采集过程，根据位置序列(i,j)对应的理想频域位置(u+u_i,j,v+v_i,j)，进行子孔径采样处理：

式中，

表示位置序列(i,j)虚拟的采样子孔径频谱；

步骤S303，将虚拟的采样子孔径频谱经二维傅里叶变换，获得对应的虚拟低分辨强度图

步骤S304，重复步骤S302～S303直至所有位置序列LED采样完毕，生成一系列虚拟低分辨强度图像

步骤S400包括如下子步骤：

步骤S401，对每个位置序列(i,j)引入独立的频域位置误差模型(Δu_i,j,Δv_i,j)，那么，引入频域位置误差模型后对应的虚拟低分辨强度图像表示为：

式中，Δu_i,j与Δv_i,j分别表示二维频域中位置序列(i,j)的位置误差分量；

步骤S402，使用二维相关函数Corr2评价位置序列(i,j)对应的虚拟低分辨强度图像与步骤S100实际采集的低分辨强度图像之间的差异性：

式中，

表示位置序列(i,j)对应的虚拟低分辨强度图像的均值，

表示位置序列(i,j)对应的步骤S100实际采集的低分辨强度图像的均值；

步骤S403，使用二维粒子群算法(PSO)搜索位置序列(i,j)在频域中的最佳频域位置参数：

(Δu_i,j,Δv_i,j)＝argmin(fitness)

fitness＝1-Corr2for(i,j)

式中，fitness表示描述粒子群的适应函数，即损失函数；当全体粒子的适应函数值收敛于最小值时，即表示所产生的虚拟低分辨强度图像与当前实际采集的低分辨强度图像之间差异最小，此时的(Δu_i,j,Δv_i,j)为位置序列(i,j)的频域校正参数；然后使用频域校正参数对位置序列(i,j)的频域位置参数进行校正，得到校正后的最佳频域位置参数。

2.根据权利要求1所述的基于傅里叶叠层显微成像的频域光源位置精确校正方法，其特征在于，步骤S404还包括：

步骤S404，使用校正后的最佳频域位置参数再次按照步骤S100的方法进行FPM采集与重建，即获得校正后的重建的带位置误差的高分辨强度与位相复振幅。

Images