一种数字全息频谱中心亚像素搜索方法

摘要

本发明公开了一种数字全息频谱中心亚像素搜索方法，该方法主要由频谱图采集、+1级频谱筛选、频谱低通滤波、二维曲面拟合、相位恢复及亚像素位移相位校准共五个步骤。本发明针对数字全息显微成像系统，基于二维曲面拟合的方法将离散的频谱矩阵连续化得到精确到亚像素的频谱正中心，并采用高斯函数频谱滤波降噪提高了数字全息的相位恢复效果和精度，对数字全息的高精度，高实时性的无标记观测起到十分重要的作用。

权利要求

1.一种数字全息频谱中心亚像素搜索方法，其特征在于，频谱中心高精度寻址过程如下：步骤一，频谱图采集：利用马赫-曾德尔干涉法采集待测样品的干涉条纹图，通过快速傅里叶变换得到对应的频谱图；

步骤二，+1级频谱筛选：根据物镜的数值孔径和照明光波长计算得到+1级频谱带宽，根据谱宽生成二维矩形窗函数实现对频谱的带通滤波，得到频谱图中载有样品信息的+1级谱；

步骤三，频谱低通滤波：将步骤二得到的+1级频谱利用二维高斯窗函数进行低通滤波，滤掉部分高频成分来降低噪声的影响；

步骤四，二维曲面拟合：通过曲面拟合方法，将经过低通滤波后的+1级离散二维频谱矩阵拟合为连续的二维曲面函数，根据得到的二维曲面函数，利用微分方法求得频谱区域极大值点位置，该位置即为精确到亚像素的频谱中心 ；

步骤五，相位恢复及亚像素位移相位校准：将步骤四得到的频谱中心位置向下取整、 ，将滤波后的+1级谱的 和 置于空矩阵中心并逆傅里叶变换到空间域，根据谱中心的亚像素偏差，利用傅里叶变换性质对恢复结果进行亚像素位移相位补偿，得到相位恢复结果；

步骤四中，以二维Hann窗函数为目标函数利用最小二乘法对+1级谱进行曲面拟合，选择拟合对象函数为Hann窗函数：其中a为方程的待定系数，N为离散化表达的拟合区间，n为像素位置，n=1,2,3,…  ,N，n是自变量，n最大不超过N。

2.根据权利要求1所述的一种数字全息频谱中心亚像素搜索方法，其特征在于，步骤三中，利用二维高斯函数进行+1级频谱滤波，滤除高频的噪声成分达到频率域降噪的效果。

3.根据权利要求1所述的一种数字全息频谱中心亚像素搜索方法，其特征在于，步骤四中，利用微分方法求得极大值点位置，方法是获取连续的二维曲面函数的各点切面，找到水平切面对应的切点，即极大值点，为所求的精度为亚像素的频谱中心。

4.根据权利要求1所述的一种数字全息频谱中心亚像素搜索方法，其特征在于，步骤五中，利用傅里叶变换性质结合高精度频谱中心进行亚像素位移以获得更好的相位恢复结果；由于只存在整数的像素点位置，将亚像素的频谱中心取整后的坐标 置于空矩阵中心，根据傅里叶变换时移和频移性质对恢复结果进行亚像素位移相位补偿，得到相位恢复结果：其中F[…]是傅里叶变换，F-1[…]是傅里叶逆变换，t为时间，w为频率，w0为频移量，f(t)为时域函数，F(w)为频域函数。

说明书

一种数字全息频谱中心亚像素搜索方法

技术领域

[0001] 本发明涉及定量相位成像与光学显微成像，特别的是针对数字全息显微成像提出的一种数字全息频谱中心亚像素搜索方法。

背景技术

[0002] 在生物学、生物医学和病理学研究领域，在细胞和亚细胞水平的细胞三维结构的定量分析是非常必要的。然而目前的光学显微成像系统大多数只能通过染色对二维的平面进行成像且无法对其进行测量，即使有个别大型显微成像设备可以实现对观测物的三维重构，但由于其耗时较长无法对观测物进行实时的动态三维测量。同时通过荧光标记的侵入式手段往往会对细胞的活性造成影响，甚至可能会杀死细胞。无标记成像技术解决了这一个问题，比如微分干涉相衬技术（DIC）。无需任何外源的标记物对细胞进行处理，可以在保持细胞活性的情况下实现长时间成像。微分干涉显微镜通过石英Nomarski棱镜使透过样品的偏振光产生干涉，将样品的相位差异转换为可以观察到的强度差。这种方法虽然可以观察显微结构但也只是一个可视化的观察结果，无法定量测量。

[0003] 数字全息（Digital Holography）显微成像将定量相位成像和光学显微技术相结合，以其高速、无损、高分辨率的成像特点对微观物体实现实时的三维形貌、透明物体厚度/折射率分布测量。同时在微光学测试方面，数字全息技术在测量具有非常紧凑和简单几何形状的小型镜头结构同样具有很大的优势。虽然近年来国内外研究人员提出了许多定量相位成像方法，比如光强传输方程（Transport-of-intensity equation, TIE），差分相衬成像（Differential phase contrast, DPC），傅里叶叠层成像（Fourier ptychographic imaging, FPI）等。在诸多方法中数字全息定量相位显微成像采用相干成像的原理，并采用数字的方法实现相位恢复和三维重构，具有更好的成像效率、重建精度和稳定性，成为细胞研究的极具前景的方法之一。

[0004] 但由于数字全息定量相位恢复的高精度要求，相位恢复的精度和对相位像差的矫正是十分必要的。因此本发明一种数字全息频谱中心亚像素搜索方法。可以不以牺牲运算速度为代价的情况下提升相位恢复精度和成像质量。整数频谱中心进行相位恢复的普通方案下会产生相位面的倾斜，这是由于去一阶像差不彻底导致的，通过曲面拟合的亚像素频谱中心寻址得到更为精确的谱中心，但因其精度为亚像素需要通过亚像素位移对倾斜的相位面进行矫正，获得更好的定量相位恢复结果和成像质量。

发明内容

[0005] 为了提高数字全息定量相位成像的成像质量，来解决数字全息的定量相位成像的成像质量和测量精度问题，实现实时动态、高分辨率、高稳定性的数字全息相位成像，本发明提出了一种数字全息频谱中心亚像素搜索方法。

[0006] 本发明的技术方案如下：一种数字全息频谱中心亚像素搜索方法，频谱中心高精度寻址过程如下：

[0007] 步骤一，频谱图采集：利用马赫-曾德尔干涉法采集待测样品的干涉条纹图，通过快速傅里叶变换得到对应的频谱图；

[0008] 步骤二，+1级频谱筛选：根据物镜的数值孔径和照明光波长计算得到+1级频谱带宽，根据谱宽生成二维矩形窗函数实现对频谱的带通滤波，得到频谱图中载有样品信息的+1级谱；

[0009] 步骤三，频谱低通滤波：将步骤二得到的+1级频谱利用二维高斯窗函数进行低通滤波，滤掉部分高频成分来降低噪声的影响；

[0010] 步骤四，二维曲面拟合：通过曲面拟合方法，将经过低通滤波后的+1级离散二维频谱矩阵拟合为连续的二维曲面函数，根据得到的二维曲面函数，利用微分方法求得频谱区域极大值点位置，该位置即为精确到亚像素的频谱中心 ；

[0011] 步骤五，相位恢复及亚像素位移相位校准：将步骤四得到的频谱中心位置向下取整 、 ，将滤波后的+1级谱的 和 置于空矩阵中心并逆傅里叶变换到空间域，根据谱中心的亚像素偏差，利用傅里叶变换性质对恢复结果进行亚像素位移相位补偿，得到相位恢复结果；

[0012] 步骤四中，以二维Hann窗函数为目标函数利用最小二乘法对+1级谱进行曲面拟合，选择拟合对象函数为Hann窗函数：

[0013]

[0014] 其中a为方程的待定系数，N为离散化表达的拟合区间，n为像素位置，n=1,2,3,… ,N，n是自变量，n最大不超过N。

[0015] 优选的，步骤三中，利用二维高斯函数进行+1级频谱滤波，滤除高频的噪声成分达到频率域降噪的效果。

[0016] 优选的，步骤四中，利用微分方法求得极大值点位置，方法是获取连续的二维曲面函数的各点切面，找到水平切面对应的切点，即极大值点，为所求的精度为亚像素的频谱中心。

[0017] 优选的，步骤五中，利用傅里叶变换性质结合高精度频谱中心进行亚像素位移以获得更好的相位恢复结果；由于只存在整数的像素点位置，将亚像素的频谱中心取整后的坐标 置于空矩阵中心，根据傅里叶变换时移和频移性质对恢复结果进行亚像素位移相位补偿，得到相位恢复结果：

[0018]

[0019]

[0020] 其中F[…]是傅里叶变换，F-1[…]是傅里叶逆变换，t为时间，w为频率，w0为频移量，f(t)为时域函数，F(w)为频域函数。

[0021] 本发明与传统方法相比，具有如下优点：（1）成像精度上得到提高，传统的方法对离散的二维矩阵进行分析得到的是频谱的最大值对应的整数坐标，这样是不够准确的。本发明提出的方法将离散问题连续化得到最大值对应的坐标更为精确。（2）成像质量上得到提高，本发明对频谱图引入的二维高斯函数滤波，对频谱的部分高频成分实现了滤除。然而图像中的噪声部分往往体现在频谱上的高频部分，因此提升成像的质量。（3）不影响实时性，本发明虽然与传统方式相比在算法上变得更加丰富，但是并没有影响数字全息的实时性过程。

附图说明

[0022] 图1为基于二维曲面拟合的数字全息频谱中心亚像素搜索方法流程图。

[0023] 图2为基于二维曲面拟合的数字全息频谱中心亚像素搜索方法系统原理图。

[0024] 图3为本发明实施例中数字全息系统干涉记录光路图。

[0025] 图4为采用二维高斯函数频谱滤波去噪效果示意图。

[0026] 图5以二维Hann窗函数为目标函数的最小二乘法曲面拟合原理图。

[0027] 图6本发明实施例的+1频谱亚像素中心搜索与传统的+1谱中心搜索在数字全息相位恢复上的效果对比图。

[0028] 图3中，1-激光器，2-空间滤波器，3-分光棱镜，4-镜子，5-样品，6-相机。

具体实施方式

[0029] 为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本申请，并不用于限定本申请。

[0030] 本发明基于最小二乘曲面拟合实现数字全息频谱中心亚像素搜索方法需要在数字全息系统下完成，搭建的全息系统利用相干性很强的激光作为记录光源，通过马赫-曾德尔干涉法对具有一定透过率的微观物体进行干涉记录采集到条纹干涉图IH(x,y)，其中H为全息Hologram的缩写，x,y分别代表采集到的干涉图的横纵坐标。具体的干涉记录光路如图3所示，其中激光器1发射激光经空间滤波器2调制后通过分光棱镜3分为物光与参考光，物光照射样品5后与参考光发生干涉得到条纹图通过相机6采集。

[0031] 得到干涉条纹图后通过数字的方法进行相位恢复，本发明主要针对+1级频谱的数字处理过程进行优化，具体的过程如下。

[0032] （1）频谱图采集：利用马赫-曾德尔干涉法采集待测样品的干涉条纹图，通过快速傅里叶变换得到对应的频谱图。

[0033] （2）+1级频谱筛选：根据离轴干涉强度分布公式：

[0034]

[0035] 其中O和R分别为物光和参考光，“*”为共轭算子，k（k=2π/λ，λ为光波长）为光波矢，i为虚数单位，θ为离轴全息夹角，x、y分别代表干涉强度分布的横纵坐标。可以发现频谱分布为零级 、+1级 、-1级 。由于+1谱的物光分布为原始物光分布并非其共轭项，因此可以产生更好的数字再现结果也更适合观测。根据物镜的数值孔径和照明光波长计算得到+1级频谱带宽生成二维矩形窗函数实现对频谱的带通滤波，得到频谱图中载有样品信息的+1级谱。

[0036] （3）二维高斯函数滤波：将得到的+1级频谱利用二维高斯窗函数进行低通滤波，滤掉部分高频成分来降低噪声的影响。

[0037] 高斯函数在统计学中是正态分布的密度函数，函数的形态特点是关于数据期望值μ对称宽度受方差σ2影响的钟形对称函数，为方便对离散的二维矩阵进行滤波操作需要将连续高斯函数离散化：

[0038]

[0039] 其中p、q为像素坐标，P、Q为像素区间，σ2为方差。将二维高斯函数离散化得到二维离散高斯函数，并将其与截取到的+1级频谱相乘实现在频谱域的噪声滤除。具体的滤除效果由图4二维高斯函数频谱滤波去噪效果示意图所示，其中频谱图的高频成分得到了抑制，从a与a’，b与b’的区域对比中可以看出高频的噪声成分得到了滤除。图4中左下框内的上方小框为a区域、下方小框为b区域，图4中右下框内的上方小框为a’ 区域、下方小框为b’区域。

[0040] （4）二维曲面拟合得到亚像素的频谱中心：通过曲面拟合方法，将经过低通滤波后的+1级离散二维频谱矩阵拟合为连续的二维曲面函数。根据得到的曲面函数利用微分方法求得频谱区域极大值点位置，该位置即为精确到亚像素的频谱中心。

[0041] 经过降噪后的+1级频谱，由于其离散性，采用传统的中心寻址方法只能得到像素为整数的谱中心位置缺乏精确性，因此利用最小二乘方法以二维Hann窗函数为拟合目标函数进行曲面拟合。这一过程的难点在于拟合目标的选取，由于经过位移变换的+1级频谱的数据形态是中心值高四周低的，因此选用函数形态为中间高四周低的二维Hann窗函数。

[0042] Hann窗又称升余弦窗，Hann窗可以看作是3个矩形时间窗函数的频谱之和，或者说是 3个sinc(t)型函数之和。具体的拟合对象函数表达为：

[0043]

[0044] 其中a为方程的待定系数，N为拟合区间，n为像素位置，n=1,2,3，…，N。

[0045] 与二维高斯函数低通滤波去噪相似，将离散的Hann窗函数通过矩阵性质变换为二维Hann窗函数。针对选取的二维Hann窗函数为拟合目标，通过最小二乘法来最小化误差的平方和寻找数据的最佳函数匹配，利用最小二乘法可以简便地求得未知数据，并使这些求得的数据与实际数据之间误差的平方和为最小。算法实现原理，首先定义观测结果为y、预测结果为y’，两者之间的误差的大小为Rss，建立一个损失函数来衡量估计值和实际之间的误差的大小，并将最小化损失函数作为一个约束条件来求出参数向量的最优解。

[0046]

[0047] 按照前面的定义，Rss是y与y’ 间的差，那么当Rss无限趋近于0的时候，则y≈y’，即所求得的预测结果就等于实际结果。根据最小二乘方法得到拟合目标的待定系数，带入原方程即可得到拟合的二维曲面方程。

[0048] 由于+1级频谱的数据形态是中间数据高四周低，因此将得到的二维曲面方程利用微分的方法求得曲面方程的各点切面，找到水平切面对应的坐标值即为+1级的频谱中心坐标，不同于传统方法得到的整数坐标，本发明的方法得到的中心坐标精度可以达到亚像素级。图5展示了以Hann窗函数为目标函数的最小二乘法曲面拟合原理，先是构建二维Hann窗函数为曲面拟合的目标函数，然后利用最小二乘法对+1级频谱连续化拟合。

[0049] （5）相位恢复及亚像素位移相位校准：将得到的频谱中心位置取整后置于空矩阵中心并逆傅里叶变换到空间域，结合谱中心的亚像素偏差，利用傅里叶变换性质对恢复结果进行亚像素位移相位补偿，得到更为准确的相位恢复结果。

[0050] 通过亚像素位移法进行频谱中心移动实现相位补偿。然而传统的方法是直接将整数的中心像素放置于空矩阵中心位置进行相位恢复，由于本发明得到的+1级谱中心是精确到亚像素的，因此通过常规的方法是无法将中心置于矩阵中心的。需要采用亚像素位移的手段实现相位补偿来得到更好的相位恢复结果。具体的过程为，首先将亚像素的频谱中心进行向下取整处理 、 ，然后将整数化后的频谱中心置于空矩阵中心位置，因为整数化的缘故产生了亚像素级的偏差，所
以需要对+1级频谱进行亚像素位移相位补偿。具体的原理要根据傅里叶变换时移/频移的性质：

[0051]

[0052]

[0053] 其中F[…]是傅里叶变换，F-1[…]是傅里叶逆变换，t为时间，w为频率，w0为频移量，f(t)为时域函数，F(w)为频域函数。频谱域的位移量w0对应在时间域的 项，因此要实现频率域的亚像素位移，首先则需要将置于空矩阵的+1级频谱逆傅里叶变换到时间域，再根据要进行移动的频移量w0在对应的时间域上乘以 项以实现亚像素位移相位补偿。

[0054] 补偿后的结果与原有结果的对比图如图6所示，左侧框内的a区域，为传统的+1谱中心搜索在数字全息相位恢复上的效果图，右侧框内的a’ 域，为本实施例的+1频谱亚像素中心搜索与传统的+1谱中心搜索在数字全息相位恢复上的效果图。从a与a’ 区域对比上可以发现，a区域颜色偏白存在相位倾斜，a’区域对相位倾斜进行了修正。

[0055] 完整的算法实现流程图的文字描述如图1所示，为方便展示采用基于最小二乘曲面拟合实现数字全息频谱中心亚像素搜索的算法效果。图2提供了系统算法原理图，该图对算法的每一步结果都加以图像效果说明，更为直观和具体的展示该发明对于数字全息相位恢复中+1级频谱中心寻址的优化效果。

[0056] 与传统的简单粗暴的寻找频谱最大值对应的整数像素坐标不同的是，本实施例不仅采用二维高斯函数进行频域去噪，同时以二维Hann窗函数为拟合目标对+1级频谱进行最小二乘法曲面拟合，将离散的频谱矩阵变为连续的曲面函数进而得到精确到亚像素的频谱中心。与传统方法相比，本实施例得到的频谱中心坐标精度可以达到亚像素因此具有更好的成像效果和测量精度。

[0057] 对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本申请。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。