MATLAB实现GS算法：从相位恢复原理到光学成像仿真实践-平芜编程栈

简介：本资源是一套面向光学专业本科生与初学者的Matlab仿真教学工具包，聚焦Gerchberg-Saxton（GS）迭代算法在光学相位恢复与波前重建中的原理实现与可视化验证，专为中国科学技术大学光学课程作业设计，解决“仅知强度分布、如何反演相位信息”这一核心教学难点。压缩包共27个文件，含10个核心Matlab脚本（如Gerchberg_Saxton_Algorithm.m、FT2Dc.m、RP.m等）、5幅算法中间过程及结果图（.fig）、5张关键输出图像（.jpg，含目标图、重建像、相位分布等），以及说明文档（.txt）、理论补充（.docx）和二进制样例数据（.bin），总大小39.03MB。已有75人学习下载，资源结构清晰：主程序模块化、参数可调、每步迭代可视化，附赠文档详述算法推导与计算全息关联，配套fig与jpg直观呈现空间域与频域交替约束效果，便于理解收敛行为与重建质量评估。

1. 项目概述：从课程作业到光学成像的核心技术

看到这个项目标题，很多光学、物理或者电子信息专业的同学应该会心一笑。这几乎是中国科学技术大学光学相关课程里一个经典的“大作业”或课程设计题目。它的核心，是利用Gerchberg-Saxton（GS）迭代算法，在MATLAB环境中仿真实现光学相位恢复与波前重建。简单来说，我们平时用相机拍照，记录的是光波的强度信息，也就是明暗，但丢失了至关重要的相位信息。而相位，决定了光波的波前形状，包含了物体深度、折射率分布等关键信息。GS算法就是一种非常巧妙的方法，它允许我们仅通过采集到的光强信息（比如两张不同位置的光斑图），通过数学迭代，反推出丢失的相位，从而重建出完整的光波场。

这个项目绝不仅仅是为了完成作业。它在现代光学的前沿领域有着举足轻重的地位。计算全息、自适应光学（用于校正大气湍流对天文望远镜的影响）、相位显微成像（用于观察透明生物样品）、光学加密以及新型的光场调控等领域，其底层核心都离不开相位恢复技术。通过亲手用MATLAB实现GS算法，你不仅能深刻理解傅里叶光学、衍射理论的核心概念，更能掌握一套解决实际逆向问题的计算工具。这个过程，是从理论公式走向实际应用的关键一步，你会遇到算法收敛性、初始猜测敏感性、噪声影响等一系列非常实际的问题，而解决它们的过程，正是能力提升的体现。

2. 核心原理拆解：GS算法如何“无中生有”出相位

要理解GS算法，我们得先搞清楚光波场是什么。一个单色光波，在空间某一点的复振幅可以表示为 U(x, y) = A(x, y) * exp(i * φ(x, y))。其中A是振幅（amplitude），直接对应我们探测到的光强 I = |A|²；φ就是相位（phase），它描述了光波在该点相对于某个参考点的振动状态。探测器（如CCD相机）是“相位盲”的，只能记录光强I，宝贵的相位φ在记录瞬间就丢失了。

GS算法的天才之处在于，它利用了两个域（输入面和输出面，通常是物面和频谱面）之间的傅里叶变换关系，以及我们对这两个面上光强分布的已知信息，通过迭代来逼近真实的相位。其核心流程是一个闭合的迭代循环，我们可以把它想象成一个不断自我修正的猜谜游戏。

算法的基本迭代步骤（一个循环）如下：

初始猜测：我们从物面（输入面）开始。已知物面的光强分布 I_object（比如一个物体的透射率）。但对于相位，我们一无所知。因此，我们随机生成一个初始相位猜测 φ_guess（通常在0到2π之间均匀分布），与已知的振幅 sqrt(I_object) 结合，构造一个初始的复振幅场：U1 = sqrt(I_object) * exp(i * φ_guess)。
正向传播（约束变换）：将构造的物面场 U1 传播到输出面（通常是焦平面或衍射场）。这个传播过程在傍轴近似下，通过一次傅里叶变换（FFT）来实现：U2 = FFT(U1)。此时，U2 同时具有振幅 A2 和相位 φ2。
施加输出面约束：这是我们拥有的另一个已知信息。在输出面（例如透镜的后焦面），我们通过实验测量得到了该处的光强分布 I_measured。GS算法的关键操作来了：我们保留计算得到的相位 φ2，但用测量得到的振幅 sqrt(I_measured) 替换掉计算出的振幅 A2。即，更新输出面场为：U2‘ = sqrt(I_measured) * exp(i * φ2)。这个步骤强制让计算出的波场满足实际的物理测量结果。
反向传播（逆变换）：将更新后的输出面场 U2‘ 传播回物面。通过逆傅里叶变换（IFFT）实现：U1‘ = IFFT(U2‘)。此时，我们得到了一个新的物面复振幅场，包含振幅 A1‘ 和相位 φ1‘。
施加输入面约束：同样，我们保留新计算出的相位 φ1‘，但用我们最初已知的物面振幅 sqrt(I_object) 替换掉计算出的振幅 A1‘。即，更新物面场为：U1_new = sqrt(I_object) * exp(i * φ1‘)。这个步骤强制让波场满足物面的已知条件。

至此，一个迭代循环完成。U1_new 将作为下一个循环的起点，重复步骤2到5。如此反复迭代，每一次循环都让计算出的波场在输入面和输出面同时满足各自的振幅（强度）约束。理论上，随着迭代次数增加，计算得到的相位 φ1‘ 会收敛到真实的物面相位。

注意：这里描述的是最基本的“两平面”GS算法。在实际应用中，根据不同的物理设置（如多平面相位恢复、基于角谱的衍射传播等），正向和反向传播算子可能不是简单的FFT/IFFT，而是更复杂的衍射积分计算。但核心的“变换-施加振幅约束-反变换”的迭代思想是完全一致的。

为什么这样迭代会收敛？直观上理解，这个算法在寻找一个解，这个解在经过特定变换（如傅里叶变换）后，其振幅分布能与我们在另一个平面上测量到的结果匹配。迭代过程就像在解一个巨大的方程组，每次施加约束都是将解向可行解集合投影。GS算法本质是一种交替投影算法，在两个约束集合（满足输入面振幅的复函数集合，和满足输出面振幅的复函数集合）之间来回投影，最终希望找到两个集合的交点，即同时满足两个约束的解。

3. MATLAB仿真系统设计与实现要点

用MATLAB搭建这个仿真系统，目的很明确：在已知“真实”相位的情况下，模拟光波的传播和探测过程，生成仿真的“测量光强”；然后，假装我们不知道这个相位，仅用“测量光强”和GS算法来恢复它，最后与“真实”相位对比，评估算法性能。这个过程让我们可以在完全可控的环境下研究算法的特性。

3.1 系统模块分解

一个结构清晰的仿真系统通常包含以下几个模块：

参数定义与初始化模块：设定仿真环境的基本参数。这包括：
- 网格参数：仿真区域的尺寸（Lx,Ly，单位：米或毫米），采样点数（Nx,Ny）。采样点数一般取2的整数次幂（如256，512，1024），以便利用FFT的高效性。网格间距dx = Lx/Nx。
- 光学参数：波长（lambda），传播距离（z，如果是角谱传播），透镜焦距（f，如果是4f系统）等。
- 目标物体：定义输入面的振幅分布A_obj（如一个圆形孔径、一个字母图案）和真实的相位分布phi_true（这是我们希望恢复的，可以是一个简单的球面波相位、一个随机相位板，或者一个复杂的像差分布）。
- 生成模拟的真实波场：U_true = A_obj .* exp(1i * phi_true)。
正向传播模拟模块：这个模块模拟真实物理过程，即光从物面传播到探测面。
- 根据实验光路，选择合适的传播模型。对于最常见的4f滤波系统或焦面探测，使用傅里叶变换：U_prop = fft2(U_true)或fftshift(fft2(fftshift(U_true)))（注意零频分量移到中心）。
- 对于非焦面的自由空间衍射，可能需要使用角谱传播或菲涅尔衍射积分。角谱传播更精确，其核心是：AS = fft2(U_true); H = exp(1i * 2*pi*z / lambda .* sqrt(1 - (lambda*fx).^2 - (lambda*fy).^2)); U_prop = ifft2(AS .* H);其中fx,fy是空间频率坐标。
- 计算探测面的光强（模拟测量结果）：I_measured = abs(U_prop).^2。为了更真实，可以在此处加入噪声，如高斯噪声：I_measured_noisy = imnoise(I_measured, 'gaussian', 0, 0.01)（添加方差为0.01的高斯噪声）。
GS算法迭代恢复模块：这是系统的核心。
- 初始化：设定最大迭代次数max_iter，收敛阈值threshold。生成随机的初始相位phi_guess = 2*pi*rand(Nx, Ny)。构造初始估计场：U_est = A_obj .* exp(1i * phi_guess)。
- 迭代循环：for iter = 1:max_iter
  - 正向传播至探测面：U_prop_est = Propagation_Model(U_est)，使用与正向模拟模块相同的传播模型。
  - 施加探测面振幅约束：A_prop_measured = sqrt(I_measured)。U_prop_est_constrained = A_prop_measured .* exp(1i * angle(U_prop_est))。angle()函数用于获取相位。
  - 反向传播回物面：U_est_back = Inverse_Propagation_Model(U_prop_est_constrained)。
  - 施加物面振幅约束：U_est_new = A_obj .* exp(1i * angle(U_est_back))。
  - 计算误差：评估当前估计与约束的匹配程度。常用误差函数是探测面振幅误差：error(iter) = sum(sum(abs(abs(U_prop_est) - A_prop_measured).^2)) / sum(sum(A_prop_measured.^2))。
  - 判断收敛：如果error(iter) < threshold，跳出循环。
  - 更新估计场：U_est = U_est_new。
- 输出：循环结束后，从最终的U_est中提取恢复的相位：phi_recovered = angle(U_est)。
结果可视化与评估模块：这是验证算法效果的关键。
- 图像显示：并排显示真实相位phi_true、恢复相位phi_recovered以及两者的差值（残差）图。由于相位是周期性的（模2π），通常使用imagesc显示，并注意相位包裹（phase wrapping）问题，可能需要使用angle()函数直接显示（范围 -π 到 π），或使用wrapToPi等函数处理。
- 误差分析：绘制迭代误差error随迭代次数的下降曲线，观察收敛速度和稳定性。
- 定量指标：计算恢复相位与真实相位之间的均方根误差（RMSE）、峰值信噪比（PSNR）或结构相似性（SSIM）。对于存在相位包裹的情况，比较前需要先进行解包裹（phase unwrapping）或计算梯度。

3.2 关键代码片段与实操心得

下面给出一些核心步骤的MATLAB代码片段和编写时的注意事项：

1. 生成带有相位的测试物体：

% 参数定义 N = 512; % 采样点数 L = 10e-3; % 物理尺寸 10mm dx = L/N; x = linspace(-L/2, L/2, N); [X, Y] = meshgrid(x, x); % 振幅物体：一个圆形孔径 radius = 1e-3; % 1mm半径 A_obj = sqrt((X.^2 + Y.^2) <= radius^2); % 振幅，1代表透光，0代表不透光 % 真实相位：一个球面波相位（模拟透镜）叠加一些像差 lambda = 632.8e-9; % He-Ne激光波长 k = 2*pi/lambda; R = 0.5; % 球面波曲率半径 phi_true = mod(k/(2*R) * (X.^2 + Y.^2) + 0.5*sin(4*pi*X/L) .* cos(2*pi*Y/L), 2*pi); % 使用mod函数将其范围限制在[0, 2π)，避免数值过大。实际仿真中，相位值范围不重要，重要的是其空间变化。 % 第二部分加入了一个简单的正弦型像差。 % 构造真实复振幅场 U_true = A_obj .* exp(1i * phi_true);

实操心得：在定义相位phi_true时，要特别注意其空间频率不能超过奈奎斯特频率（即每个采样点间隔内，相位变化不能超过π），否则会产生混叠，导致仿真失真。对于快速变化的相位，需要确保采样足够密集。

2. 角谱传播函数（用于自由空间衍射模拟）：

function Uout = AngularSpectrumPropagation(Uin, lambda, dx, dy, z) % Uin: 输入场 % lambda: 波长 % dx, dy: 输入面采样间隔 % z: 传播距离 % Uout: 输出场 [Ny, Nx] = size(Uin); % 生成空间频率坐标 fx = linspace(-1/(2*dx), 1/(2*dx), Nx); fy = linspace(-1/(2*dy), 1/(2*dy), Ny); [FX, FY] = meshgrid(fx, fy); % 计算角谱传递函数 H = exp(1i * 2*pi/lambda * z .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 对于频率分量满足 (lambda*fx)^2 + (lambda*fy)^2 > 1 的部分，sqrt结果为虚数，对应倏逝波，通常直接置零或做衰减处理。 H(sqrt((lambda*FX).^2 + (lambda*FY).^2) >= 1) = 0; % 进行角谱传播 Uin_f = fft2(fftshift(Uin)); % 将输入场频谱零频移到角落，便于与H相乘 Uout_f = Uin_f .* H; Uout = ifftshift(ifft2(Uout_f)); % 反变换并移回中心 end

注意：角谱传播是精确的标量衍射模型，但计算时要注意传递函数H中根号内的项为负的情况（对应高频倏逝波），处理不当会引入数值不稳定。通常的做法是将其滤除（置零），如上面代码所示。这相当于一个低通滤波，是物理合理的。

3. GS算法迭代核心循环：

max_iter = 200; threshold = 1e-6; error_history = zeros(max_iter, 1); % 初始猜测：随机相位 phi_guess = 2*pi * rand(N, N); U_est = A_obj .* exp(1i * phi_guess); % 已知的测量面振幅（模拟得到，已加噪声） A_measured = sqrt(I_measured_noisy); for iter = 1:max_iter % 1. 正向传播 (使用角谱传播或FFT，需与正向模拟一致) U_prop_est = AngularSpectrumPropagation(U_est, lambda, dx, dx, z); % 2. 施加测量面振幅约束 A_prop_est = abs(U_prop_est); U_prop_est_constrained = A_measured .* exp(1i * angle(U_prop_est)); % 3. 反向传播 U_est_back = AngularSpectrumPropagation(U_prop_est_constrained, lambda, dx, dx, -z); % 注意距离取负 % 4. 施加物面振幅约束 U_est_new = A_obj .* exp(1i * angle(U_est_back)); % 5. 计算误差（测量面振幅误差） error_history(iter) = sum(sum((A_prop_est - A_measured).^2)) / sum(sum(A_measured.^2)); % 6. 判断收敛 if error_history(iter) < threshold fprintf('在 %d 次迭代后收敛。\n', iter); break; end % 7. 更新估计 U_est = U_est_new; end phi_recovered = angle(U_est);

实操心得：反向传播函数理论上应该是正向传播的逆运算。对于角谱传播，传递函数是exp(-i*...)，即距离取负-z。对于FFT模型（4f系统），反向传播就是IFFT。务必保持正向模拟和迭代恢复中使用完全相同的传播模型和参数，否则算法可能无法收敛到正确解。

4. 影响算法性能的关键因素与调优策略

GS算法原理简洁，但实际应用中其性能受到多种因素影响。理解这些因素并学会调优，是完成高质量相位恢复的关键。

4.1 初始相位猜测

GS算法本质上是一个非线性优化过程，初始值会影响收敛速度和最终结果，尤其是在解不唯一或存在噪声的情况下。

随机相位：最常用的方法。2*pi*rand()生成均匀分布的随机相位。优点是简单，有助于算法跳出局部极小值。缺点是每次运行结果可能略有不同，收敛所需的迭代次数不稳定。
常数相位（零相位）：zeros(N, N)。作为最简单的猜测，在物体相位变化平缓且信噪比较高时，可能收敛很快。但如果真实相位变化剧烈，零相位猜测可能引导算法收敛到一个错误的解（局部极小值）。
先验信息：如果对物体相位有部分了解（例如，知道它是一个缓慢变化的像差），可以将其作为初始猜测，能显著加速收敛并提高准确性。

策略：在仿真中，可以尝试多种初始猜测，观察收敛结果的一致性。在实际实验中，如果没有先验信息，通常采用随机相位，并多次独立运行算法，选取误差最小的结果，或对多次结果进行平均，以增强鲁棒性。

4.2 振幅约束的施加与支持域（Support Constraint）

在基础GS算法中，我们对物面的振幅约束是严格的A_obj。但在很多情况下，我们只知道物体的大致范围（支持域），而不知道其内部精确的振幅分布。例如，在晶体学中，只知道分子的大致形状。

紧支持域：定义一个二值掩膜support，在物体可能存在的区域为1，其他区域为0。物面约束变为：U_est_new = (A_obj .* exp(1i*angle(U_est_back))) .* support + (U_est_back .* (1-support)) * beta。这里beta是一个松弛因子（0<beta<1），用于在非支持域缓慢衰减场，而不是直接置零，这有助于提高收敛性，这就是著名的混合输入输出（HIO）算法的核心思想。
松弛振幅约束：对于物面振幅，如果不是精确已知，可以施加一个范围约束，而不是固定值。

策略：如果你的项目涉及非均匀振幅物体或支持域信息，强烈建议实现HIO或其变种（如OSP、RAAR算法），它们比纯GS算法具有更强的抗噪声能力和更广的收敛域。

4.3 噪声的影响与处理

实验测量的光强I_measured必然包含噪声（散粒噪声、读出噪声等）。噪声会破坏振幅约束的有效性，导致GS算法振荡或不收敛到正确解。

现象：在噪声较大时，误差曲线error_history不会平稳下降到零，而是在某个值附近波动。恢复的相位图中会出现明显的随机斑点或条纹。
处理方法：
1. 多次独立运行取平均：用不同的随机初始相位多次运行GS算法，将最终恢复的相位进行平均，可以有效抑制随机噪声带来的波动。
2. 正则化：在迭代过程中加入正则化项，惩罚相位的高频变化（即假设相位是平滑的）。例如，在每次迭代后对恢复的相位进行轻微的高斯滤波。
3. 使用更鲁棒的算法：如前所述的HIO算法对噪声的容忍度比GS更高。
4. 数据预处理：对测量的光强图像进行平滑滤波或降噪处理（如小波去噪、非局部均值去噪），但要注意不能过度平滑而损失细节。

4.4 收敛性与停滞问题

GS算法可能陷入停滞（stagnation），即误差不再下降，但恢复的结果明显不正确。这通常是因为算法陷入了局部极小值。

诊断：观察误差曲线，如果曲线在迭代几十次后早早地进入平台期，且恢复的相位与真实相位相差甚远，可能就是陷入了停滞。
解决方案：
1. 改变初始猜测：换一个随机种子重新开始。
2. 引入扰动：在迭代过程中，每隔一定次数，人为地对当前估计的相位加入一个小的随机扰动，帮助算法跳出局部极小。这被称为“扰动GS”或“随机梯度”思想的一种简单实现。
3. 切换算法：当GS停滞时，切换到HIO算法迭代几十次，然后再切换回GS，往往能打破僵局。这种GS-HIO交替迭代的策略非常有效。
4. 多平面相位恢复：如果条件允许，采集多个不同传播距离（或离焦量）下的光强图像。对每个平面的光强同时施加约束，大大增加了算法的约束条件，能显著改善收敛性和唯一性。这就是著名的多平面迭代相位恢复技术。

5. 仿真结果分析、可视化与报告撰写

完成了算法实现和迭代，如何呈现你的工作成果至关重要。一份优秀的课程报告或项目总结，需要清晰、直观地展示仿真过程和结论。

5.1 结果可视化技巧

相位显示：直接使用imagesc(phi_recovered)显示相位，由于相位通常被包裹在[-π, π]或[0, 2π]区间，会呈现出明显的2π跳变条纹。为了更直观地看到相位的连续变化，需要进行相位解包裹。

% 使用MATLAB的 unwrap函数（适用于1D，对2D效果有限）或第三方解包裹算法 % 一个简单的2D解包裹思路（质量图引导）： % phi_wrapped = angle(U_est); % 包裹相位 % 可以使用基于最小二乘的解包裹算法，如 `unwrapPhaseMesh` (需要自己实现或找工具箱) % 或者，直接显示相位梯度（如x和y方向的差分），梯度图不受2π跳变影响，能反映相位的局部变化。 dx_phase = diff(phi_recovered, 1, 2); % x方向差分 dy_phase = diff(phi_recovered, 1, 1); % y方向差分

对于课程作业，如果真实相位本身变化范围不超过2π，也可以直接比较包裹后的相位。但务必在报告中说明这一点。

误差分析图：
- 迭代误差曲线：semilogy(1:iter, error_history(1:iter))。使用对数坐标可以更清楚地看到误差前期的快速下降和后期的缓慢收敛。
- 相位残差图：imagesc(wrapToPi(phi_recovered - phi_true))。显示恢复相位与真实相位的差值（包裹到[-π, π]）。理想的恢复结果应该是一张均匀的灰度图（差值接近0）。
- 剖面线对比：在相位图上选取一条水平或垂直线，绘制真实相位和恢复相位沿该线的变化曲线，进行直接对比。
```
line_index = N/2; % 中间一行 figure; plot(x, phi_true(line_index, :), 'b-', 'LineWidth', 1.5, 'DisplayName', '真实相位'); hold on; plot(x, phi_recovered(line_index, :), 'r--', 'LineWidth', 1.5, 'DisplayName', '恢复相位'); xlabel('位置 (m)'); ylabel('相位 (rad)'); legend; title('相位剖面线对比');
```

定量指标计算：

% 1. 均方根误差 (RMSE) - 需解包裹或处理跳变 phase_diff = phi_recovered - phi_true; % 处理2π跳变：将差值映射到[-π, π]区间 phase_diff_wrapped = mod(phase_diff + pi, 2*pi) - pi; rmse = sqrt(mean(phase_diff_wrapped(:).^2)); fprintf('相位RMSE: %.4f rad\n', rmse); % 2. 相关系数 (Correlation Coefficient) C = corrcoef(phi_true(:), phi_recovered(:)); correlation = C(1,2); fprintf('相位相关系数: %.4f\n', correlation); % 3. 对于振幅恢复，可以计算相对误差 amp_recovered = abs(U_est); amp_error = norm(amp_recovered(:) - A_obj(:)) / norm(A_obj(:)); fprintf('振幅相对误差: %.4f\n', amp_error);

5.2 探索性仿真实验设计

为了深入理解算法，不要只满足于一个成功运行的例子。设计一系列对比实验，能让你的报告脱颖而出：

噪声水平影响实验：固定其他参数，改变添加到I_measured上的高斯噪声方差（例如从0到0.1），观察恢复相位的RMSE和迭代收敛曲线如何变化。绘制“噪声水平-恢复误差”关系图。
初始猜测敏感性实验：对同一个物体和噪声水平，分别用随机相位、零相位和一种有偏的初始猜测（例如，使用一个离焦的球面波相位）运行GS算法。比较它们的最终误差、收敛速度和恢复结果的主观质量。
支持域约束实验：如果物体是一个简单形状（如矩形），尝试在物面约束中引入支持域。比较使用精确振幅约束、宽松支持域约束（HIO算法）以及无支持域信息时，算法的恢复能力。可以尝试恢复一个在支持域外也有微弱信号的物体。
欠采样与混叠实验：故意减少采样点数N，或者增加相位物体的空间频率，使其超过奈奎斯特频率。观察混叠如何导致恢复失败，并解释原因。

5.3 项目报告与代码组织心得

最后，将你的工作整理成一份可交付的作业或项目报告。

代码组织建议：

使用脚本（.m脚本）作为主程序，调用多个函数文件（.m函数）。
将核心功能模块化：create_test_object.m,forward_propagation.m,gs_algorithm.m,hio_algorithm.m,plot_results.m。
使用清晰的变量名，添加充分的注释。特别是对于光学参数和坐标定义，注释其物理意义和单位。
在主脚本开头，集中定义所有可调参数，方便进行不同的实验。
使用savefig或exportgraphics将生成的图表高质量保存，便于插入报告。

报告撰写要点：

引言：简要说明相位恢复的意义、GS算法的背景和本项目目标。
原理：清晰阐述GS算法的数学原理和迭代步骤，最好配以流程图。
方法：详细描述你的MATLAB仿真模型，包括光学设置、参数、噪声模型以及你实现的GS算法细节（是否包含HIO？如何处理支持域？）。
结果与分析：这是核心部分。展示关键实验结果图（原始物体、测量光强、恢复相位、误差曲线、剖面对比等）。结合图表，分析算法性能，讨论前面提到的关键因素（初始猜测、噪声、支持域）的影响。展示你设计的探索性实验的结果和结论。
讨论与结论：总结GS算法的优缺点，你在实现过程中遇到的主要挑战和解决方案，以及对未来改进的设想（如尝试更先进的算法：Fienup算法、PIE、ePIE等）。
参考文献：引用关键的教材、论文或网络资源。

通过这样一个从理论到实践，从实现到分析，从基础到拓展的完整过程，你完成的就不仅仅是一个课程作业，而是一个扎实的光学计算研究小课题。这份经验对你理解更复杂的计算成像算法，乃至从事相关领域的研究或开发工作，都将是一块重要的基石。

本文还有配套的精品资源，点击获取