一种面向大视场高分辨率显微成像的图像迭代重构方法与流程

本发明属于光学显微成像技术，特别是一种面向大视场高分辨率显微成像的图像迭代重构方法。

背景技术：

在显微成像领域，更高的分辨率一直是追求的目标，但是在提高分辨率的同时存在一个关键性问题，那就是并没有随分辨率一起提高的显微镜的空间带宽积。从成像系统角度看，为了实现高分辨率，必须增加显微物镜的数值孔径，但空间分辨率的提高与视场的扩大往往是一对难以调和的矛盾。简言之，就是在低倍镜下可以看到被检物体的全貌，换成高倍物镜时，就只能看到被检物体的很小一部份。目前，为了突破分辨率与视场大小难以同时兼顾的矛盾，常见的方法是采用常规显微镜系统配合高精度机械扫描和后期空域图像拼接方法将多个小视场高分辨率图像拼接融合生成一幅大视场高分辨率图像([1]2013205777012，适用于结核杆菌抗酸染色图像拼接的装置)。但是由于引入了机械移动装置，所以系统成像时的稳定性和成像速度又成为一对难以调和的矛盾，提高扫描速度必将影响成像稳定性。所以，想要突破分辨率与视场大小难以同时兼顾的矛盾又不引入了机械移动装置，必须采用近年来提出的计算成像的方法，比如基于合成孔径的成像方法。

基于合成孔径成像原理的扫描成像方法最早是由Hoppe为了研究晶体结构所提出的，并通过研究晶体和非晶体的扫描透射电子衍射显微成像，验证了此方法的有效性。Rodenburg和Faulkner等结合相位恢复算法将此方法多次改进，目前这种成像方法已在可见光域、X射线、电子显微镜等不同波段得到了实验证实，并发展出若干种技术以提高成像质量以及分辨率，该技术显示了在大幅面成像和高分辨成像方面的巨大潜力。传统的合成孔径成像技术是通过移动一个全透的小孔(或待测样品本身)使入射平面波照射到待测样品的不同部位，即由小孔控制照明光束尺寸、几何形状及位置，并利用由此得到的一系列衍射强度图样重构出待测样品的振幅与位相信息([2]王雅丽等人.可见光域叠层成像中照明光束的关键参量研究[J].物理学报,2013,Vol.62,No.6.064206-1-9)。合成孔径成像术的关键在于：每次照射待测样品的一个“子孔径”也就是待测样品的某一部分时，都要和至少一个其他的“子孔径”发生交叠。这样就可建立一种重构算法，在分别重构每“子孔径”的复振幅时也要同时满足其他“子孔径”衍射分布的约束，使得最后的待测样品的整体复振幅信息是所有“子孔径”的共同解，从而由各个“子孔径”拼接合成一幅大视场高分辨率的待测样品的图像。合成孔径成像术可以说是一种稳健而简约的显微成像技术，但其一直缺乏一种稳健的图像重构算法。特别是当所采集的图像信噪比较低时，往往难以得到理想的重构图像，因此如何提高重构质量与信噪比成为了合成孔径成像技术必须克服的一个技术难题。

技术实现要素：

本发明的目的在于提供一种面向大视场高分辨率显微成像的图像迭代重构方法，以解决显微系统分辨率与视场大小难以同时兼顾的矛盾，提升显微系统的空间带宽积。

实现本发明目的的技术解决方案为：一种面向大视场高分辨率显微成像的图像迭代重构方法，步骤如下：

步骤一，图像采集：LED阵列作为显微镜的照明光源，顺次点亮其中每一个LED元素，照射样品后采集相对应的图像；

步骤二，初始化：利用LED阵列中位于中心的LED元素照射样品所拍摄到的低分辨率图像来初始化高分辨率图像的振幅与相位；

步骤三，迭代重构：采用增量梯度法将所采集的每一幅图像在频域中逐一进行合成孔径运算；

步骤四，增量梯度迭代系数更新：以代价函数值为判据对增量梯度迭代系数进行更新；

步骤五，停止迭代判断：当增量梯度迭代系数小于一个给定的阈值时，停止迭代，此时的高分辨率图像的振幅与相位就是最终得到的大视场高分辨率显微图像。

本发明还可以在步骤一中的照明光源分别为红光、绿光或蓝光的单色光来照明物体，对于每种照明波长分别采用上述五个步骤进行图像重构，将各自重构的三组图像分别作为最终的真彩色图像的红、绿、蓝分量合成得到重构真彩色图像。

本发明与现有技术相比，其显著优点：(1)利用LED阵列产生不同角度的入射光照射待测样品，再用相机拍摄多幅低分辨率图像，最后采用频域合成孔径技术将多幅大视场低分辨率图像合成一幅大视场高分辨率图像，无需任何机械扫描装置，系统成像稳定性高。与常规显微镜系统的高精度机械扫描和后期空域图像拼接方法相比，大大降低了系统的成本。(2)无需复杂的参数调节，并对采集图像中的噪声具有很强的抵御能力，能够非常稳定并且准确地重建出大视场高分辨率图像。

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

附图说明

图1是基于可编程LED阵列的显微镜光路示意图。

图2是LED阵列中每个像素点的坐标系示意图。

图3为本发明面向大视场高分辨率显微成像的图像迭代重构方法的流程示意图。

图4(a)是4倍物镜(数值孔径0.1)下利用LED阵列中位于中心的LED元素照射1951USAF分辨率板样品所拍摄到的低分辨率图像；选取图4(a)方框中的图像直接进行放大，得到图4(b)。从中再选择更小的区域直接进行放大，又得到图4(c)。

图5(a)为采用本发明方法对1951USAF分辨率板重构的大视场高分辨率显微图像；选取图5(a)方框中的图像直接放大，得到图5(b)。从中再选择更小的区域进行放大，又得到图5(c)。

具体实施方式

本发明面向大视场高分辨率显微成像的图像迭代重构方法所依赖的硬件平台是基于可编程LED阵列的显微镜。结合图1，基于可编程LED阵列的显微镜主要包括显微物镜1、样品2、载物台3、LED阵列4组成。LED阵列4作为显微镜的照明光源，其被直接安置在样品载物台3下方，其距离载物台的上表面距离H一般在20-100mm之间，并且LED阵列4的中心处于显微物镜1的光轴上。LED阵列1中包括若干个LED元素5，它们规则排布形成一个二维阵列。其中单个LED元素均为红绿蓝三色LED，其典型波长为红光635nm、绿光525nm和蓝光475nm。每个LED元素之间中心间距d典型值1-10mm。LED阵列4并不需要进行单独加工，在市场上可直接购置。其包含呈阵列排列的一组多个LED，这些LED通过固定基板实现物理与电路连接，如表1给出了一个市面上可购置的LED阵列的产品参数。在此LED阵列中，LED元素共有32行、32列，一共1024个，单个LED的亮度在2000cd/m²以上。

表1 LED阵列的物理参数

LED阵列1中每个LED元素均可通过实现单独点亮，点亮LED元素的具体方法为现有技术，实现电路可以采用(但不限于)单片机、ARM、或者可编程逻辑器件等现有技术即可实现，具体实现方法可参考相关文献，如郭宝增,邓淳苗：基于FPGA的LED显示屏控制系统设计[J].液晶与显示,2010,25(3):424-428。

本发明面向大视场高分辨率显微成像的图像迭代重构方法实施之前，首先对LED阵列4中每LED元素所对应的照明光的空间频率进行标记，具体方法如下：结合图2，建立坐标系。其中矩形区域代表LED阵列4有效区域，坐标原点位于成像系统的光轴中央。对于任意一个LED元素P，其位置坐标为(P_x,P_y)，首先计算该LED元素所对应的照明光的空间频率矢量：

$u_i=(u_x,u_y)=\frac{2\mathrm{\π}}{\mathrm{\λ}}(\frac{P_x}{\sqrt{P_x^2+P_y^2+H^2}},\frac{P_y}{\sqrt{P_x^2+P_y^2+H^2}})$

类似地，针对LED阵列4中每个LED元素，都可以计算得到其相对应的空间频率矢量，记作u_i，下标i＝1,2,...,N其中N为LED阵列4中LED元素的总个数，λ为单色照明光(红、绿或蓝)的波长。u_i为第i个LED的照明光空间频率矢量。

结合图3，基于上述技术，本发明面向大视场高分辨率显微成像的图像迭代重构方法，其步骤包括：

步骤一：图像采集。LED阵列4作为显微镜的照明光源，顺次点亮其中每一个LED像素，照射样品后采集相对应的图像。由于整个LED阵列中共包含N个LED像素，那么共计拍摄N幅低分辨率图像，记作I_i(r)，i＝1,2,...,N，r为实空间的二维坐标r＝(x,y)。理想情况下，这些拍摄到的低分辨率图像I_i(r)可以被表示为

其中为二维傅里叶变换，u为与r所对应的空间频率坐标，为二维逆傅里叶变换，O(u)为待求的高分辨率物函数的傅里叶变换，即P(u)为显微镜的孔径函数。可以看出，物函数的傅里叶变换首先被第i个LED的照明光空间频率矢量u_i所移动，然后受显微镜的孔径函数P(u)限制后，经过傅里叶逆变换后在实空间形成低分辨率图像I_i(r)，这也就对应着相机所拍摄到的N幅图像I_i(r)，i＝1,2,...,N。

求解高分辨率物函数的傅里叶变换O(u)实质上等价于求解如下的最优化问题

ε为代价函数，其衡量所拍摄到的图像与计算得到的图像之间差异的大小。上述公式表明，待求的高分辨率物函数的傅里叶变换O(u)是令代价函数ε最小化的解。以下为了便于叙述，将上述公式与变量全部表示为向量形式：代表将拍摄到的低分辨率图像I_i(r)按像素顺序排布成的列向量，该列向量中元素个数与图像的像素数相等，均为M。I_i为低分辨率图像向量。代表高分辨率物函数的傅里叶变换O(u)按像素顺序排布成的列向量，该列向量中元素个数为L，O为高分辨率物函数的傅里叶变换向量。由于重构图像的分辨率往往远高于原始图像的分辨率，所以L>>M。P_i是一个由孔径函数P(u)所决定的M×L的矩阵，称之为孔径函数矩阵，其作用就是在整个高分辨率物函数的傅里叶变换向量的M个原素中抽取由孔径所限制的L个元素。通过上述向量表示后，代价函数ε可以被表示为

$\mathrm{\ϵ}=\underset{i}{\Σ}\left|\right|\sqrt{I_i}-\left|F^{-1}{P}_{i}O\right||{|}^2=\underset{i}{\Σ}||\sqrt{I_i}-|g_i\left|\right|{|}^2$

其中||·||代表欧几里得范数，F是离散傅里叶变换的矩阵表达形式，g_i≡F^-1P_iO.根据Parseval定理，代价函数ε可以被化简为

$\mathrm{\ϵ}=\underset{i}{\Σ}\left|\right|{\mathrm{\Ψ}}_i-P_{i}O|{|}^2$

其中

${\mathrm{\Ψ}}_i={\mathrm{\Π}}_{m\→\sqrt{I_i}}\left(P_{i}O\right)$

代表更新后的子频谱；代表实空间模投影操作

${\mathrm{\Π}}_{m\→\sqrt{I_i}}\left(P_{i}O\right)=FDiag\left(g_i\right|g_i{|}^{-1})\sqrt{I_i}$

其代表将gi的幅度部分采用测量值所替换，而保持其相位部分不变。Diag()代表对角矩阵。

步骤二：初始化。利用LED阵列中位于中心的LED元素照射样品所拍摄到的低分辨率图像来初始化高分辨率图像的振幅与相位，如下式所示：

$O_0=F\sqrt{\↑I_{center}}$

式中，I_center为LED阵列4中位于中心的LED元素照射样品所拍摄到的低分辨率图像向量；↑代表图像上采样，即将图像由原始的M像素插值为L个像素。这里的插值方法可直接采用现有技术，如最邻近插值法、双线性插值法或双三次插值法等等，插值方法的选取并不会影响最终的重构结果。增量梯度迭代系数α⁰初始化为α⁰＝1，代价函数值ε⁰初始化为这里O₀的下标0代表内循环迭代次数，α⁰与ε⁰的上标0代表外循环迭代次数。

步骤三：迭代重构。采用增量梯度法将所采集的每一幅图像在频域中逐一进行合成孔径运算。为了最小化代价函数ε，将代价函数写为增量形式，即：

ε_i＝||Ψ_i-P_iO||²

因此，代价函数的增量梯度可以表示为

$\▿{\mathrm{\ϵ}}_i=-{P_i}^{*}({\mathrm{\Ψ}}_i-P_{i}O)$

然后采用增量梯度法对高分辨率物函数的傅里叶变换向量进行更新，这里主要可细分为三个子步骤：

第一步，采用增量梯度法对高分辨率物函数的傅里叶变换向量进行如下更新：

$O_{i+1}^k=O_i^k+{\mathrm{\α}}^k{W}_i\▿{\mathrm{\ϵ}}_i^k=O_i^k-{\mathrm{\α}}^k{W}_i{P_i}^{*}({\mathrm{\Ψ}}_i^k-P_i{O}_i^k)$

其中，α^k为增量梯度迭代系数；代表被I_i更新后的子频谱

${\mathrm{\Ψ}}_i^k={\mathrm{\Π}}_{m\→\sqrt{I_i}}\left(P_i{O}_i^k\right)$

上面的公式中的下标i代表内循环迭代次数，上标k代表外循环迭代次数。W_i为权重矩阵

$W_i=Diag\left(\right|P_{m,i}{|}_{max}^{-2})$

Diag()代表对角矩阵，其对角线元素依次为|P_m,i|_max代表孔径函数模的最大值。W_i的作用是对孔径函数内部的频率成分进行加权，并排除孔径函数外部的频率成分。

第二步，高分辨率物函数的傅里叶变换向量更新完毕之后，内循环迭代次数加1，即i←i+1，切换到下一幅低分辨率图像向量I_i，i＝1,2,...,N，重复第一步。

第三步，当所有图像遍历过一次后，即每幅低分辨率图像向量I_i，i＝1,2,...,N均参与迭代运算直到i＝N，步骤三完成。

步骤四：增量梯度迭代系数更新，分为三个子步骤：

1、计算当前高分辨率物函数的傅里叶变换向量所对应的代价函数值ε^k所对应的误差函数：

${\mathrm{\ϵ}}^k=\underset{i}{\Σ}\left|\right|\sqrt{I_i}-\left|F^{-1}{P}_i{O}_N^k\right||{|}^2$

2、以代价函数值ε^k为判据对增量梯度迭代系数进行更新，更新公式如下：

上述公式表明，增量梯度迭代系数α^k在代价函数值ε^k不再减少时就减小到原来的一半。这样可以保证算法的自适应收敛，并有效提高算法对噪声的鲁棒性。

3、外循环迭代次数加1，即k←k+1，内循环迭代次数清0，即i←0。

这里的核心在于增量梯度迭代系数α^k的趋势在于逐步减小以保证算法的稳定收敛。

步骤五：停止迭代判断。当增量梯度迭代系数α^k小于一个给定的阈值时(建议值为0.001)，停止迭代。此时对迭代得到的高分辨率物函数的傅里叶变换向量进行逆傅里叶变换变换到空域中

$O_{result}=F^{-1}{O}_N^k$

即得到了待测样品大视场高分辨率的振幅分布|O_result|和相位分布angle(O_result)，其中angle()代表取复数的幅角。如果停止条件不满足，则返回步骤三，重新进行一次内循环，直到满足停止条件为止。

本发明物体焦面自动确定方法为：在一种面向大视场高分辨率显微成像的图像迭代重构方法的步骤一：图像采集之前。必须保证载物台3的高度正确，即样品2能够准确的位于显微物镜1的焦面上，对于显微物镜1是无穷远物镜，即保证物体严格位于显微物镜1的前焦面上，此时显微物镜1可以对样品2进行清晰成像。为了保证载物台3的高度正确，按照以下步骤进行载物台高度的调节：

①LED阵列4作为显微镜的照明光源，点亮其中所有LED像素。

②将载物台3的高度调节到最小值z₀，采集相对应的图像

③按如下公式求解所对应的离焦测度函数值

${\mathrm{Laplace}}_{z_0}(x,y)=|2I_{z_0}(x,y)-I_{z_0}(x-step,y)-I_{z_0}(x+step,y)|+|2I_{z_0}(x,y)-I_{z_0}(x,y-step)-I_{z_0}(x,y+step)|$

$L_{z_0}={\mathrm{\Σ\ΣLaplace}}_{z_0}(x,y),{\mathrm{for\Σ\ΣLaplace}}_{z_0}(x,y)\>T$

其中step为梯度求解步长，一般取1-5像素。T为梯度阈值，一般取图像最大灰度值的3％。

④将载物台3的高度提高到z₁，采集相对应的图像

⑤类似地，求解所对应的离焦测度函数

⑥重复上述过程，直至找到所有范围内最小的离焦测度函数值此时所对应的载物台3的高度即是正确值。此时可保证显微物镜1对样品2进行清晰成像。

本发明的步骤三中迭代重构的图像顺序选择：这里有两种备选方案，一种是LED从最中央(光轴上)到外围(远离光轴)顺次进行更新。还有一种是按图像的强度值||I_i||降序进行更新。两种更新顺序都可以得到良好的重构效果。

上述步骤仅仅针对单色照明的情况(如红、绿、蓝等)，若想重构真彩色图像。则每个LED元素需分别采用红光、绿光、蓝光来照明物体，然后对于每种照明波长分别采用上述五个步骤进行图像重构，重构的三组图像分别作为最终的真彩色图像的红、绿、蓝分量合成即可。

为了测试大视场高分辨率显微成像的图像迭代重构方法，我们选取了1951USAF分辨能力测试板进行了成像测试。实验中，使用的LED阵列包含21行21列共441个LED，并利用其产生225个不同角度的照明光，LED之间间距为2.5mm，发出的红光波长为632.8nm。系统所采用的显微物镜的数值孔径为0.1，放大倍率为4x。利用LED阵列中位于中心的LED元素照射样品所拍摄到的低分辨率图像如图4所示。选取图4(a)方框中的图像直接进行放大，得到图4(b)。从中再选择更小的区域直接进行放大，又得到图4(c)。如图可见，最小可分辨的特征为第七组第三个元素，根据1951USAF分辨能力测试板的物理参数(见表2)可知成像系统的原始成像分辨率约为6.2μm/线对。这与成像系统的瑞利衍射极限公式推断的结果吻合良好。而采用本发明重构后的高分辨率图像如图5所示，选取图5(a)方框中的图像直接进行放大，得到图5(b)。从中再选择更小的区域直接进行放大，又得到图5(c)。分辨率板中最小的特征均可以分辨(第九组第三个元素)。通过表2可知，经过合成后的成像系统的分辨率优于1.56μm/线对，即分辨率高于1000um/(645*2)＝0.775um。对比其中图5(c)和图4(c)可以明显看出本发明方法能够在低数值孔径大视场情况下实现大视场高分辨率成像，且重构图像具有良好的信噪比。

表2 1951USAF分辨能力测试板的物理参数