1. 简介

Gerchberg-Saxton (GS) 算法是一种常用于相位恢复和光学成像的迭代算法。该算法最初由R.W. Gerchberg和W.O. Saxton于1972年提出 [1]，主要用于从强度测量数据中恢复相位信息。以下是该算法的基本步骤：

1. 初始化：
   • 准备输入图像的幅度信息（通常是从实验数据中获得的强度图像，取其平方根得到幅度）。
   • 初始化相位信息，通常可以设置为随机相位或者全零相位。
2. 频域迭代：
   • 对当前的图像（包含初始相位信息）进行傅里叶变换，得到频域图像。
   • 用实验测得的频域幅值替换频域图像的幅度，保留原相位信息。
   • 对替换后的频域图像进行逆傅里叶变换，回到空间域。
3. 空间域迭代：
   • 用实验测得的空间域幅值替换空间域图像的幅度，保留逆傅里叶变换后得到的相位信息。
   • 将更新后的图像再次进行傅里叶变换，进入下一次迭代。
4. 收敛判断：
   • 判断相位是否收敛，即相位变化是否在一定阈值范围内。
   • 如果达到收敛条件，则输出最终的相位信息；否则，返回步骤2继续迭代。

2. 算法描述

输入：

• $|x[n,m]|$：空间域的幅度信息
• $|X[k,l]|$​： 频域的幅度信息
• $ϵ$：误差阈值

输出：

• $z[n,m]$ -满足以下幅度约束的二维数组，即：$|z[n,m]|=|x[n,m]|$ 且 $|Z[k,l]|=|X[k,l]|$, 其中 $Z[k,l]$ 是 $z[n,m]$ 的离散傅里叶变换(DFT)。

初始化：

• 选择初始 $z_0[n,m]=|x[n,m]|\exp (j\varphi [n,m])$ (例如，使用随机相位$\varphi [n,m]$)。

一般步骤 (第$i$次迭代，$i=1,2,\dots )$ ：

1. ​ 对 $z_i[n,m]$ 进行二维傅里叶变换，得到 $Z_i[k,l]$。

2. ​ 保持当前傅里叶相位不变，但施加频域幅度约束：$Z_i^{\prime }[k,l]=|X[k,l]|\cdot Z_i[k,l]/|Z_i[k,l]|$。

3. ​ 对$Z_i^{\prime }[k,l]$进行二维逆傅里叶变换，得到$z_i^{\prime }[n,m]$。

4. ​ 保持当前空间相位不变，但施加空间域幅度约束：$z_{i+1}[n,m]=|x[n,m]|\cdot z_i^{\prime }[n,m]/|z_i^{\prime }[n,m]|\text{。}$

5. ​ 回到步骤1。

终止条件：直到误差
$$E_i=\sum _{k,l}{\Vert \left|Z_i[k,l]\right|-\left|X[k,l]\right|\Vert }^2\le ϵ$$

3. 混合输入输出（Hybrid Input-Output, HIO）算法

​ 混合输入输出（Hybrid Input-Output, HIO）算法是Gerchberg-Saxton (GS) 算法的一种改进版本，用于加快相位恢复的收敛速度，并解决GS算法容易陷入局部最优解的问题。HIO算法由Fienup于1982年提出 [2]，通过在迭代过程中引入一种新颖的更新策略，改善了相位恢复的性能。

3.1 HIO算法步骤

1. 初始化：

   • 与GS算法相同，准备输入图像的幅度信息和初始相位信息（随机相位或全零相位）。

2. 频域迭代：

   • 对当前的图像（包含当前相位信息）进行傅里叶变换，得到频域图像。
   • 用实验测得的频域幅值替换频域图像的幅度，保留原相位信息。
   • 对替换后的频域图像进行逆傅里叶变换，回到空间域。

3. 空间域更新：

   • 使用以下更新公式对空间域图像进行修正：
     $$z_{n+1}(x)=\{\begin{array}{ll}z(x)& \text{if}x\in \Omega \\ z_n(x)-\beta z(x)& \text{if}x\notin \Omega \end{array}$$
     其中，$z(x)$ 是逆傅里叶变换后的图像，$z_n(x)$ 是上一次迭代的图像，$\beta$ 是一个控制参数，通常取值在0.5到1之间，$\Omega$ 是已知幅值信息的区域。

4. 收敛判断：

   • 判断相位是否收敛，即相位变化是否在一定阈值范围内。
   • 如果达到收敛条件，则输出最终的相位信息；否则，返回步骤2继续迭代。

3.2 HIO算法的优势

1. 更快的收敛速度：HIO算法通过引入混合输入输出的更新策略，有效避免了GS算法的局部最优解问题，通常能够更快地收敛到全局最优解。
2. 更好的鲁棒性：由于HIO算法可以在不满足约束条件的区域进行负反馈修正，因此在处理复杂相位恢复问题时表现更为稳定。

3.3 算法描述

输入：

• $|x[n,m]|$：空间域的幅度信息
• $|X[k,l]|$： 频域的幅度信息
• $ϵ$：误差阈值
• $\beta$ ：HIO算法的反馈参数

输出：

• $z[n,m]$ - 满足以下幅度约束的二维数组，即：$|z[n,m]|=|x[n,m]|$ 且 $|Z[k,l]|=|X[k,l]|$, 其中 $Z[k,l]$ 是 $z[n,m]$ 的离散傅里叶变换(DFT)。

初始化：

• 选择初始 $z_0[n,m]=|x[n,m]|\exp (j\varphi [n,m])$ (例如，使用随机相位$\varphi [n,m]$)。

一般步骤 (第$i$次迭代，$i=1,2,\dots ):$

1. 对$z_i[n,m]$进行二维傅里叶变换，得到$Z_i[k,l]$。

2. 保持当前傅里叶相位不变，但施加频域幅度约束：$Z_i^{\prime }[k,l]=|X[k,l]|\cdot Z_i[k,l]/|Z_i[k,l]|$。

3. 对$Z_i^{\prime }[k,l]$进行二维逆傅里叶变换，得到$z_i^{\prime }[n,m]$。

4. 保持当前空间相位不变，但施加空间域幅度约束，且对不满足约束的区域进行反馈修正：
   $$z_{i+1}[n,m]=\{\begin{array}{ll}|x[n,m]|\cdot z_i^{\prime }[n,m]/|z_i^{\prime }[n,m]|& \text{if}|x[n,m]|\ne 0\\ z_i[n,m]-\beta z_i^{\prime }[n,m]& \text{if}|x[n,m]|=0\end{array}$$

5. 回到步骤1。

终止条件：直到误差
$$E_i=\sum _{k,l}\Vert |Z_i[k,l]|-|X[k,l]|{\Vert }^2\le ϵ$$

4. 算法实现与对比

import numpy as np
import matplotlib.pyplot as plt
from skimage.data import camera, checkerboard
from skimage.transform import resize

# 加载skimage库中的示例图像
amplitude_image = camera()
phase_image = checkerboard()

# 如果振幅图像和相位图像的大小不同，调整相位图像的大小
if amplitude_image.shape != phase_image.shape:
    phase_image = resize(phase_image, amplitude_image.shape, anti_aliasing=True)

# 对振幅和相位图像进行填充
padding_width = 50
amplitude_image = np.pad(amplitude_image, pad_width=padding_width, mode='constant', constant_values=0)
phase_image = np.pad(phase_image, pad_width=padding_width, mode='constant', constant_values=0)

# 将图像归一化到[0, 1]范围内
amplitude_image = amplitude_image / np.max(amplitude_image)
phase_image = phase_image / np.max(phase_image)

# 生成复数对象
complex_obj = amplitude_image * np.exp(1j * phase_image)

# 进行傅里叶变换
ft_obj = np.fft.fft2(complex_obj)
abs_ft_obj = np.abs(ft_obj)

# 初始化相位图像
initial_phase = np.angle(np.exp(1j * np.random.rand(*amplitude_image.shape)))


# Gerchberg-Saxton (GS) 算法
def gerchberg_saxton(amplitude, initial_phase, num_iterations, epsilon):
    z = amplitude * np.exp(1j * initial_phase)
    errors = []
    for i in range(num_iterations):
        # 正向傅里叶变换
        Z = np.fft.fft2(z)

        # 计算误差并检查终止条件
        error = np.sum((np.abs(Z) - abs_ft_obj) ** 2)
        errors.append(error)
        if error <= epsilon:
            print(f'GS algorithm converged after {i + 1} iterations with error {error:.6f}')
            break

        # 在傅里叶空间中施加振幅约束
        Z = abs_ft_obj * Z / np.abs(Z)
        # 逆傅里叶变换
        z_prime = np.fft.ifft2(Z)
        # 在实空间中施加振幅约束
        z = amplitude * z_prime / np.abs(z_prime)

    final_phase = np.angle(z)
    return final_phase, errors


# Hybrid Input-Output (HIO) 算法
def hio(amplitude, initial_phase, num_iterations, beta, epsilon):
    z = amplitude * np.exp(1j * initial_phase)
    errors = []
    for i in range(num_iterations):
        # 正向傅里叶变换
        Z = np.fft.fft2(z)

        # 计算误差并检查终止条件
        error = np.sum((np.abs(Z) - abs_ft_obj) ** 2)
        errors.append(error)
        if error <= epsilon:
            print(f'HIO algorithm converged after {i + 1} iterations with error {error:.6f}')
            break

        # 在傅里叶空间中施加振幅约束
        Z = abs_ft_obj * Z / np.abs(Z)
        # 逆傅里叶变换
        z_prime = np.fft.ifft2(Z)
        # 在实空间中施加振幅约束，并引入反馈机制
        mask = (amplitude != 0)
        z = np.where(mask,
                     amplitude * z_prime / np.abs(z_prime),
                     z - beta * z_prime)

    final_phase = np.angle(z)
    return final_phase, errors


# 参数设置
num_iterations = 50
beta = 0.9
epsilon = 1e-6

# 应用GS算法
gs_phase, gs_errors = gerchberg_saxton(amplitude_image, initial_phase, num_iterations, epsilon)

# 应用HIO算法
hio_phase, hio_errors = hio(amplitude_image, initial_phase, num_iterations, beta, epsilon)

# 设置绘图的字体和字号
font = {'family': 'serif',
        'weight': 'normal',
        'size': 20}
plt.rc('font', **font)

# 显示结果
plt.figure(figsize=(12, 6))

plt.subplot(1, 3, 1)
plt.title('Original Phase')
plt.imshow(phase_image, cmap='gray')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.title('GS Phase')
plt.imshow(gs_phase, cmap='gray')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.title('HIO Phase')
plt.imshow(hio_phase, cmap='gray')
plt.axis('off')

plt.tight_layout()
plt.show()


# 绘制误差下降曲线
plt.figure(figsize=(12, 6))
plt.plot(gs_errors, label='GS Errors')
plt.plot(hio_errors, label='HIO Errors')
plt.xlabel('Iteration')
plt.ylabel('Error')
plt.yscale('log')
plt.title('Error Convergence')
plt.legend()
plt.show()

循环50次的结果对比如下：

5. 总结

GS算法和HIO算法都是相位检索的有效方法，各有优缺点。GS算法简单易实现，但在复杂情况下收敛较慢且易陷入局部极小值。HIO算法通过引入反馈机制加速收敛，并增强鲁棒性，但其实现相对复杂且对参数选择敏感。在实际应用中，可以根据具体问题的需求选择合适的算法，或者结合两种算法的优点进行混合使用。

参考文献

1. R. Gerchberg and W. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik, vol. 35, p. 237, 1972.
2. Fienup J R. Phase retrieval algorithms: a comparison[J]. Applied optics, 1982, 21(15): 2758-2769.