Gerchberg-Saxton (GS) 和混合输入输出（Hybrid Input-Output, HIO）算法_gerchberg-saxton算法

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 ] = { ∣ x [ n , m ] ∣ ⋅ z i ′ [ n , m ] / ∣ z i ′ [ n , m ] ∣ if ∣ x [ n , m ] ∣ ≠ 0 z i [ n , m ] − β z i ′ [ n , m ] if ∣ x [ n , m ] ∣ = 0 z_{i+1}[n,m]=

   $$\begin{cases}|x[n,m]|\cdot z_i'[n,m]/|z_i'[n,m]|&\text{if}|x[n,m]|\neq0\\z_i[n,m]-\beta z_i'[n,m]&\text{if}|x[n,m]|=0\end{cases}$$ zi+1​[n,m]={∣x[n,m]∣⋅zi′​[n,m]/∣zi′​[n,m]∣zi​[n,m]−βzi′​[n,m]​if∣x[n,m]∣=0if∣x[n,m]∣=0​

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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

循环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.