MATLAB实战：手把手教你用GS和TIE算法恢复丢失的图像相位（附完整代码）

MATLAB实战：从衍射强度图到相位恢复的完整实现指南

在光学成像和数字图像处理领域，相位恢复是一个经典而富有挑战性的问题。当我们只能获取物体的衍射强度图而丢失了相位信息时，如何准确重建原始图像？本文将深入探讨两种核心算法——Gerchberg-Saxton（GS）算法和传输强度方程（TIE）的MATLAB实现，带您一步步完成从理论到代码的完整过程。

1. 相位恢复基础与环境准备

相位恢复问题的本质可以概括为：已知光波场的强度分布，如何反推出其相位分布。这个问题在X射线晶体学、电子显微镜和光学成像等领域有着广泛的应用。我们先来搭建MATLAB实验环境。

1.1 必要工具与数据准备

首先确保您的MATLAB安装了Image Processing Toolbox。我们将使用以下测试图像：

% 下载测试图像 img_url = 'https://example.com/test_image.tif'; % 替换为实际图像URL websave('test_image.tif', img_url); % 基础参数设置 lambda = 632.8e-6; % 波长(mm) d = 20; % 衍射距离(mm) pixel_size = 8e-3; % 像素大小(mm)

1.2 角谱传播理论实现

角谱传播是理解相位恢复的关键理论基础。以下是其MATLAB实现：

function [Uz, I] = angular_spectrum_propagation(U0, lambda, d, L) [N, ~] = size(U0); fx = ([0:fix(N/2), ceil(N/2)-1:-1:1])/L; [fx, fy] = meshgrid(fx, fx); H = exp(1j*2*pi*d*sqrt(1/lambda^2 - (fx.^2 + fy.^2))); Uz_fft = fft2(U0) .* H; Uz = ifft2(Uz_fft); I = abs(Uz).^2; I = I/max(I(:)); % 归一化 end

2. GS算法实现与优化

Gerchberg-Saxton算法是最经典的相位恢复算法，其核心思想是通过在空域和频域之间反复迭代，逐步逼近真实解。

2.1 基础GS算法实现

function [phase, loss_history] = basic_GS(I_measured, iterations) [N, ~] = size(I_measured); phase = 2*pi*rand(N); % 随机初始相位 loss_history = zeros(iterations, 1); for k = 1:iterations % 正向传播 E = sqrt(I_measured) .* exp(1j*phase); E_f = fft2(E); % 频域约束 E_f_constrained = abs(E_f) .* exp(1j*angle(E_f)); % 反向传播 E_new = ifft2(E_f_constrained); phase = angle(E_new); % 记录损失 loss_history(k) = sum(abs(abs(E_new).^2 - I_measured), 'all')/N^2; end end

2.2 GS算法的常见问题与解决方案

在实际应用中，我们发现基础GS算法存在几个典型问题：

1. 收敛速度慢：迭代数百次才能得到较好结果
2. 陷入局部最优：特别是对于复杂图像
3. 对噪声敏感：实测数据中的噪声会显著影响结果

改进策略：

• 引入松弛因子加速收敛：

  alpha = 0.8; % 松弛因子 phase = alpha*phase_new + (1-alpha)*phase_old;

• 混合输入输出（HIO）方法：

  beta = 0.9; % HIO参数 E_support = (abs(E_new) > threshold); E_corrected = E_support .* E_new + (1-E_support) .* (E_prev - beta*E_new);

3. TIE算法实现与技巧

传输强度方程（Transport of Intensity Equation）提供了另一种相位恢复的思路，特别适合处理聚焦和离焦图像对。

3.1 TIE核心算法实现

function phase = solve_TIE(I1, I2, I3, dz, lambda, L) % I1, I2, I3: 分别对应-dz, 0, +dz位置的强度图 D = (I3 - I1)/(2*dz); [N, ~] = size(I2); % 频率坐标 fx = ([0:fix(N/2), ceil(N/2)-1:-1:1])/L; [fx, fy] = meshgrid(fx, fx); q = fx.^2 + fy.^2; % 求解TIE numerator = fft2(D./I2); denominator = -4*pi^2*(q + eps); % 避免除以零 phase = real(ifft2(numerator./denominator))/lambda; phase = phase - min(phase(:)); % 平移相位到正值 end

3.2 TIE实践中的关键参数

注意：TIE求解需要良好的图像对齐。建议使用以下方法进行亚像素对齐：

[output, Greg] = dftregistration(fft2(I1), fft2(I2), 100);

4. 混合算法与性能提升

结合GS和TIE的优势，我们可以构建更强大的混合算法。以下是两种典型方案：

4.1 TIE初始化+GS精修

% 先用TIE获取初始相位估计 phase_TIE = solve_TIE(I1, I2, I3, dz, lambda, L); % 然后用GS算法进行精修 [phase_final, loss] = improved_GS(I2, phase_TIE, 200); function [phase, loss] = improved_GS(I_target, initial_phase, iterations) phase = initial_phase; for k = 1:iterations % 加入动量项的GS更新 if k > 1 momentum = 0.9*(phase - phase_prev); else momentum = 0; end E = sqrt(I_target) .* exp(1j*phase); E_f = fft2(E); E_f_constrained = abs(E_f) .* exp(1j*angle(E_f)); E_new = ifft2(E_f_constrained); phase_prev = phase; phase = angle(E_new) + momentum; loss(k) = calculate_loss(E_new, I_target); end end

4.2 多距离数据融合方法

当有多个传播距离的强度图时，可以构建更稳健的恢复流程：

1. 数据预处理流程：

   • 对所有强度图进行归一化
   • 使用亚像素方法精确对齐
   • 应用非局部均值去噪

2. 分层恢复策略：

   % 先用大距离数据恢复低频成分 phase_low = solve_TIE(I_d30, I_d20, I_d10, 10, lambda, L); % 再用小距离数据恢复高频细节 phase_high = basic_GS(I_d10, 50); % 频域融合 phase_low_f = fft2(phase_low); phase_high_f = fft2(phase_high); cutoff = 0.3; % 截止频率 [fx, fy] = meshgrid(linspace(-0.5,0.5,N)); mask = (sqrt(fx.^2 + fy.^2) < cutoff); phase_fused_f = phase_low_f.*mask + phase_high_f.*(1-mask); phase_final = real(ifft2(phase_fused_f));

5. 结果评估与可视化

良好的评估体系可以帮助我们判断算法效果并指导参数调整。

5.1 定量评估指标

function [mse, psnr, ssim] = evaluate_results(phase_est, phase_gt) % 均方误差 mse = immse(phase_est, phase_gt); % 峰值信噪比 max_val = max(phase_gt(:)); psnr = 10*log10(max_val^2/mse); % 结构相似性 ssim = ssim_index(phase_est, phase_gt); end

5.2 可视化技巧

多图对比显示：

figure; subplot(2,2,1); imshow(phase_gt, []); title('真实相位'); subplot(2,2,2); imshow(phase_GS, []); title('GS恢复结果'); subplot(2,2,3); imshow(phase_TIE, []); title('TIE恢复结果'); subplot(2,2,4); imshow(phase_hybrid, []); title('混合算法结果'); % 添加彩色条统一刻度 for i = 1:4 subplot(2,2,i); caxis([0 2*pi]); colorbar; end

迭代过程动画：

figure; h = imshow(phase_init, []); for k = 1:iterations % 迭代更新phase... set(h, 'CData', phase_current); title(['迭代次数: ' num2str(k)]); drawnow; pause(0.05); end

6. 实战经验与调试技巧

在实际项目中，我们积累了一些宝贵经验：

1. 数据归一化的艺术：

   • 强度图建议归一化到[0,1]范围
   • 相位图保持在[0,2π]或[-π,π]范围
   • 迭代过程中定期进行归一化防止数值溢出

2. 迭代停止条件：

   tolerance = 1e-6; for k = 1:max_iter % 迭代计算... if k > 1 && abs(loss(k)-loss(k-1)) < tolerance break; end end

3. 并行计算加速：

   parfor d_idx = 1:num_distances results{d_idx} = process_single_distance(data{d_idx}); end

4. 常见问题排查表：

在完成一个实际生物样本的相位恢复项目时，我们发现将TIE的初始估计与GS的精修相结合，配合适当的多距离数据融合策略，能够在保持细节的同时获得更稳定的恢复结果。特别是在处理低信噪比数据时，加入TV正则项的改进GS算法表现尤为突出。