CN103778650A - 主成分分析相位恢复算法 - Google Patents

Abstract

本发明公开了一种主成分分析相位恢复算法，包括重组干涉图、获得背景分量m_x、计算协方差矩阵C、计算对角化协方差矩阵U、得到主成分分量y、通过反正切函数求解相位6个步骤。该算法使用多张干涉图恢复出生物样品的定量相位，并且是一种时域处理方法，仅需要使用矩阵运算便可以得到定量相位分布，并且该算法需要的计算时间少于多数频率域算法，能够适应大数据量图像处理等需求。

Description

主成分分析相位恢复算法

技术领域

本发明涉及一种算法，特别是涉及一种主分成分析相位恢复算法。 

背景技术

基于正则化光学流场算法的定量相位流式细胞仪，使用的正则化光学流场算法的特点是仅需要两张干涉图便可以得到生物样品的定量相位分布，该算法提供了一种准实时测量观察生物细胞动态的方法和手段。然而该方法还是有一定的缺点。首先，该方法只使用两张移相干涉图，从其中的信息中可以获得定量的相位信息，然而，这些信息量不足以对噪声等外部影响起到足够的抑制作用。其次，使用三维精密微移平台操纵细胞样品在视场中进行运动时，其在视场中运动的过程中，使用CCD摄像机能够得到十数张干涉图。由于使用三维精密微移平台，导致细胞并不会有大范围的运动，因此，很多图像因此冗余，对最后的数据处理和分析体现的作用不大。最后，该方法还是必须使用频域处理的方法才能够得到最后的定量相位分布，频率域的处理往往会带来频率泄露，振铃效应等不良影响。 

发明内容

本发明所要解决的技术问题是，克服现有技术的缺点，提供一种属于时域处理，且仅需要使用矩阵运算便可以得到定量相位分布的算法。为了解决以上技术问题，本发明提供一种主分成分析相位恢复算法，包括以下步骤： 

步骤1重组干涉图：将每一张移相干涉图重组成行向量，然后将这N张移相干涉图的行向量按照列组合组成如下的x矩阵：x＝[x₁，x₂，...，x_N]^T，其中每一行都是每一张移相干涉图重组得到的一维数据，其长度为M＝N_x×N_y，N_x和N_y分别为移相干涉图x轴与y轴所占的像素点数，T代表为矩阵的转置，因此矩阵x为一N×M的矩阵； 

步骤2获得背景分量m_x：矩阵m_x具有和矩阵x同样的维度，而矩阵m_x中所有元素的值均一致，其是矩阵x中所有元素的平均值； 

步骤3计算协方差矩阵C：从矩阵x得到协方差矩阵C：C＝(x-m_x)(x-m_x)^T； 

步骤4计算对角化协方差矩阵U：由于协方差矩阵C是一实对称N×N矩阵，因此该协方差矩阵可以如下式实现对角化：D＝UCU^T，矩阵D为对角化协方差矩阵而矩阵U是一个正交变换的矩阵，其大小为都为N×N； 

步骤5得到主成分分量y：主成分便可以由矩阵D，矩阵x和矩阵m_x得到：y＝U(x-m_x)，其中，矩阵y的第一列和第二列便分别代表了主成分的正交特征值，其分别为I_c和I_s； 

步骤6通过反正切函数求解相位：Φ＝arctan(I_s/I_c)。 

本发明的有益效果是：该算法使用多张干涉图恢复出生物样品的定量相位，并且是一种时域处理方法，仅需要使用矩阵运算便可以得到定量相位分布，并且该算法需要的计算时间还少于多数频率域算法，能够适应大数据量图像处理等需求。 

具体实施方式

实施例1 

本实施例提供了一种主分成分析相位恢复算法，包括以下步骤： 

步骤1重组干涉图：将每一张移相干涉图重组成行向量，然后将这N张移相干涉图的行向量按照列组合组成如下的x矩阵：x＝[x₁，x₂，...，x_N]^T，其中每一行都是每一张移相干涉图重组得到的一维数据，其长度为M＝N_x×N_y，N_x和N_y分别为移相干涉图x轴与y轴所占的像素点数，T代表为矩阵的转置； 

步骤2获得背景分量m_x：矩阵m_x具有和矩阵x同样的维度，而矩阵m_x中所有元素的值均一致，其是矩阵x中所有元素的平均值； 

步骤3计算协方差矩阵C：从矩阵x得到协方差矩阵C：C＝(x-m_x)(x-m_x)^T； 

步骤4计算对角化协方差矩阵U：协方差矩阵C是一实对称N×N矩阵，因此该协方差矩阵可以如下式实现对角化：D＝UCU^T，矩阵D为对角化协方差矩阵而矩阵U是一个正交 变换的矩阵，其大小为都为N×N； 

步骤5得到主成分分量y：主成分便可以由矩阵D，矩阵x和矩阵m_x得到：y＝U(x-m_x)，其中，矩阵y的第一列和第二列便分别代表了主成分的正交特征值，其分别为I_c和I_s； 

步骤6通过反正切函数求解相位：Φ＝arctan(I_s/I_c)。 

随后，将相位重组成为二维图，然后对相位图进行解包消倾斜后，可以得到视场内的连续相位分布。 

主成分分析相位恢复算法的数学原理： 

为了说明主成分分析相位恢复算法的原理，以下使用数学推导的方式对主成分分析相位恢复算法对移相干涉图的处理进行分析。 

移相干涉图可以使用公式表达为： 

I_n(x，y)＝S(x，y)+B(x，y)cos(Φ(x，y)+δ_n)   (1.6) 

其中，I_n(x，y)，S(x，y)，B(x，y)分别是干涉图，背景光，干涉调制项的强度分布，而Φ(x，y)为待求解的相位分布，δ_n为移相相位。 

如果将以上干涉图重组成为列向量，公式(1.6)可以写为： 

I_n，k＝S_k+B_kcos(Φ_k+δ_n)   (1.7) 

式中，k＝[1，2，...，M]代表移相干涉图的像素点数，而n＝[1，2，...，N]代表不同移相角度的干涉图。参考公式(1.1)，矩阵x可以写为： 

$x=\left(\begin{array}{cccc}{S}_1+B_1\cos ({\mathrm{\Φ}}_1+{\mathrm{\δ}}_1)& S_2+B_2\cos ({\mathrm{\Φ}}_2+{\mathrm{\δ}}_1)& ...& S_M+B_M\cos ({\mathrm{\Φ}}_M+{\mathrm{\δ}}_1)\\ ...& ...& ...& ...\\ S_1+B_1\cos ({\mathrm{\Φ}}_1+{\mathrm{\δ}}_N)& S_2+B_2\cos ({\mathrm{\Φ}}_2+{\mathrm{\δ}}_N)& ...& S_M+B_M\cos ({\mathrm{\Φ}}_M+{\mathrm{\δ}}_N)\end{array}\right)---\left(1.8\right)$

每一张干涉图都为一M×1的行向量。 

假设背景光为一平滑分布，则矩阵m_x实际上代表了背景光的强度： 

$\{m_x{\}}_k=\frac{1}{N}\underset{n}{\mathrm{\Σ}}{I}_{n,k}\≈S_k---\left(1.1\right)$

设矩阵

则其可以表示为： 

$\tilde{x}=\left(\begin{array}{cccc}{B}_1\cos ({\mathrm{\Φ}}_1+{\mathrm{\δ}}_1)& B_2\cos ({\mathrm{\Φ}}_2+{\mathrm{\δ}}_1)& ...& B_M\cos ({\mathrm{\Φ}}_M+{\mathrm{\δ}}_1)\\ ...& ...& ...& ...\\ B_1\cos ({\mathrm{\Φ}}_1+{\mathrm{\δ}}_N)& B_2\cos ({\mathrm{\Φ}}_2+{\mathrm{\δ}}_N)& ...& B_M\cos ({\mathrm{\Φ}}_M+{\mathrm{\δ}}_N)\end{array}\right)---(1.10)$

由三角函数公式，矩阵

可以写为： 

$\tilde{x}=a_n{u}_k+b_n{v}_k---\left(1.11\right)$

其中： 

α_n＝cosδ_n   (1.11a) 

b_n＝-sinδ_n   (1.11b) 

u_k＝B_kcosΦ_k   (1.11c) 

v_k＝B_ksinΦ_k   (1.11d) 

那么协方差矩阵C可以写为： 

$C_{i,j}=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}(a_i{u}_k+b_i{v}_k)(a_j{u}_k+b_j{v}_k)---\left(1.12\right)$

将公式(1.12)展开，可以得到： 

$C_{i,j}=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{a}_i{a}_j{u}_k{u}_k+a_i{b}_j{u}_k{v}_k+a_j{b}_i{u}_k{v}_k+b_i{b}_j{v}_k{v}_k---\left(1.13\right)$

由于协方差矩阵C为实对称矩阵，其可以化简为： 

$C_{i,j}=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{a}_i{a}_j{u}_k{u}_k+2a_i{b}_j{u}_k{v}_k+b_i{b}_j{v}_k{v}_k---\left(1.14\right)$

在这里，令A_ij＝a_ia_j，E_ij＝2a_ib_j以及F_ij＝b_ib_j，上式可以化简为： 

$C_{i,j}=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{u}_k{u}_k+E_{\mathrm{ij}}{u}_k{v}_k+F_{\mathrm{ij}}{v}_k{v}_k---\left(1.15\right)$

如果在整个视场中有多条条纹，那么有以下近似： 

$\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{u}_k{v}_k=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{B}_k^2\cos {\mathrm{\&Phi;}}_k\sin {\mathrm{\&Phi;}}_k<<\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{u}_k{u}_k=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{B}_k^2\cos {\mathrm{\&Phi;}}_k^2---\left(1.16\right)$

$\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{u}_k{v}_k=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{B}_k^2\cos {\mathrm{\&Phi;}}_k\sin {\mathrm{\&Phi;}}_k<<\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{v}_k{v}_k=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{B}_k^2\sin {\mathrm{\&Phi;}}_k^2---\left(1.17\right)$

在该近似条件下，协方差矩阵C可以写为： 

C＝αA+βF   (1.18) 

其中矩阵A和矩阵F皆为N×N矩阵，并具有以下形式： 

A＝[cosδ₁，...，cosδ_N]^T[cosδ₁，...，cosδ_N]   (1.19) 

F＝[sinδ₁，...，sinδ_N]^T[sinδ₁，...，sinδ_N]   (1.20) 

而α和β可以使用下式表示： 

$\mathrm{\&alpha;}=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{u}_k{u}_k---\left(1.21\right)$

$\mathrm{\&beta;}=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{v}_k{v}_k---\left(1.22\right)$

因为矩阵A和矩阵F都来源于单一一维向量的乘积，因此它们的秩都为1，从而，它们只有单一的特征值以及特征向量。矩阵A和矩阵F的特征值为： 

${\mathrm{\&lambda;}}_A=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\cos }^2{\mathrm{\&delta;}}_i---\left(1.23\right)$

${\mathrm{\&lambda;}}_F=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\sin }^2{\mathrm{\&delta;}}_i---\left(1.24\right)$

而特征向量为： 

${\mathrm{\&omega;}}_A=\frac{[\cos {\mathrm{\&delta;}}_1,...,\cos {\mathrm{\&delta;}}_N]}{{\mathrm{\&lambda;}}_A}---\left(1.25\right)$

${\mathrm{\&omega;}}_F=\frac{[\sin {\mathrm{\&delta;}}_1,...,\sin {\mathrm{\&delta;}}_N]}{{\mathrm{\&lambda;}}_F}---\left(1.26\right)$

特征值与特征向量的检验方式可以通过Aω_A＝λ_Aω_A与Fω_F＝λ_Fω_F来进行检验。由公式(1.18)可知，协方差矩阵C的秩为2，并且有两个特征值并有与其对应的两个特征向量。如果相移的角度满足其至少相移了2π的范围，因此，下式可近似成立： 

$\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\cos {\mathrm{\&delta;}}_n\sin {\mathrm{\&delta;}}_n\&cong;0---\left(1.27\right)$

而且还满足： 

Aω_F＝0   (1.28) 

Fω_A＝0   (1.29) 

在这种情况下，协方差矩阵C的两个特征值分别为ω₁＝ω_A以及ω₂＝ω_F，其对应的特征向量为：λ₁＝αλ_A和λ₂＝αλ_F。 

基于协方差矩阵C的两个特征值以及特征向量，可以得到对角协方差矩阵D以及正交矩阵U。而对角协方差矩阵D的两个非零元素为：D₁₁＝λ₁以及D₂₂＝λ₂。此外，正交矩阵U的第一和第二行对应于协方差矩阵C的两个特征值ω₁和ω₂。在得到正交矩阵U之后，主成分便可以由以下两式给出： 

$y_{1,k}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{B}_k\cos ({\mathrm{\&Phi;}}_k+{\mathrm{\&delta;}}_n)\cos {\mathrm{\&delta;}}_n---\left(1.30\right)$

$y_{2,k}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{B}_k\cos ({\mathrm{\&Phi;}}_k+{\mathrm{\&delta;}}_n)\sin {\mathrm{\&delta;}}_n---\left(1.31\right)$

根据公式(1.27)给出的近似条件，可以得到： 

$y_{1,k}\&cong;B_k\cos {\mathrm{\&Phi;}}_k---\left(1.32\right)$

$y_{2,k}\&cong;-B_k\sin {\mathrm{\&Phi;}}_k---\left(1.33\right)$

最后，待求解的带有包裹的相位可以表示为： 

${\mathrm{\&Phi;}}_k=\arctan (-\frac{y_{2,k}}{y_{1,k}})---\left(1.34\right)$

公式(1.34)与使用传统的经典最小二乘法一致，然而该方法需要知道每一张干涉图的相移δ_n。在主成分分析中，相移可以通过协方差矩阵C的特征向量得到。 

以上介绍的方法包括两个步骤，首先，通过协方差矩阵得到相移估计值，之后，使用最小二乘解法得到相位值。该方法主要依赖于公式(1.9)，(1.16)，(1.17)和(1.27)的近似关系。在实验中，这些近似很难全部满足。在这种情况下，协方差矩阵的秩将不再是2，而将更大一些。这意味着，将有超过两个互不相关的主成分分量，其会携带相移变化的信息。 

主成分分析相位恢复算法用于定量相位流式细胞仪的原理和模拟 

当保存生物样品的样品室由三维精密微移平台在视场内进行扫描时，会得到一系列的离轴干涉图。然而，如果仅关注以某一生物样品为中心的视场，该范围的干涉条纹可以看作为移相干涉图。然后可以便可以使用主成分分析相位恢复算法用于定量相位恢复。该方法的主要特点就是无需知道移相的具体角度，非常适合在定量相位流式细胞仪中使用，并且其为时域算法，所用时间较短，并且避免了频域算法中的误差。 

以下，本节通过模拟的方式首先说明使用主成分分析相位恢复算法用于定量相位流式细胞仪的原理和模拟过程。 

在模拟中，使用血红细胞的模型作为生物样品，其形态可以使用如下公式进行表示： 

$z\left(\mathrm{\&rho;}\right)=[1-{\left(\frac{\mathrm{\&rho;}}{a}\right)}^2{]}^{1/2}[0.72+4.152{\left(\frac{\mathrm{\&rho;}}{a}\right)}^2-3.426{\left(\frac{\mathrm{\&rho;}}{a}\right)}^4]---\left(1.35\right)$

血红细胞为一凹圆盘状形貌，其直径为7.65微米，而且其内部构成较为均匀，可以看作为一相位物体。 

其相位与厚度的关系可以使用下式表示： 

其中L为血红细胞的厚度，λ为探测用波长，Δn为血红细胞折射率与折射率液(生理盐水)的折射率差。一般血红细胞的折射率在～1.4左右，而生理盐水的折射率约为～1.34。从而使用公式(1.36)可以得到血红细胞的相位分布。 

当血红细胞模型在视场内水平移动时，使用采集装置可以得到一系列的干涉图，如图1，所示。 

图1中，(a)-(c)血红细胞在视场内扫描得到一系列的离轴干涉图；(d)-(f)仅以某一生物样品为中心的视场，该范围的干涉条纹可以看作为移相干涉图；(g)模拟中使用的血红细胞模型；(h)使用基于快速傅里叶变换的相位恢复算法得到的相位；(i)使用10幅干涉图，利用主成分分析相位恢复算法得到的相位；(j)利用不同数量的干涉图，使用主成分分析相位恢复算法所消耗的时间。 

然而，如果仅关注以血红细胞模型中心为视场，通过通过截取该有效视场，可以得到如图1中(d)-(f)所示的干涉条纹。从图中可以看出，这些干涉图可以看作为一系列的移相干涉图。使用主成分分析相位恢复算法无需确定每幅干涉图的移相角，便可以恢复出生物样品的定量相位分布。图1中(g)为模拟中使用的血红细胞模型的相位分布，图1中(h)为使用主成分分析相位恢复算法得到的相位分布。为了比较这种方法与传统的相位恢复算法，图1中(i)展示了使用基于快速傅里叶变换相位恢复算法得到的血红细胞的定量相位分布。使用不同的相位恢复算法都能够清楚的得到血红细胞凹圆盘状的结构。因此，为了定量说明方法的优劣，这里，使用了相关系数来评价两种方法。其中，使用主成分分析相位恢复算法得到的结果为0.9999，而使用基于快速傅里叶变换相位恢复算法得到的结果仅有0.9868。结果说明该算法有较高的计算精度。此外，这里还比较了使用该算法与其它算法的计算效率。使用同一计算平台(Dell Inspiron，Core i5，2.40GHz，4G RAM)，计算100组移相干涉图，每组10幅移相干涉图(像素数256×256)，使用主成分分析相位恢复算法的时间为0.7729秒，而使用基于快速傅里叶变 换相位恢复算法则需要0.9984秒。可见，使用主成分分析相位恢复算法可以将计算效率提高约20％。此外，在模拟过程中，我们还计算了对于不同张数的移相干涉图，使用主成分分析相位恢复算法的计算时间，如图1中(j)所示，从图中可以看出，计算时间与使用的移相干涉图张数成正比。 

Claims (1)

1.一种主分成分析相位恢复算法，其特征在于：包括以下步骤：

步骤1重组干涉图：将每一张移相干涉图重组成行向量，然后将这N张移相干涉图的行向量按照列组合组成如下的x矩阵：x＝[x₁，x₂，...，x_N]^T，其中每一行都是每一张移相干涉图重组得到的一维数据，其长度为M＝N_x×N_y，N_x和N_y分别为移相干涉图x轴与y轴所占的像素点数，T代表为矩阵的转置；

步骤2获得背景分量m_x：矩阵m_x具有和矩阵x同样的维度，而矩阵m_x中所有元素的值均一致，其是矩阵x中所有元素的平均值；

步骤3计算协方差矩阵C：从矩阵x得到协方差矩阵C：C＝(x-m_x)(x-m_x)^T；

步骤4计算对角化协方差矩阵U：协方差矩阵C是一实对称N×N矩阵，因此该协方差矩阵可以如下式实现对角化：D＝UCU^T，矩阵D为对角化协方差矩阵而矩阵U是一个正交变换的矩阵，其大小为都为N×N；

步骤5得到主成分分量y：主成分便可以由矩阵D，矩阵x和矩阵m_x得到：y＝U(x-m_x)，其中，矩阵y的第一列和第二列便分别代表了主成分的正交特征值，其分别为I_c和I_s；

步骤6通过反正切函数求解相位：Φ＝arctan(I_s/I_c)。