浅谈计算成像：计算相位成像（二）

前言：

计算成像这一个领域我断断续续关注已有一段时间，但是实操去进行coding的经验还相对有限，本文及后续系列文章旨在分享我对计算成像领域的认识、相关算法总结，相关器件和光学系统设计等。长期关注我的读者可能已经注意到，我的专业领域主要集中在器件、物理、硬件层面，虽然计算成像我断断续续关注已有一段时间，对于全息重建、ptycho、TIE，横向剪切干涉、单像素等有一些摸索，但是实操经验还相对有限，对于算法和系统层面的深入理解也有所欠缺，如果本文存在任何疏漏，恳请各位不吝赐教！也希望相关同行能多多交流~

/00  计算相位成像算法概述/

要实现相位信息获取，最直接的方法是通过干涉将相位信息转化为强度信息，包括干涉全息（离轴干涉和同轴干涉）、横向剪切干涉等，此外还可通过非干涉的方式实现相位复原，计算相位成像算法中最常见的非干涉算法包括迭代法（典型代表是GS算法，ptychography PIE算法）和直接法（TIE算法）。

常见的相位成像技术，参考文献9

本系列文章会对主要的的相位成像算法进行一个简要介绍并分享网上的开源代码供读者学习，包括 TIE,ptychography，Fourier ptycho，GS，inline holography等，这些算法都是计算相位成像领域的关键技术，它们用于从各种测量数据中恢复物体的相位信息。以下是每种技术的简要介绍：

TIE (Transport of Intensity Equation)：

• 基于光强传输方程，利用了光波在传播过程中的相位与强度之间的关系实现对相位进行直接求解。

Ptychography：

• 一种基于扫描成像的相位恢复技术，通过移动样品或探测器来实现多次扫描，并控制实现每次扫描区域保证一定的重叠度，多张扫描图在重建的迭代过程中关联起来，通过分析这些强度图，从而可以重建样品的高分辨相位信息。

Inline Holography：

• 即同轴全息术，是一种全息成像技术，其中参考光束和物光束在同一轴线上。无透镜和共路干涉系统中用到较多。

Off-axis Holography：

• 非同轴全息术，与同轴全息术不同，参考光束和物光束在空间中有一定的角度偏移。这种偏移有助于避免重建过程中的零级衍射和双图像问题。

QWLSI (Quantitative Wide-field Light Sheet Imaging)：

• 四波横向剪切技术，通过光栅将波前复制四分并”横向位移“后与其余复制波前两两干涉。从而通过这一自干涉实现相位信息获取。

DPC (Differential Phase Contrast)：

• 差分相位成像，与QWLSI有一定的类似性。

Fourier Ptychography：

在傅里叶面进行扫描的叠层成像。

/01 G-S算法/

Gerchberg-Saxton（GS）算法是一种通过强度测量获取相位信息的迭代式相位恢复算法，这种算法最初由Roy Gerchberg和Warren Saxton在1972年提出。

通过以两个不同平面上(可以是入射光场面和成像面，或者图像平面和远场傅里叶平面，最常选择的是后者，即以空间域和频率域作为两个平面)获得的强度信息作为已知和约束条件，通过正向传播和反向传播，迭代恢复波前的相位，从而获得复振幅图像信息。经典的GS算法是利用傅里叶变换及其逆变换来计算光场在两个平面之间的传播，其原理是基于角谱传播。

GS算法的计算过程如图所示，（1）先通过估算或者随机生成给出一个初始相位分布，（2）以入射面振幅信息作为已知，乘以相位分布获得入射面复振幅信息；（3）将入射面复振幅信息通过正向传播到输出面，获得输出面复振幅信息，（4）将输出面振幅信息作为已知约束，替换掉传播得到的复振幅分布中的实部振幅项，（5）以修改后的复振幅图像反向传播到入射面获得新的复振幅图像，并重复（2-5）,以此反复迭代，实现收敛。可以看到，GS求解是迭代式计算过程，因此是一种误差减小算法。

GS算法的伪代码如下，其很好的概括了GS算法的核心过程：

  A = IFT(Target)
  while error criterion is not satisfied
    B = Amplitude(Source) * exp(i*Phase(A))
    C = FT(B)
    D = Amplitude(Target) * exp(i*Phase(C))
   A = IFT(D)
  end while
  Retrieved_Phase = Phase(A)

GS算法是很多迭代式计算相位成像算法的基础，在编码成像，叠层成像、全息成像等里也能看到GS算法的影子。虽然 GS 迭代算法可以简单的方式得到相位分布信息，但是由于算法本身固有的缺陷，GS算法在迭代过程中可能会陷入局部最优解。Gerchberg-Saxton (GS) 算法的衍生算法是为了解决原始GS算法在某些情况下收敛速度慢或可能不收敛的问题而开发的。这些衍生算法通过引入新的迭代策略、优化误差函数、在更新复振幅时引入权重函数等方式来改进相位恢复的过程。一些GS衍生算法包括：weighted Gerchberg-Saxton，Adaptive-Adaptive算法，Fienup算法，HIO 算法等

加权GS算法：会先进行数轮GS算法迭代，之后会将重建图像平面中的限制条件替换为目标图像的振幅乘以一个分数形式的权重系数，然后继续运算直到满足一定的误差条件。以原始图像与重建图像间的归一 化标准差作为评判标准，加权ＧＳ算法比ＧＳ算法取得了一定程度的提升且收敛速率会有所提高。

Fienup算法：在重建图像平面上的限制条件有所不同，不再是将该平面上的光波振幅直接限制为目标图像振幅，而是设为上次迭代结果与原始图像振幅的一个函数与上次迭代结果之和。

AA算法：AA算法（Adaptive-Adaptive method）是在重建图像平面上的限制条件是将振幅分布设为原始图像与上次迭代结果的线性组合。

HIO算法：Hybrid Input-Output (HIO) 算法在更新过程引入控制权重参数调控约束和上次一迭代的占比，从而使得约束的限制作用得以”放松“，核心思想是在迭代过程中引入反馈动态调整输入和输出约束，以加速收敛并提高相位恢复的准确性和鲁棒性。

HIO算法原理

广义的Gerchberg-Saxton算法可以用于任何已知部分约束（以测量数据或先验信息的形式）在两个域中的每一个中的问题，这两个域通常是物体（或图像）域和傅里叶域。只需在这两个域之间来回变换，在一个域中满足约束后再返回到另一个域。此外还可以把两个平面延申到多个平面的迭代。

Multiple-Phase Retrieval 算法，来自参考文献5

Multiple-Phase Retrieval 算法，来自参考文献5

除了上述提及的算法外，GS还有很多衍生算法，这里就暂时不一一列举了。后面有机会在总结介绍一下。

相关GS算法代码

目前网上有很多GS算法代码，大家可以通过代码演示感受GS算法的实现原理。

python代码：1.https://github.com/kurokuman/Gerchberg-Saxton-algorithm

2.https://blog.csdn.net/xy_optics/article/details/139079428

3.https://github.com/madisousa/Phase-Retrieval-Algorithms-

matlab代码：1.https://ww2.mathworks.cn/matlabcentral/fileexchange/65979-gerchberg-saxton-algorithm

/02 TIE算法/

TIE（transport of intensity Equation） 算法基于光强传输方程，利用了光波在传播过程中的相位与强度之间的关系实现对相位进行直接求解，光强传输方程是一个描述光波的轴向光强变化量和相位定量耦合的偏微分方程。TIE方程指出，光波的相位信息向强度信息的转化不仅仅依赖于干涉，光波自身的传播效应就是一种自发的光波-相位转化过程。光强传输方程可以表示为：

光强传输方程是旁轴近似下的推导近似，其可以通过将光波复振幅代入亥姆霍兹方程，然后分离实部虚部得到；也可通过菲涅耳定律再小传播距离下得到还可通过傍轴近似下坡印廷定理推导得到。

TIE算法原理，来自参考文献9

从物理本质上讲，光强传输方程本质是能量守恒定律，相位信息表示了波前信息，波前分布决定了场在沿着垂直传播方向上的分布情况。根据光强传输方程可以看到，其光强在传播方向上的变化情况与其波前相位相联系，也就与其横向场分布情况相联系，也就是说一个波束在传播方向上光强的减少必然伴随着其能量在垂直传播方向上的再分布，也就是相位的变化。

TIE的求解方法较多，包括格林函数法，FFT,DCT，迭代DCT，US-TIE等。左超教授在其文章中做了很多相关算法的归纳总结。

光强传输方程求解算法对比

由于TIE算法可以直接求解相位，其也可和迭代法结合，以TIE直接求解的相位信息作为初始量，再通过迭代提高恢复的相位信息的准确性。

iTIE+角谱迭代的混合算法

参考文献：

1、Gerchberg R W. A practical algorithm for the determination of plane from image and diffraction pictures[J]. Optik, 1972, 35(2): 237-246.

2、Fienup J R. Phase retrieval algorithms: a comparison[J]. Applied optics, 1982, 21(15): 2758-2769.

3、Di Leonardo, Roberto, Francesca Ianni, and Giancarlo Ruocco. "Computer generation of optimal holograms for optical trap arrays." Optics Express 15.4 (2007): 1913-1922.

4、同轴数字全息中的相位恢复算法[J]. Chinese Journal of Lasers, 2014, 41(2): 209006-1.

5、Modified Gerchberg–Saxton (G-S) Algorithm and Its Application

6、Phase hologram optimization with bandwidth constraint strategy for speckle-free optical reconstruction

7、卜浩祯,焦述铭.纯相位全息图优化算法[J].液晶与显示,2021,36(06):810-826.

8、Phase retrieval algorithms: a comparison

9、 Transport of intensity equation: a tutorial

10、《光强传输方程》 陈钱、左超著

11、Jialin Zhang, Qian Chen, Jiasong Sun, Long Tian, and Chao Zuo, "On a universal solution to the transport-of-intensity equation," Opt. Lett. 45, 3649-3652 (2020)

12、一种混合迭代算法在大离焦距离TIE相位恢复中的应用，光学学报，2016