CN105157561B - 基于降采样相位畸变主成分分析的快速相位畸变补偿方法 - Google Patents

Abstract

本发明公开了一种基于降采样相位畸变主成分分析的快速相位畸变补偿方法，首先利用相机拍摄得到全息图，得到全息图的频谱后提取出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央；从频谱中提取出相位畸变成分频谱，逆傅立叶变换并提取出物体降采样的相位图；奇异值分解提取出相位图的第一主成分，最小二乘法拟合得出x，y方向上的抛物线方程；排除错误点后再次拟合得出x，y方向上的准确的抛物线方程；对x，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图；物体原始相位图减去畸变相位图，即得到补偿后的相位图。本发明既能降低主成分分析的计算时间，提高了相位畸变补偿速度，又能准确拟合出相位畸变函数。

Description

基于降采样相位畸变主成分分析的快速相位畸变补偿方法

技术领域

本发明属于光学成像校正技术，特别是一种基于降采样相位畸变主成分分析的快速相位畸变补偿方法。

背景技术

数字全息作为一种新型相干测量与成像技术，其最大的优势在于能够同时独立地获取物体的定量振幅信息和相位信息。当对生物细胞等几乎透明的物体进行定量检测时，相位信息显得尤为重要。然而，要准确地获得物体的相位像，就必须要首先对再现像中的相位畸变进行校正，这就需要知道实验中的各个参数，如记录距离、物参夹角、显微物镜的放大倍率等。同时，相位畸变校正对记录参数的精度要求非常高。但是，实际的测量值通常与理论值有偏差，这就会导致再现像的相位分布发生相位畸变，从而无法准确获得物体的相位信息。因此，近几年相位畸变校正已经成为国内外研究者的关注热点。

按实现方法分类，相位畸变校正可分为两类：一类是通过硬件方法即在实验记录中通过设计相应的系统光路来消除相位畸变，比如较为典型的意大利研究组P.Ferraro等人提出一种有效的两步曝光法([1]P.Ferraro，et al.“Compensation ofthe InherentWave Front Curvature in Digital Holographic Coherent Microscopy forQuantitative Phase-contrast Imaging”.Appl.Opt.2003，42(11)：1938～1946)，这种方法的思路是分别拍摄有样品和无样品时的两幅全息图再进行相位相减，可以一次去掉所有畸变。然而，两步曝光法需记录两幅全息图，对系统稳定性要求较高。美国的M.K.Kim研究组提出一种物理补偿的方法([2]M.K.Kim.“Applications of Digital Holography inBiomedical Microscopy”.J.Opt.Soc.Korea.2010，14(2)：77～89)，在参考光路适当位置放置一个与物光路相同的显微物镜，即使球面参考光的曲率与物光的曲率相同，这种方法可以在直接在实验记录中消除二次相位畸变，但是，这种方法要求在实验中精确调整参考光路中的显微物镜位置，这在实验中是很难实现的。另一类是通过软件方法，即在计算机中，通过数值再现进行畸变校正，瑞士研究组T.Colomb等人提出了三种消除相位畸变的方法([3]T.Colomb，et a1.“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，etal.“Phase Aberration Compensation of Digital Holographic Microscopy based onLeast Squares Surface Fitting”.Opt.Commun..2009，(282)：3873～3877)，只需一幅全息图，便消除相位畸变。但是上述通过后期计算进行相位畸变补偿的方法大多运算量大，不论是泽尼克多项式拟合还是最小二乘曲面拟合等拟合方法的计算时间都随拍摄得到的全息图尺寸大小的增大而急剧增大，如何实现精度又高速度又快的相位畸变补偿成为了数字全息显微成像中一项技术难题。

发明内容

本发明的目的在于提供一种基于降采样相位畸变主成分分析的快速相位畸变补偿方法，实现了快速的数字全息显微成像中相位畸变补偿。

实现本发明目的的技术解决方案为：一种基于降采样相位畸变主成分分析的快速相位畸变补偿方法，步骤如下：

步骤一，利用相机拍摄得到全息图，然后利用傅立叶变换得到全息图的频谱，并提取出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央；

步骤二，从频谱中提取出相位畸变成分频谱；

步骤三，逆傅立叶变换并提取出物体降采样的相位图；

步骤四，奇异值分解提取出相位图的第一主成分；

步骤五，最小二乘法拟合得出x，y方向上的抛物线方程；

步骤六，排除错误点后再次拟合得出x，y方向上的准确的抛物线方程；

步骤七，对x，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图；

步骤八，物体原始相位图减去畸变相位图，即得到补偿后的相位图。

本发明与现有技术相比，其显著优点：(1)本发明只对小尺寸的畸变频谱进行主成分分析，而没有对大尺寸的整幅频谱进行主成分分析，这样显著降低了主成分分析的计算时间，提高了相位畸变补偿速度。(2)本发明对降采样带来的错误点进行了排除，又对拟合得到的相位畸变进行了亚像素插值，这样既能准确拟合出相位畸变函数，又不会出现由频谱截断带来的振铃效应，保证了相位畸变补偿精度。

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

附图说明

图1为本发明基于降采样相位畸变主成分分析的快速相位畸变补偿方法的流程图。

图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)是相位畸变补偿以后的物体的相位图。

图3为不进行升采样的相位补偿方法得到的相位畸变图。

具体实施方式

结合图1，本发明基于降采样相位畸变主成分分析的快速相位畸变补偿方法，首先利用数字全息显微系统拍摄得到全息图。然后利用傅立叶变换得到全息图的频谱。再提取出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央。之后在+1级谱中提取出相位畸变主成分频谱。并利用逆傅立叶变换并提取出物体降采样的相位图。然后利用奇异值分解提取出相位图的第一主成分。再用最小二乘法拟合得出x，y方向上的抛物线方程。接下来排除错误点后再次拟合得出x，y方向上的准确的抛物线方程。之后对x，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图。最后将物体原始相位图减去畸变相位图，即得到补偿后的相位图，具体实现步骤如下：

步骤一：利用相机拍摄得到全息图，然后利用傅立叶变换得到全息图的频谱，并提取出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央。

步骤二：从+1级频谱中提取出相位畸变成分频谱。从频谱中提取出相位畸变成分频谱，可以通过手动的从频谱中框选出绝对值较大的部分作为相位畸变成分频谱，也可以利用如下公式进行自动选取：

F_abs＝abs(F)

[row，col]＝find(F_abs＞T)

row_max＝max(row)

row_min＝min(rpw)

col_max＝max(col)

col_min＝min(col)

F_A＝F(row_min：row_max，col_min：col_max)

其中，F为经过步骤一至三处理后的物体原始频谱，abs(·)为求绝对值的函数，F_abs表示物体原始频谱F的绝对值，max(·)为求最大值的函数，T是人为设定的阈值，本发明中将阈值设定为物体原始频谱绝对值F_abs中最大值的千分之一。find(F_abs＞T)表示找出F_abs中绝对值大于阈值T的所有点，并得到它们的横坐标和纵坐标，row，col分别表示求得的横坐标和纵坐标。min(·)为求最小值的函数，row_max，row_min，col_max，col_min分别表示求得的横坐标的最大值，横坐标的最小值，纵坐标的最大值和纵坐标的最小值。F(row_min：row_max，col_min：col_max)表示从F中提取出的第row_min行到第row_max行、第col_min列到第col_max列的频谱成分，F_A表示提取出的相位畸变成分频谱。

步骤三：逆傅立叶变换并提取出物体降采样的相位图。逆傅立叶变换并提取出物体降采样的相位图的公式为：

W_A＝ifft2(F_A)

P_A＝angle(W_A)

其中，ifft2(·)为求二维逆傅立叶变换的函数，WA表示相位畸变成分的复振幅，angle(·)为复振幅中相位的函数，P_A为相位畸变成分。

步骤四：奇异值分解提取出相位图的第一主成分。奇异值分解提取出相位图的第一主成分的公式为：

[SV，D]＝svd(P_A)

P_A1＝S(1，1)V(：，1)D(1，：)

其中，svd(·)为求矩阵的奇异值分解函数，S表示奇异值分解所得的主值矩阵，V表示奇异值分解所得的列主成分矩阵，D表示奇异值分解所得的行主成分矩阵。S(1，1)表示主值矩阵S中第一行第一列的数值，即第一主成分的主值，V(；，1)表示列主成分矩阵V的第一列，即第一主成分的列向量，D(1，：)表示行主成分矩阵D的第一行，即第一主成分的行向量，P_A1表示提取出的相位图P_A的第一主成分，它是由第一主成分的主值S(1，1)乘以第一主成分的列向量V(：，1)再乘以第一主成分的行向量D(1，：)所得。

步骤五：最小二乘法拟合得出x(横向)和y(纵向)两个方向上的抛物线方程。最小二乘法拟合得出x和y两个方向上的抛物线方程的公式为：

[a₂，a₁，a₀]＝polyfit[V(：，1)，2]

[b₂，b₁，b₀]＝polyfit[D(1，：)，2]

X₁＝1：1：hig_A

Y₁＝1：1：wid_A

其中，polyfit(·，2)为最小二乘法拟合抛物线方程的函数，a₂，a₁，a₀分别表示V(：，1)拟合的抛物线方程的二次项、一次项、常数项的系数，b₂，b₁，b₀分别表示D(1，：)拟合的抛物线方程的二次项、一次项、常数项的系数。hig_A为P_A的行数，wid_A为P_A的列数，X₁表示从1到hig_A间隔为1的列向量，Y₁表示从1到wid_A间隔为1的行向量，V₁表示拟合出的列向量抛物线方程，D₁表示拟合出的行向量抛物线方程。

步骤六：排除错误点后再次拟合得出x，y方向上的准确的抛物线方程，步骤为：

第一步，将V(：，1)与V₁相减并取绝对值，找出其中大于阈值T_VD的点记为错误点，然后排除V(：，1)中的错误点，将剩余向量记为V₂。阈值T_VD的选取范围一般在之间，本发明中取阈值

第二步，将D(：，1)与D₁相减并取绝对值，找出其中大于阈值T_VD的点记为错误点，然后排除D(：，1)中的错误点，将剩余向量记为D₂。

第三步，按照如下公式拟合出准确的x，y方向上的抛物线方程

[c₂，c₁，c₀]＝polyfit[V₂，2]

[d₂，d₁，d₀]＝polyfit[D₂，2]

其中，c₂，c₁，c₀分别表示V₂拟合的抛物线方程的二次项、一次项、常数项的系数，d₂，d₁，d₀分别表示D₂拟合的抛物线方程的二次项、一次项、常数项的系数。

步骤七：对x，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图。对x，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图的公式为：

P_A2＝S(1，1)V₃D₃

其中，hig为F的行数，wid为F的列数，X₂表示从1到hig_A间隔为的亚像素列向量，Y₂表示从1到wid_A间隔为的亚像素行向量，V₃表示升采样的列向量抛物线方程，D₃表示升采样的行向量抛物线方程，P_A2表示升采样的畸变相位图。

步骤八：物体原始相位图减去畸变相位图，即得到补偿后的相位图。

为了测试基于降采样相位畸变主成分分析的快速相位畸变补偿方法的有效性，选取人体巨噬细胞进行数字全息显微成像时的快速相位畸变补偿。如图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)是相位畸变补偿以后的物体的相位图。对比图2(i)和图2(j)可以看出物体原始的相位图中含有大量的相位畸变，严重影响了对物体的识别和理解，而相位畸变补偿后的相位图背景均匀，细胞的相位信息突出而且清晰，说明本方法能够有效的补偿相位畸变。对比其他相位补偿算法，例如不进行生采样的相位补偿方法得到的畸变相位图，如图3所示，可以发现不进行生采样的相位补偿方法会造成图像右侧严重的频谱截断效应，导致相位补偿出现误差。此外，本方法进行主成分分析和抛物线拟合的计算时间总共仅仅需要不到0.2s，而其他整幅图像进行主成分分析或者曲面拟合的方法所需的计算时间都超过20s，可见本发明实现了快速相位畸变补偿，补偿速度是其他方法的100倍以上。

Claims (5)

1.一种基于降采样相位畸变主成分分析的快速相位畸变补偿方法，其特征在于步骤如下：

步骤一，利用相机拍摄得到全息图，然后利用傅立叶变换得到全息图的频谱，并提取出+1级谱，去除其余频谱，并将+1级谱平移到频谱中央；

步骤二，从频谱中提取出相位畸变成分频谱；

步骤三，逆傅立叶变换并提取出物体降采样的相位图；

步骤四，奇异值分解提取出相位图的第一主成分；

步骤五，最小二乘法拟合得出x，y方向上的抛物线方程；

步骤六，排除错误点后再次拟合得出x，y方向上的准确的抛物线方程；

步骤七，对x，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图；

步骤八，物体原始相位图减去畸变相位图，即得到补偿后的相位图；

在步骤六中，排除错误点后再次拟合得出x，y方向上的准确的抛物线方程的步骤为：

第一步，将V(:,1)与V₁相减并取绝对值，找出其中大于阈值T_VD的点记为错误点，然后排除V(:,1)中的错误点，将剩余向量记为V₂，阈值T_VD的选取范围在之间；

第二步，将D(:,1)与D₁相减并取绝对值，找出其中大于阈值T_VD的点记为错误点，然后排除D(:,1)中的错误点，将剩余向量记为D₂；

第三步，按照如下公式拟合出准确的x，y方向上的抛物线方程

[c₂,c₁,c₀]＝polyfit[V₂,2]

[d₂,d₁,d₀]＝polyfit[D₂,2]

其中，c₂,c₁,c₀分别表示V₂拟合的抛物线方程的二次项、一次项、常数项的系数，d₂,d₁,d₀分别表示D₂拟合的抛物线方程的二次项、一次项、常数项的系数；

在步骤七中，对x，y方向上的抛物线进行亚像素插值，得到升采样的畸变相位图的公式为：

P_A2＝S(1,1)V₃D₃

其中，hig为F的行数，wid为F的列数，X₂表示从1到hig_A间隔为的亚像素列向量，Y₂表示从1到wid_A间隔为的亚像素行向量，V₃表示升采样的列向量抛物线方程，D₃表示升采样的行向量抛物线方程,P_A2表示升采样的畸变相位图。

2.根据权利要求1所述的基于降采样相位畸变主成分分析的快速相位畸变补偿方法，其特征在于步骤二中，从频谱中提取出相位畸变成分频谱，通过手动的从频谱中框选出绝对值较高的部分作为相位畸变成分频谱，或者利用如下公式进行自动选取：

F_abs＝abs(F)

[row,col]＝find(F_abs>T)

row_max＝max(row)

row_min＝min(row)

col_max＝max(col)

col_min＝min(col)

F_A＝F(row_min:row_max,col_min:col_max)

其中，F为经过处理后的物体原始频谱，abs(·)为求绝对值的函数，F_abs表示物体原始频谱F的绝对值，max(·)为求最大值的函数，T是设定的阈值，find(F_abs>T)表示找出F_abs中绝对值大于阈值T的所有点，并得到它们的横坐标和纵坐标，row,col分别表示求得的横坐标和纵坐标；min(·)为求最小值的函数，row_max,row_min,col_max,col_min分别表示求得的横坐标的最大值，横坐标的最小值，纵坐标的最大值和纵坐标的最小值；F(row_min:row_max,col_min:col_max)表示从F中提取出的第row_min行到第row_max行、第col_min列到第col_max列的频谱成分，F_A表示提取出的相位畸变成分频谱。

3.根据权利要求1所述的基于降采样相位畸变主成分分析的快速相位畸变补偿方法，其特征在于步骤三中，逆傅立叶变换并提取出物体降采样的相位图的公式为：

W_A＝ifft2(F_A)

P_A＝angle(W_A)

其中，ifft2(·)为求二维逆傅立叶变换的函数，W_A表示相位畸变成分的复振幅，angle(·)为复振幅中相位的函数，P_A为相位畸变成分。

4.根据权利要求1所述的基于降采样相位畸变主成分分析的快速相位畸变补偿方法，其特征在于步骤四中，奇异值分解提取出相位图的第一主成分的公式为：

[S,V,D]＝svd(P_A)

P_A1＝S(1,1)V(:,1)D(1,:)

其中，svd(·)为求矩阵的奇异值分解函数，S表示奇异值分解所得的主值矩阵，V表示奇异值分解所得的列主成分矩阵，D表示奇异值分解所得的行主成分矩阵，S(1,1)表示主值矩阵S中第一行第一列的数值，即第一主成分的主值，V(:,1)表示列主成分矩阵V的第一列，即第一主成分的列向量，D(1,:)表示行主成分矩阵D的第一行，即第一主成分的行向量，P_A1表示提取出的相位图P_A的第一主成分，它是由第一主成分的主值S(1,1)乘以第一主成分的列向量V(:,1)再乘以第一主成分的行向量D(1,:)所得。

5.根据权利要求1所述的基于降采样相位畸变主成分分析的快速相位畸变补偿方法，其特征在于步骤五中，最小二乘法拟合得出x，y方向上的抛物线方程的公式为：

[a₂,a₁,a₀]＝polyfit[V(:,1),2]

[b₂,b₁,b₀]＝polyfit[D(1,:),2]

X₁＝1:1:hig_A

Y₁＝1:1:wid_A

其中，polyfit(·,2)为最小二乘法拟合抛物线方程的函数，a₂,a₁,a₀分别表示V(:,1)拟合的抛物线方程的二次项、一次项、常数项的系数，b₂,b₁,b₀分别表示D(1,:)拟合的抛物线方程的二次项、一次项、常数项的系数, hig_A为P_A的行数，wid_A为P_A的列数，X₁表示从1到hig_A间隔为1的列向量，Y₁表示从1到wid_A间隔为1的行向量，V₁表示拟合出的列向量抛物线方程，D₁表示拟合出的行向量抛物线方程。