基于降采样相位畸变主成分分析的快速相位畸变补偿方法
【技术领域】
[0001] 本发明属于光学成像校正技术，特别是一种基于降米样相位畸变主成分分析的快 速相位畸变补偿方法。
【背景技术】
[0002] 数字全息作为一种新型相干测量与成像技术，其最大的优势在于能够同时独立地 获取物体的定量振幅信息和相位信息。当对生物细胞等几乎透明的物体进行定量检测时， 相位信息显得尤为重要。然而，要准确地获得物体的相位像，就必须要首先对再现像中的相 位畸变进行校正，这就需要知道实验中的各个参数，如记录距离、物参夹角、显微物镜的放 大倍率等。同时，相位畸变校正对记录参数的精度要求非常高。但是，实际的测量值通常与 理论值有偏差，这就会导致再现像的相位分布发生相位畸变，从而无法准确获得物体的相 位信息。因此，近几年相位畸变校正已经成为国内外研究者的关注热点。
[0003] 按实现方法分类，相位畸变校正可分为两类：一类是通过硬件方法即在实验记录 中通过设计相应的系统光路来消除相位畸变，比如较为典型的意大利研究组P. Ferraro 等人提出一种有效的两步曝光法（[1]P. Ferraro，et al. "Compensation ofthe Inherent Wave Front Curvature in Digital Holographic Coherent Microscopy for Quantitative Phase-contrast Imaging".Appl. Opt. 2003,42(11) :1938 ~1946)，这种方 法的思路是分别拍摄有样品和无样品时的两幅全息图再进行相位相减，可以一次去掉所有 畸变。然而，两步曝光法需记录两幅全息图，对系统稳定性要求较高。美国的M.K. Kim研 究组提出一种物理补偿的方法（[2]Μ· K. Kim. "Applications of Digital Holography in Biomedical Microscopy". J. Opt. Soc. Korea. 2010,14(2) :77 ~89)，在参考光路适当位置 放置一个与物光路相同的显微物镜，即使球面参考光的曲率与物光的曲率相同，这种方法 可以在直接在实验记录中消除二次相位畸变，但是，这种方法要求在实验中精确调整参考 光路中的显微物镜位置，这在实验中是很难实现的。另一类是通过软件方法，即在计算机 中，通过数值再现进行畸变校正，瑞士研究组T. Colomb等人提出了三种消除相位畸变的方 法（[3]T. Colomb，et al. "Numerical Parametric Lens for Shifting, Magnification, and Complete Aberration Compensation in Digital Holographic Microscopy'，· J. Opt. Soc. Am. A. 2006,23(12) :3177~3190)，第一种是自动相位掩膜法，通过多次曲线拟合法 来自动确定重建参数来校正相位畸变，第二种是提出利用参考共辄全息图来校正相位畸 变，第三种是选取再现视场中无物体的平坦区域，通过泽尼克多项式拟合来获得畸变相 位。在国内，西北工业大学的赵建林等人提出了基于最小二乘曲面拟合法（[4]J.L.Di，et al. "Phase Aberration Compensation of Digital Holographic Microscopy based on Least Squares Surface Fitting". Opt. Commun. · 2009,（282) :3873 ~3877)，只需一幅全 息图，便消除相位畸变。但是上述通过后期计算进行相位畸变补偿的方法大多运算量大，不 论是泽尼克多项式拟合还是最小二乘曲面拟合等拟合方法的计算时间都随拍摄得到的全 息图尺寸大小的增大而急剧增大，如何实现精度又高速度又快的相位畸变补偿成为了数字 全息显微成像中一项技术难题。
【发明内容】
[0004] 本发明的目的在于提供一种基于降采样相位畸变主成分分析的快速相位畸变补 偿方法，实现了快速的数字全息显微成像中相位畸变补偿。
[0005] 实现本发明目的的技术解决方案为：一种基于降采样相位畸变主成分分析的快速 相位畸变补偿方法，步骤如下：
[0006] 步骤一，利用相机拍摄得到全息图，然后利用傅立叶变换得到全息图的频谱，并提 取出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央；
[0007] 步骤二，从频谱中提取出相位畸变成分频谱；
[0008] 步骤三，逆傅立叶变换并提取出物体降采样的相位图；
[0009] 步骤四，奇异值分解提取出相位图的第一主成分；
[0010] 步骤五，最小二乘法拟合得出X，y方向上的抛物线方程；
[0011] 步骤六，排除错误点后再次拟合得出X，y方向上的准确的抛物线方程；
[0012] 步骤七，对X，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图；
[0013] 步骤八，物体原始相位图减去畸变相位图，即得到补偿后的相位图。
[0014] 本发明与现有技术相比，其显著优点：（1)本发明只对小尺寸的畸变频谱进行主 成分分析，而没有对大尺寸的整幅频谱进行主成分分析，这样显著降低了主成分分析的计 算时间，提高了相位畸变补偿速度。（2)本发明对降采样带来的错误点进行了排除，又对拟 合得到的相位畸变进行了亚像素插值，这样既能准确拟合出相位畸变函数，又不会出现由 频谱截断带来的振铃效应，保证了相位畸变补偿精度。
[0015] 下面结合附图对本发明作进一步详细描述。
【附图说明】
[0016] 图1为本发明基于降采样相位畸变主成分分析的快速相位畸变补偿方法的流程 图。
[0017] 图2为测试对象为对人体巨噬细胞进行数字全息显微成像时的快速相位畸变补 偿结果；图2(a)是数字全息显微镜拍摄到的原始干涉图；图2(b)是原始干涉图2(a)经过 傅立叶变换的频谱，图中用小框框出的是+1级谱；图2(c)是+1级谱平移到频谱中央后的 结果，即物体的原始频谱；图2(d)是+1级谱的放大图像，图中用小框框出的是畸变频谱； 图2(e)是降采样的相位畸变的第一主成分图；图2(f)是相位畸变的第一主成分列向量进 行抛物线拟合的曲线图；图2(g)是排除错误点后进行抛物线拟合的曲线图；图2(h)是升 采样的相位畸变图；图2 (i)是包含相位畸变的物体的相位图；图2 (j)是相位畸变补偿以 后的物体的相位图。
[0018] 图3为不进行升采样的相位补偿方法得到的相位畸变图。
【具体实施方式】
[0019] 结合图1，本发明基于降采样相位畸变主成分分析的快速相位畸变补偿方法，首先 利用数字全息显微系统拍摄得到全息图。然后利用傅立叶变换得到全息图的频谱。再提取 出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央。之后在+1级谱中提取出相位畸 变主成分频谱。并利用逆傅立叶变换并提取出物体降采样的相位图。然后利用奇异值分解 提取出相位图的第一主成分。再用最小二乘法拟合得出X，y方向上的抛物线方程。接下来 排除错误点后再次拟合得出X，y方向上的准确的抛物线方程。之后对X，y方向上的抛物线 进行亚像素插值，得到升采样的畸变相位图。最后将物体原始相位图减去畸变相位图，即得 到补偿后的相位图，具体实现步骤如下：
[0020] 步骤一：利用相机拍摄得到全息图，然后利用傅立叶变换得到全息图的频谱，并提 取出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央。
[0021] 步骤二：从+1级频谱中提取出相位畸变成分频谱。从频谱中提取出相位畸变成分 频谱，可以通过手动的从频谱中框选出绝对值较大的部分作为相位畸变成分频谱，也可以 利用如下公式进行自动选取：
[0022]
[0023]
[0024]
[0025]
[0026]
[0027]
[0028]
[0029]
[0030] 其中，F为经过步骤一至三处理后的物体原始频谱，abs (·)为求绝对值的函数， Fabs表示物体原始频谱F的绝对值，max( ·)为求最大值的函数，T是人为设定的阈值，本发 明中将阈值设定为物体原始频谱绝对值Fabs中最大值的千分之一。find(F abs> T)表示找出 Fabs中绝对值大于阈值T的所有点，并得到它们的横坐标和纵坐标，row，col分别表示求得 的横坐标和纵坐标。min( ·)为求最小值的函数，rowmax，rowmin，colmax，Col min分别表示求得 的横坐标的最大值，横坐标的最小值，纵坐标的最大值和纵坐标的最小值。F (rowmin:row _， colmin:col _)表示从F中提取出的第rowmin行到第rowmax行、第col min列到第col _列的频 谱成分，Fa表示提取出的相位畸变成分频谱。
[0031] 步骤三：逆傅立叶变换并提取出物体降采样的相位图。逆傅立叶变换并提取出物 体降采样的相位图的公式为：
[0032]
[0033]
[0034] 其中，ifft2 (·)为求二维逆傅立叶变换的函数，WA表示相位畸变成分的复振幅， angle( ·)为复振幅中相位的函数，Pa为相位畸变成分。
[0035] 步骤四：奇异值分解提取出相位图的第一主成分。奇异值分解提取出相位图的第 一主成分的公式为：
[0036] [SV, D] = svd (Pa)
[0037] Pai= S(1，1)V( :，1)D(1，：）
[0038] 其中，svd( ·)为求矩阵的奇异值分解函数，S表示奇异值分解所得的主值矩阵，V 表不奇异值分解所得的列主成分矩阵，D表不奇异值分解所得的行主成分矩阵。S(l，1)表 不主值矩阵S中第一行第一列的数值，即第一主成分的主值，V( ;，1)表不列主成分矩阵V 的第一列，即第一主成分的列向量，D(l，：）表不行主成分矩阵D的第一行，即第一主成分的 行向量，Piu表不提取出的相位图Pa的第一主成分，它是由第一主成分的主值S (1，1)乘以第 一主成分的列向量V( :，1)再乘以第一主成分的行向量D(1，：）所得。
[0039] 步骤五：最小二乘法拟合得出x(横向）和y (纵向）两个方向上的抛物线方程。 最小二乘法拟合得出X和y两个方向上的抛物线方程的公式为：
[0040] [a2, a!，a。] = polyfit[V( :，1)，2]<