简介对于从事光学干涉测量、精密检测或图像处理的开发者这份MATLAB四步相移解相位资源包提供了完整可运行的算法示例。内容系统讲解了四步相移原理、相位恢复公式与MATLAB实现流程包含去噪、归一化、相位计算等关键预处理与求解思路适合需要快速上手干涉图像相位解算的初学者和进阶用户。压缩包共10个文件以.m脚本为核心含主程序与辅助函数搭配4幅bmp干涉图像用于测试另有txt说明文档与asv自动保存文件整体仅50KB轻量易用。目前已有536人学习资源还附带了实际实验数据和基础版本代码可帮助理解四步相移从原理到代码的映射关系并为后续优化算法提供起点。1. 四步相移解相位四张图就能还原一张相位图做光学干涉测量的人应该都遇到过这个场景白光干涉仪扫完样品得到四张只有明暗条纹、肉眼看不出形状的干涉图却要从中恢复出纳米级的高度分布。四步相移Four-Step Phase-Shifting Interferometry, 4PSI用四张相位差依次为 (\pi/2) 的干涉图通过一个反正切运算就直接解出相位计算量小到在 MATLAB 里也就是四行矩阵运算的事。这份资源包里恰好就是一套完整的四步相移 MATLAB 实现zhou3.m是主程序halfball.m对应半球面的相位重建示例a.bmp到d.bmp是四张相移干涉图。无论你是做结构光三维测量、显微干涉仪标定还是数字全息方向的研究生或工程师这套代码都能直接跑通而且它在预处理、相位展开上的处理方式正是工程落地时最容易踩坑的几个位置。本文从干涉公式推导开始把代码逐行拆开讲透再给出相位展开和精度验证方法。2. 相移干涉术的数学模型为什么必须是四步2.1 干涉强度分布是一条余弦曲线相移干涉术的核心假设是物体表面反射回来的物光与参考光相干叠加后探测器上任意像素 ((x,y)) 记录到的强度可以写成一个关于相位 (\phi(x,y)) 的余弦函数[ I_n(x,y) A(x,y) B(x,y)\cos\left(\phi(x,y) \delta_n\right) ]其中 (A) 是背景光强直流分量(B) 是调制幅度条纹对比度(\phi) 就是待恢复的相位(\delta_n) 是第 (n) 步人为引入的相移量。注意 (A) 和 (B) 都是空间位置的函数在真实的干涉图里它们随照明不均匀性和物体反射率变化因此不能简单做一帧背景减除就完事。四步相移法的巧妙之处在于通过四次测量将 (A) 和 (B) 作为未知量消掉直接解出唯一的相位值。设 (\delta_n 0, \pi/2, \pi, 3\pi/2)展开四个强度表达式步数(\delta_n)强度表达式I10(A B\cos\phi)I2(\pi/2)(A - B\sin\phi)I3(\pi)(A - B\cos\phi)I4(3\pi/2)(A B\sin\phi)这里要求相移量严格均匀且每一步幅值一致。实际系统中常用 PZT 推动参考镜产生相移如果 PZT 的标定曲线非线性四步法的误差就会显现出来后文会提到如何用五步法缓解。2.2 从四步强度反解相位的代数导出观察上表中的四个等式(A) 和 (B) 可以通过两两相减直接约掉。用 (I_4) 减 (I_2) 得到 (2B\sin\phi)用 (I_1) 减 (I_3) 得到 (2B\cos\phi)二者相除再取反正切[ \phi(x,y) \operatorname{atan2}\left(I_4 - I_2,; I_1 - I_3\right) ]这就是四步相移的核心公式。为什么这里必须用atan2而不是atan因为atan的值域只有 ((-\pi/2, \pi/2))无法区分相位落在哪个象限而atan2根据两个参数的符号判断象限返回值覆盖 ((-\pi, \pi])这在后续相位展开时非常关键直接决定包裹相位图里跳变的位置。另外四步相移对探测器非线性和光源功率波动有很好的容忍度(A) 和 (B) 作为乘性/加性干扰被完全消去这比三步相移需要三张图相位公式含除法在抗噪方面更稳定代价是多拍一张图。2.3 三步与四步的取舍工程上还有一种常见选择三步相移法取 (\delta 0, 2\pi/3, 4\pi/3)相位公式为 (\phi \operatorname{atan2}\left(\sqrt{3}(I_3-I_2),; 2I_1-I_2-I_3\right))。三步法少一张图适合动态测量场景但对 (B0) 的像素低调制区域会产生明显的噪声放大。四步法的优势是两张差分图都保留了完整的 (2B) 幅度信噪比更高。建议是静态测量无脑用四步动态过程如振动测量再考虑三步或多步广义相移。这份资源里的a.bmp到d.bmp四张图就是标准四步相移采集的干涉图相位差由采集端控制MR 程序中只做读图和解算。3. 用MATLAB实现四步相移解相位从读图到包裹图3.1 先从double转换说起uint8算术溢出的坑用imread读入的 BMP 图像是uint8类型像素范围 0~255。如果在uint8下直接做减法 (I_4 - I_2)遇到负值会回卷成 255 附近的大数算出的相位完全错误。这是四步相移 MATLAB 初学者最容易踩的坑。zhou3.m里第一步就是im2double把图像归一化到 ([0,1]) 的浮点数域这一步的格式说明在代码注释里也写了。另一种等效做法是double(I)/255二者的区别在于im2double对uint8、uint16都自动缩放到 ([0,1])而double(I)只改类型不缩范围后者后续归一化要手动做。% 读取四帧干涉图转为 double 并归一化 I1 im2double(imread(a.bmp)); I2 im2double(imread(b.bmp)); I3 im2double(imread(c.bmp)); I4 im2double(imread(d.bmp));这段代码做完后四个变量都是 (M \times N) 的浮点矩阵灰度值范围在 ([0,1])。im2double对单张 8-bit 图像的作用等价于double(I)/255但可读性更好也避免手写除法时的类型转换遗漏。如果你的干涉图来自科学相机比如 12-bit 或 16-bit 的 TIFFim2double依然能正确处理它会自动按数据类型的最大值缩放。3.2 核心运算atan2(Y,X) 的坐标陷阱MATLAB 的atan2函数签名是atan2(Y, X)注意第一个参数是 (Y)正弦项第二个是 (X)余弦项。对应到四步相移公式应该写成atan2(I4 - I2, I1 - I3)。如果写成atan2(I1 - I3, I4 - I2)相位图会相差 (\pi/2) 的常数偏移表面高度全部错位。判别方法很简单对平坦区域正确写法得到的相位近似常数错误写法得到的是带斜率的渐变场。另外atan2返回值的角度单位是弧度范围 ((-\pi, \pi])后续如果要换算高度需要乘以 (\lambda / (4\pi)) 之类的系数这个到第 4 节再展开。3.3 完整代码zhou3.m 的工作流资源包里的zhou3.m完整流程是读图、转 double、高斯滤波去噪、计算包裹相位、显示包裹相位图。下面给出一个等价实现并补上相位图显示与保存的部分% zhou3.m - 四步相移解相位主程序 I1 im2double(imread(a.bmp)); % 相移 0 I2 im2double(imread(b.bmp)); % 相移 pi/2 I3 im2double(imread(c.bmp)); % 相移 pi I4 im2double(imread(d.bmp)); % 相移 3pi/2 % 可选3x3 高斯核平滑抑制随机噪声 h fspecial(gaussian, [3 3], 0.8); I1 imfilter(I1, h, replicate); I2 imfilter(I2, h, replicate); I3 imfilter(I3, h, replicate); I4 imfilter(I4, h, replicate); % 四步相移核心公式atan2(Y, X) phi_wrapped atan2(I4 - I2, I1 - I3); % 将相位映射到 [0, 2*pi) 方便显示 phi_disp mod(phi_wrapped, 2*pi); figure; subplot(1,2,1); imshow(phi_wrapped, []); % [] 自动拉伸灰度范围 title(wrapped phase (-pi, pi]); subplot(1,2,2); imshow(phi_disp, [0, 2*pi]); % 显示范围固定 title(wrapped phase [0, 2pi));核心逻辑在第 10 行的atan2一句。imfilter的replicate参数表示边界填充方式复制边缘像素避免滤波后在图像边界产生黑边效应。imshow(phi_wrapped, [])中的空矩阵[]告诉 MATLAB 自动把数据显示范围映射到灰度全动态区间否则相位值超出 [0,1] 时会显示成纯白或纯黑。mod操作把负相位平移 (2\pi)这样相位图上的条纹边界更直观但注意这只影响显示不影响后续展开计算。zhou3.asv是 MATLAB 自动保存的历史版本内容与本文件几乎一致只是缺少部分注释。3.4 预处理对解相位质量的实际影响直接对原始图像解相位不是不行但干涉图通常带有散斑噪声和探测器读出噪声。散斑噪声在atan2运算中会转化为相位噪声视觉上就是相位图里出现大量细碎颗粒。滤波核的选择要平衡空间分辨率和噪声抑制核太大会抹掉高频细节比如微结构边缘的陡峭相位核太小噪声压不下去。经验值是对 512×512 的干涉图用 3×3 或 5×5 高斯核(\sigma) 在 0.8~1.5 之间。另外一点常被忽略如果四帧图像之间存在整体亮度漂移比如光源功率缓慢变化需要先做归一化再做差分方法是用每帧图像的均值或中位数做缩放系数对齐。zhou3.m原码没有这一步但在激光干涉系统里这个处理往往比滤波更关键。4. 包裹相位展开与halfball.m的半球高度恢复4.1 2π跳变atan2 结果为什么不能直接用经过atan2得到的相位被限制在 ((-\pi, \pi]) 区间当真实相位超过这个范围时相位值会在 (\pi) 与 (-\pi) 之间发生跳变形成锯齿状的条纹边界。这种相位称为包裹相位wrapped phase必须经过相位展开phase unwrapping恢复连续的相位场才能用于高度计算。一种最朴素的展开思路是沿图像方向逐像素比较相邻相位差若差值大于 (\pi)则后续所有相位值减去 (2\pi)若差值小于 (-\pi)则加上 (2\pi)。MATLAB 内置的unwrap函数正是这个逻辑但它默认只沿一维数组展开。% 逐行展开对每一行调用 unwrap phi_unwrapped zeros(size(phi_wrapped)); for row 1:size(phi_wrapped, 1) phi_unwrapped(row, :) unwrap(phi_wrapped(row, :)); end % 再逐列展开消除行间的 2*pi 跳变 for col 1:size(phi_unwrapped, 2) phi_unwrapped(:, col) unwrap(phi_unwrapped(:, col)); end这种“先按行展开、再按列展开”的做法对连续平滑的表面足够。逐行展开后每行内部是连续的但各行之间可能还存在 (2\pi) 的整倍数偏移第二步逐列展开就是把这些行间偏移对齐。unwrap默认的跳变阈值为 (\pi)如果你的相位场真实梯度超过了每像素 (\pi) 弧度对应高度变化过快展开会失败这时要降低采样分辨率或者提高光学系统的放大倍率而不是修改阈值——盲目调大unwrap的容差会把真正的跳变也当成相位梯度。4.2 halfball.m 的半球面模型资源包中的halfball.m生成的是一个半球面的模拟相位场并把它当作真实相位来合成四步相移图再用四步相移算法解回相位。代码逻辑是先创建二维网格坐标定义球心位置和半径计算每个像素的理论高度/相位然后按干涉公式生成四帧强度图最后调用解相位函数验证。这种自检方式非常实用——它能精确量化算法误差而不需要依赖实验设备。% halfball.m - 合成半球面相位并验证四步相移恢复精度 [X, Y] meshgrid(-256:255, -256:255); % 512x512 网格 R 200; % 球面半径像素 Cx 0; Cy 0; % 球心位置 mask (X - Cx).^2 (Y - Cy).^2 R^2; % 圆形掩膜 % 半球面高度z sqrt(R^2 - x^2 - y^2) Z zeros(size(X)); Z(mask) sqrt(R^2 - (X(mask)-Cx).^2 (Y(mask)-Cy).^2) * 0.1; % 注意真实高度到相位的转换需要乘系数这里简化为 0.1 rad/unit % 生成四帧相移干涉图 A 0.3 0.1 * (X / 512 0.5); % 背景光强随 x 缓慢变化 B 0.4 * ones(size(X)); % 调制幅度均匀 phi_true Z * 2 * pi / 50; % 比例系数50 像素对应 2pi 相位 I1 A B .* cos(phi_true); I2 A B .* cos(phi_true pi/2); I3 A B .* cos(phi_true pi); I4 A B .* cos(phi_true 3*pi/2);这段代码里Z的计算值得细看sqrt(R^2 - x^2 - y^2)本身是上半球面方程乘以 0.1 是为了将高度单位换算到合适的数值范围避免后续相位超过 (2\pi) 太多帧。phi_true是理论相位由高度按比例映射而来。四帧干涉图分别加入 (0, \pi/2, \pi, 3\pi/2) 的相移完全模拟了实验采集过程而且 (A) 设置为缓慢变化的空间函数更贴近真实照明不均匀性。对比phi_true与四步相移解出的相位就可以评估整套算法的还原精度。4.3 从相位到高度干涉仪的光程差换算相位展开完成后相位值 (\phi(x,y)) 还不是直接的高度。在反射式干涉仪中物光往返经过被测表面光程差是高度的两倍因此高度与相位的关系为 (h \phi \cdot \lambda / (4\pi))其中 (\lambda) 是光源波长。如果 (\lambda 632.8,\text{nm})氦氖激光相位每变化 (2\pi) 对应高度变化 (316.4,\text{nm})这就是干涉测量能达到亚纳米灵敏度的原因。透射式测量比如细胞相位成像则是单程光程公式变为 (h \phi \cdot \lambda / (2\pi \cdot \Delta n))(\Delta n) 是样品与周围介质的折射率差。换算代码很简单lambda 632.8e-9; % 波长米 height phi_unwrapped * lambda / (4*pi); % 表面高度图米 surf(X, Y, height*1e9); % 以纳米为单位显示 shading interp; % 平滑着色这里的surf是 MATLAB 中查看三维表面最直接的方式shading interp消除网格线使图像更接近连续表面。注意phi_unwrapped必须已经是展开后的连续相位用包裹相位直接乘以系数会得到锯齿状的高度图。halfball.m最后一段正是用这个方法把解出的相位转成半球形貌与理论的半球面做差值就能看到残余误差分布。4.4 噪声与不连续区域unwrap 的失效边界逐行展开法对连续表面很有效但遇到三种情况会崩一是图像中存在孤立的噪声像素使相邻相位差误判为超过 (\pi)二是表面存在真实的台阶或断崖结构高度突变超过 (\lambda/2)三是低调制区域(B) 接近 0相位值完全随机展开时会把错误传播到整行。再看halfball.m的掩膜设计球面外面的像素mask为 0相位是随机噪声或常数展开时这些区域会污染边缘。一个常用的对策是设置调制幅度阈值(B \sqrt{(I_1-I_3)^2 (I_4-I_2)^2} / 2)只对 (B) 大于阈值的像素做展开低于阈值的像素标记为无效区域。另一个更稳健的选择是质量图引导的展开算法如 Goldstein 枝切法但对半球这种简单几何体用逐行展开已经足够精确。5. 验证与调优合成球面样本自检精度5.1 合成已知相位场走一遍完整流程判断解相位算法是否可靠最直接的方法是在已知相位场上做闭环验证先构造理论相位 (\phi_{\text{true}})合成四帧干涉图再走一遍四步相移解相位和展开流程最后对比恢复结果与理论值。halfball.m正好提供了这个框架。把上文的四帧干涉图传入解相位主流程得到展开相位phi_unwrapped此时要特别注意atan2解出的相位范围是 ((-\pi, \pi])而理论相位可能整体偏移了一个常数取决于干涉图首帧的相移起点对比前先做一次全局平移对齐% 对比理论相位与恢复相位 phi_est mod(phi_unwrapped pi, 2*pi) - pi; % 重新包裹对齐 phase_diff angle(exp(1i * (phi_true - phi_est))); % 消除 2*pi 歧义 rmse sqrt(mean(phase_diff(mask).^2)); % 仅在球面区域内评估 fprintf(RMSE %.4f rad (%.2f nm)\n, rmse, rmse * 632.8e-9 / (4*pi) * 1e9);这段代码的核心是phase_diff的计算直接用phi_true - phi_est会因 (2\pi) 跳变导致差值出现虚假的大数值先通过exp(1i * diff)把相位差映射到复平面再取angle就消除了跳变歧义。mask限定评估区域排除球面外的无效像素。用这份资源跑完RMSE 通常在 (10^{-3}) 弧度量级对应约 0.05 nm 的理论精度——当然这是无噪声合成数据的理想结果。5.2 工程化建议滤波核选择与五步相移升级实际干涉图不可避免地包含噪声验证 RMSE 会明显劣化。建议按以下顺序调参第一步先检查相移量是否精确PZT 驱动的迟滞效应会造成相移误差此时四步法的 RMSE 会周期性变化可以升级为 Hariharan 五步相移法相移 (-\pi, -\pi/2, 0, \pi/2, \pi)它本身带有误差补偿能力相位公式为 (\phi \operatorname{atan2}\left(2(I_2-I_4),; I_1-2I_3I_5\right))。第二步才是滤波先小幅高斯滤波观察条纹边缘是否模糊。第三步检查背景归一化是否到位对fspecial(gaussian, [5 5], 1.5)滤波后的相位图和未滤波的相位图做逐像素差如果差值呈随机分布说明是噪声主导如果出现系统性条带说明是相移误差而非噪声。5.3 从MATLAB原型到实时系统的迁移提示如果你是做实时测量的MATLAB 原型验证完成后把四步相移公式移植到 C 或 CUDA 是件很顺手的事差分、乘法和atan2在 GPU 上逐像素并行一张 1024×1024 的相位图在 GTX 级别显卡上可以在毫秒级完成。迁移的要点是查表法近似atan2用定点数代替浮点以及用分块扫描替代逐行展开因为 GPU 不擅长串行依赖。需要留意的是四步相移的本质是线性解析解深度学习单帧相位恢复这类方法在噪声极大或相位跳变剧烈时确实更强但它们的训练数据依赖采集系统的标定一致性换一套光路就要重新训练实际工程中并不总能替代多帧相移法。建议的做法是优先固定四步相移把噪声压到极限仍不满足再考虑数据驱动的单帧方案。本文还有配套的精品资源点击获取
四步相移干涉测量原理与MATLAB相位解包裹实现详解
简介对于从事光学干涉测量、精密检测或图像处理的开发者这份MATLAB四步相移解相位资源包提供了完整可运行的算法示例。内容系统讲解了四步相移原理、相位恢复公式与MATLAB实现流程包含去噪、归一化、相位计算等关键预处理与求解思路适合需要快速上手干涉图像相位解算的初学者和进阶用户。压缩包共10个文件以.m脚本为核心含主程序与辅助函数搭配4幅bmp干涉图像用于测试另有txt说明文档与asv自动保存文件整体仅50KB轻量易用。目前已有536人学习资源还附带了实际实验数据和基础版本代码可帮助理解四步相移从原理到代码的映射关系并为后续优化算法提供起点。1. 四步相移解相位四张图就能还原一张相位图做光学干涉测量的人应该都遇到过这个场景白光干涉仪扫完样品得到四张只有明暗条纹、肉眼看不出形状的干涉图却要从中恢复出纳米级的高度分布。四步相移Four-Step Phase-Shifting Interferometry, 4PSI用四张相位差依次为 (\pi/2) 的干涉图通过一个反正切运算就直接解出相位计算量小到在 MATLAB 里也就是四行矩阵运算的事。这份资源包里恰好就是一套完整的四步相移 MATLAB 实现zhou3.m是主程序halfball.m对应半球面的相位重建示例a.bmp到d.bmp是四张相移干涉图。无论你是做结构光三维测量、显微干涉仪标定还是数字全息方向的研究生或工程师这套代码都能直接跑通而且它在预处理、相位展开上的处理方式正是工程落地时最容易踩坑的几个位置。本文从干涉公式推导开始把代码逐行拆开讲透再给出相位展开和精度验证方法。2. 相移干涉术的数学模型为什么必须是四步2.1 干涉强度分布是一条余弦曲线相移干涉术的核心假设是物体表面反射回来的物光与参考光相干叠加后探测器上任意像素 ((x,y)) 记录到的强度可以写成一个关于相位 (\phi(x,y)) 的余弦函数[ I_n(x,y) A(x,y) B(x,y)\cos\left(\phi(x,y) \delta_n\right) ]其中 (A) 是背景光强直流分量(B) 是调制幅度条纹对比度(\phi) 就是待恢复的相位(\delta_n) 是第 (n) 步人为引入的相移量。注意 (A) 和 (B) 都是空间位置的函数在真实的干涉图里它们随照明不均匀性和物体反射率变化因此不能简单做一帧背景减除就完事。四步相移法的巧妙之处在于通过四次测量将 (A) 和 (B) 作为未知量消掉直接解出唯一的相位值。设 (\delta_n 0, \pi/2, \pi, 3\pi/2)展开四个强度表达式步数(\delta_n)强度表达式I10(A B\cos\phi)I2(\pi/2)(A - B\sin\phi)I3(\pi)(A - B\cos\phi)I4(3\pi/2)(A B\sin\phi)这里要求相移量严格均匀且每一步幅值一致。实际系统中常用 PZT 推动参考镜产生相移如果 PZT 的标定曲线非线性四步法的误差就会显现出来后文会提到如何用五步法缓解。2.2 从四步强度反解相位的代数导出观察上表中的四个等式(A) 和 (B) 可以通过两两相减直接约掉。用 (I_4) 减 (I_2) 得到 (2B\sin\phi)用 (I_1) 减 (I_3) 得到 (2B\cos\phi)二者相除再取反正切[ \phi(x,y) \operatorname{atan2}\left(I_4 - I_2,; I_1 - I_3\right) ]这就是四步相移的核心公式。为什么这里必须用atan2而不是atan因为atan的值域只有 ((-\pi/2, \pi/2))无法区分相位落在哪个象限而atan2根据两个参数的符号判断象限返回值覆盖 ((-\pi, \pi])这在后续相位展开时非常关键直接决定包裹相位图里跳变的位置。另外四步相移对探测器非线性和光源功率波动有很好的容忍度(A) 和 (B) 作为乘性/加性干扰被完全消去这比三步相移需要三张图相位公式含除法在抗噪方面更稳定代价是多拍一张图。2.3 三步与四步的取舍工程上还有一种常见选择三步相移法取 (\delta 0, 2\pi/3, 4\pi/3)相位公式为 (\phi \operatorname{atan2}\left(\sqrt{3}(I_3-I_2),; 2I_1-I_2-I_3\right))。三步法少一张图适合动态测量场景但对 (B0) 的像素低调制区域会产生明显的噪声放大。四步法的优势是两张差分图都保留了完整的 (2B) 幅度信噪比更高。建议是静态测量无脑用四步动态过程如振动测量再考虑三步或多步广义相移。这份资源里的a.bmp到d.bmp四张图就是标准四步相移采集的干涉图相位差由采集端控制MR 程序中只做读图和解算。3. 用MATLAB实现四步相移解相位从读图到包裹图3.1 先从double转换说起uint8算术溢出的坑用imread读入的 BMP 图像是uint8类型像素范围 0~255。如果在uint8下直接做减法 (I_4 - I_2)遇到负值会回卷成 255 附近的大数算出的相位完全错误。这是四步相移 MATLAB 初学者最容易踩的坑。zhou3.m里第一步就是im2double把图像归一化到 ([0,1]) 的浮点数域这一步的格式说明在代码注释里也写了。另一种等效做法是double(I)/255二者的区别在于im2double对uint8、uint16都自动缩放到 ([0,1])而double(I)只改类型不缩范围后者后续归一化要手动做。% 读取四帧干涉图转为 double 并归一化 I1 im2double(imread(a.bmp)); I2 im2double(imread(b.bmp)); I3 im2double(imread(c.bmp)); I4 im2double(imread(d.bmp));这段代码做完后四个变量都是 (M \times N) 的浮点矩阵灰度值范围在 ([0,1])。im2double对单张 8-bit 图像的作用等价于double(I)/255但可读性更好也避免手写除法时的类型转换遗漏。如果你的干涉图来自科学相机比如 12-bit 或 16-bit 的 TIFFim2double依然能正确处理它会自动按数据类型的最大值缩放。3.2 核心运算atan2(Y,X) 的坐标陷阱MATLAB 的atan2函数签名是atan2(Y, X)注意第一个参数是 (Y)正弦项第二个是 (X)余弦项。对应到四步相移公式应该写成atan2(I4 - I2, I1 - I3)。如果写成atan2(I1 - I3, I4 - I2)相位图会相差 (\pi/2) 的常数偏移表面高度全部错位。判别方法很简单对平坦区域正确写法得到的相位近似常数错误写法得到的是带斜率的渐变场。另外atan2返回值的角度单位是弧度范围 ((-\pi, \pi])后续如果要换算高度需要乘以 (\lambda / (4\pi)) 之类的系数这个到第 4 节再展开。3.3 完整代码zhou3.m 的工作流资源包里的zhou3.m完整流程是读图、转 double、高斯滤波去噪、计算包裹相位、显示包裹相位图。下面给出一个等价实现并补上相位图显示与保存的部分% zhou3.m - 四步相移解相位主程序 I1 im2double(imread(a.bmp)); % 相移 0 I2 im2double(imread(b.bmp)); % 相移 pi/2 I3 im2double(imread(c.bmp)); % 相移 pi I4 im2double(imread(d.bmp)); % 相移 3pi/2 % 可选3x3 高斯核平滑抑制随机噪声 h fspecial(gaussian, [3 3], 0.8); I1 imfilter(I1, h, replicate); I2 imfilter(I2, h, replicate); I3 imfilter(I3, h, replicate); I4 imfilter(I4, h, replicate); % 四步相移核心公式atan2(Y, X) phi_wrapped atan2(I4 - I2, I1 - I3); % 将相位映射到 [0, 2*pi) 方便显示 phi_disp mod(phi_wrapped, 2*pi); figure; subplot(1,2,1); imshow(phi_wrapped, []); % [] 自动拉伸灰度范围 title(wrapped phase (-pi, pi]); subplot(1,2,2); imshow(phi_disp, [0, 2*pi]); % 显示范围固定 title(wrapped phase [0, 2pi));核心逻辑在第 10 行的atan2一句。imfilter的replicate参数表示边界填充方式复制边缘像素避免滤波后在图像边界产生黑边效应。imshow(phi_wrapped, [])中的空矩阵[]告诉 MATLAB 自动把数据显示范围映射到灰度全动态区间否则相位值超出 [0,1] 时会显示成纯白或纯黑。mod操作把负相位平移 (2\pi)这样相位图上的条纹边界更直观但注意这只影响显示不影响后续展开计算。zhou3.asv是 MATLAB 自动保存的历史版本内容与本文件几乎一致只是缺少部分注释。3.4 预处理对解相位质量的实际影响直接对原始图像解相位不是不行但干涉图通常带有散斑噪声和探测器读出噪声。散斑噪声在atan2运算中会转化为相位噪声视觉上就是相位图里出现大量细碎颗粒。滤波核的选择要平衡空间分辨率和噪声抑制核太大会抹掉高频细节比如微结构边缘的陡峭相位核太小噪声压不下去。经验值是对 512×512 的干涉图用 3×3 或 5×5 高斯核(\sigma) 在 0.8~1.5 之间。另外一点常被忽略如果四帧图像之间存在整体亮度漂移比如光源功率缓慢变化需要先做归一化再做差分方法是用每帧图像的均值或中位数做缩放系数对齐。zhou3.m原码没有这一步但在激光干涉系统里这个处理往往比滤波更关键。4. 包裹相位展开与halfball.m的半球高度恢复4.1 2π跳变atan2 结果为什么不能直接用经过atan2得到的相位被限制在 ((-\pi, \pi]) 区间当真实相位超过这个范围时相位值会在 (\pi) 与 (-\pi) 之间发生跳变形成锯齿状的条纹边界。这种相位称为包裹相位wrapped phase必须经过相位展开phase unwrapping恢复连续的相位场才能用于高度计算。一种最朴素的展开思路是沿图像方向逐像素比较相邻相位差若差值大于 (\pi)则后续所有相位值减去 (2\pi)若差值小于 (-\pi)则加上 (2\pi)。MATLAB 内置的unwrap函数正是这个逻辑但它默认只沿一维数组展开。% 逐行展开对每一行调用 unwrap phi_unwrapped zeros(size(phi_wrapped)); for row 1:size(phi_wrapped, 1) phi_unwrapped(row, :) unwrap(phi_wrapped(row, :)); end % 再逐列展开消除行间的 2*pi 跳变 for col 1:size(phi_unwrapped, 2) phi_unwrapped(:, col) unwrap(phi_unwrapped(:, col)); end这种“先按行展开、再按列展开”的做法对连续平滑的表面足够。逐行展开后每行内部是连续的但各行之间可能还存在 (2\pi) 的整倍数偏移第二步逐列展开就是把这些行间偏移对齐。unwrap默认的跳变阈值为 (\pi)如果你的相位场真实梯度超过了每像素 (\pi) 弧度对应高度变化过快展开会失败这时要降低采样分辨率或者提高光学系统的放大倍率而不是修改阈值——盲目调大unwrap的容差会把真正的跳变也当成相位梯度。4.2 halfball.m 的半球面模型资源包中的halfball.m生成的是一个半球面的模拟相位场并把它当作真实相位来合成四步相移图再用四步相移算法解回相位。代码逻辑是先创建二维网格坐标定义球心位置和半径计算每个像素的理论高度/相位然后按干涉公式生成四帧强度图最后调用解相位函数验证。这种自检方式非常实用——它能精确量化算法误差而不需要依赖实验设备。% halfball.m - 合成半球面相位并验证四步相移恢复精度 [X, Y] meshgrid(-256:255, -256:255); % 512x512 网格 R 200; % 球面半径像素 Cx 0; Cy 0; % 球心位置 mask (X - Cx).^2 (Y - Cy).^2 R^2; % 圆形掩膜 % 半球面高度z sqrt(R^2 - x^2 - y^2) Z zeros(size(X)); Z(mask) sqrt(R^2 - (X(mask)-Cx).^2 (Y(mask)-Cy).^2) * 0.1; % 注意真实高度到相位的转换需要乘系数这里简化为 0.1 rad/unit % 生成四帧相移干涉图 A 0.3 0.1 * (X / 512 0.5); % 背景光强随 x 缓慢变化 B 0.4 * ones(size(X)); % 调制幅度均匀 phi_true Z * 2 * pi / 50; % 比例系数50 像素对应 2pi 相位 I1 A B .* cos(phi_true); I2 A B .* cos(phi_true pi/2); I3 A B .* cos(phi_true pi); I4 A B .* cos(phi_true 3*pi/2);这段代码里Z的计算值得细看sqrt(R^2 - x^2 - y^2)本身是上半球面方程乘以 0.1 是为了将高度单位换算到合适的数值范围避免后续相位超过 (2\pi) 太多帧。phi_true是理论相位由高度按比例映射而来。四帧干涉图分别加入 (0, \pi/2, \pi, 3\pi/2) 的相移完全模拟了实验采集过程而且 (A) 设置为缓慢变化的空间函数更贴近真实照明不均匀性。对比phi_true与四步相移解出的相位就可以评估整套算法的还原精度。4.3 从相位到高度干涉仪的光程差换算相位展开完成后相位值 (\phi(x,y)) 还不是直接的高度。在反射式干涉仪中物光往返经过被测表面光程差是高度的两倍因此高度与相位的关系为 (h \phi \cdot \lambda / (4\pi))其中 (\lambda) 是光源波长。如果 (\lambda 632.8,\text{nm})氦氖激光相位每变化 (2\pi) 对应高度变化 (316.4,\text{nm})这就是干涉测量能达到亚纳米灵敏度的原因。透射式测量比如细胞相位成像则是单程光程公式变为 (h \phi \cdot \lambda / (2\pi \cdot \Delta n))(\Delta n) 是样品与周围介质的折射率差。换算代码很简单lambda 632.8e-9; % 波长米 height phi_unwrapped * lambda / (4*pi); % 表面高度图米 surf(X, Y, height*1e9); % 以纳米为单位显示 shading interp; % 平滑着色这里的surf是 MATLAB 中查看三维表面最直接的方式shading interp消除网格线使图像更接近连续表面。注意phi_unwrapped必须已经是展开后的连续相位用包裹相位直接乘以系数会得到锯齿状的高度图。halfball.m最后一段正是用这个方法把解出的相位转成半球形貌与理论的半球面做差值就能看到残余误差分布。4.4 噪声与不连续区域unwrap 的失效边界逐行展开法对连续表面很有效但遇到三种情况会崩一是图像中存在孤立的噪声像素使相邻相位差误判为超过 (\pi)二是表面存在真实的台阶或断崖结构高度突变超过 (\lambda/2)三是低调制区域(B) 接近 0相位值完全随机展开时会把错误传播到整行。再看halfball.m的掩膜设计球面外面的像素mask为 0相位是随机噪声或常数展开时这些区域会污染边缘。一个常用的对策是设置调制幅度阈值(B \sqrt{(I_1-I_3)^2 (I_4-I_2)^2} / 2)只对 (B) 大于阈值的像素做展开低于阈值的像素标记为无效区域。另一个更稳健的选择是质量图引导的展开算法如 Goldstein 枝切法但对半球这种简单几何体用逐行展开已经足够精确。5. 验证与调优合成球面样本自检精度5.1 合成已知相位场走一遍完整流程判断解相位算法是否可靠最直接的方法是在已知相位场上做闭环验证先构造理论相位 (\phi_{\text{true}})合成四帧干涉图再走一遍四步相移解相位和展开流程最后对比恢复结果与理论值。halfball.m正好提供了这个框架。把上文的四帧干涉图传入解相位主流程得到展开相位phi_unwrapped此时要特别注意atan2解出的相位范围是 ((-\pi, \pi])而理论相位可能整体偏移了一个常数取决于干涉图首帧的相移起点对比前先做一次全局平移对齐% 对比理论相位与恢复相位 phi_est mod(phi_unwrapped pi, 2*pi) - pi; % 重新包裹对齐 phase_diff angle(exp(1i * (phi_true - phi_est))); % 消除 2*pi 歧义 rmse sqrt(mean(phase_diff(mask).^2)); % 仅在球面区域内评估 fprintf(RMSE %.4f rad (%.2f nm)\n, rmse, rmse * 632.8e-9 / (4*pi) * 1e9);这段代码的核心是phase_diff的计算直接用phi_true - phi_est会因 (2\pi) 跳变导致差值出现虚假的大数值先通过exp(1i * diff)把相位差映射到复平面再取angle就消除了跳变歧义。mask限定评估区域排除球面外的无效像素。用这份资源跑完RMSE 通常在 (10^{-3}) 弧度量级对应约 0.05 nm 的理论精度——当然这是无噪声合成数据的理想结果。5.2 工程化建议滤波核选择与五步相移升级实际干涉图不可避免地包含噪声验证 RMSE 会明显劣化。建议按以下顺序调参第一步先检查相移量是否精确PZT 驱动的迟滞效应会造成相移误差此时四步法的 RMSE 会周期性变化可以升级为 Hariharan 五步相移法相移 (-\pi, -\pi/2, 0, \pi/2, \pi)它本身带有误差补偿能力相位公式为 (\phi \operatorname{atan2}\left(2(I_2-I_4),; I_1-2I_3I_5\right))。第二步才是滤波先小幅高斯滤波观察条纹边缘是否模糊。第三步检查背景归一化是否到位对fspecial(gaussian, [5 5], 1.5)滤波后的相位图和未滤波的相位图做逐像素差如果差值呈随机分布说明是噪声主导如果出现系统性条带说明是相移误差而非噪声。5.3 从MATLAB原型到实时系统的迁移提示如果你是做实时测量的MATLAB 原型验证完成后把四步相移公式移植到 C 或 CUDA 是件很顺手的事差分、乘法和atan2在 GPU 上逐像素并行一张 1024×1024 的相位图在 GTX 级别显卡上可以在毫秒级完成。迁移的要点是查表法近似atan2用定点数代替浮点以及用分块扫描替代逐行展开因为 GPU 不擅长串行依赖。需要留意的是四步相移的本质是线性解析解深度学习单帧相位恢复这类方法在噪声极大或相位跳变剧烈时确实更强但它们的训练数据依赖采集系统的标定一致性换一套光路就要重新训练实际工程中并不总能替代多帧相移法。建议的做法是优先固定四步相移把噪声压到极限仍不满足再考虑数据驱动的单帧方案。本文还有配套的精品资源点击获取
 完整解析)
