简介：面向光学测量、结构光三维成像与相位重建领域开发者的MATLAB实现资源，聚焦四步相移法提取包裹相位，并结合最小二乘法完成相位解包裹，解决反正切运算引起相位跳变、难以直接还原连续相位的问题。资源已通过运行验证，算法稳定性较好，适合相关方向研究生、工程师进行原理学习、复现实验与算法二次开发。压缩包共7个文件、约526KB，含2个.m脚本、4个bmp测试图像及1个辅助db文件；脚本分别覆盖四步相移相位计算与最小二乘解包裹主流程，bmp图像提供标准条纹样本，便于直接运行、观察处理前后效果。已有1335人学习下载，程序结构紧凑、注释清楚，也可作为结构光投影测量、机器视觉三维感知等实验教学的参考实现。
1. 四步相移与最小二乘相位解包裹：条纹投影三维测量里避不开的两段代码
四步相移法程序和最小二乘法相位解包裹程序，这两个名字放在一起，基本就是在做条纹投影三维测量或干涉计量中“相位到深度”的最后一公里。拿到四张相移条纹图，先通过四步相移算出(−π, π]的包裹相位；但物体形貌对应的是连续相位，中间那一个个2π跳变必须用解包裹算法接起来。最小二乘解包裹就是其中应用最广、对噪声容忍度最好的一类。这篇文章想帮三类人：正在搭结构光3D测量系统的工程师、做显微干涉或全息重建的研究生、以及要维护光学计量代码的你——把这两段程序的原理、可跑通的代码和最容易翻车的坑，一次讲清楚。
2. 四步相移法程序：四张条纹图算出包裹相位的Python实现与参数细节
2.1 四步相移公式：为什么差分能消掉背景，又为什么要用atan2
四步相移的基本思路是：向被测表面投射正弦条纹，用相机采集四张有固定相移的调制条纹图。每张图像的像素强度可以统一写成：
I_i(x, y) = A(x, y) + B(x, y)cos[φ(x, y) + δ_i]，i = 1, 2, 3, 4
其中A是背景光强（环境光加直流分量），B是条纹调制幅度，φ是物面相位，δ是人为附加的相移量。四步法选取 δ 为 0、π/2、π、3π/2 四个步进。把I2和I4相减：
I4 − I2 = Bcos(φ + 3π/2) − Bcos(φ + π/2) = 2Bsin(φ)
再看I1和I3的差：
I1 − I3 = Bcos(φ) − Bcos(φ + π) = 2Bcos(φ)
两个式子一比，背景项A和调制项B全部消掉，只剩下正切关系。理论上 φ = atan[(I4−I2)/(I1−I3)]，但工程代码里几乎没人这么写。原因有二：分母接近0时商会被噪声放大到失真；更重要的是 atan 的值域只有(−π/2, π/2)，会把第二、第三象限的相位判错。用 atan2(I4−I2, I1−I3) 则能根据分子分母的符号判定象限，输出完整的(−π, π]区间，也就是常说的“包裹相位”。
这里就带出包裹相位的核心特征：反正切结果永远被截断在(−π, π]，真实连续相位一旦超出这个范围，就会出现从+π到−π的跳变。物体表面越陡、条纹频率越高，跳变越密集，这也是后面必须做解包裹的根本原因。四步法不是步数越多越好，三步法少采集一张但对噪声更敏感，五步Hariharan法能校正线性相移误差但多花一帧时间；工程上四步是灵敏度和鲁棒性的折中，也是我默认的起点。
2.2 最小可跑通的Python示例：从模拟条纹图到包裹相位
我习惯先用模拟数据验证算法正确性，再上真实相机，这样能把算法问题与设备问题分开。下面这段代码生成一个带噪声的相位面，完整走一遍四步相移提取包裹相位：
import numpy as np # 模拟一个 256x256 的相位面：倾斜 + 抛物线弯曲 h, w = 256, 256 x, y = np.meshgrid(np.arange(w), np.arange(h)) true_phase = 0.003 * x + 0.004 * (y - 128) ** 2 # 生成四步相移条纹图，加高斯噪声模拟相机暗区 frames = [] for delta in [0, np.pi / 2, np.pi, 3 * np.pi / 2]: img = 128 + 90 * np.cos(true_phase + delta) noisy = img + 3 * np.random.randn(h, w) frames.append(noisy.astype(np.float64)) I1, I2, I3, I4 = frames # 四步相移公式：提取包裹相位 wrapped_phase = np.arctan2(I4 - I2, I1 - I3) # 调制幅度：像素点的条纹对比度，可用来生成掩码 modulation = 0.5 * np.sqrt((I1 - I3) ** 2 + (I4 - I2) ** 2)说明两点。第一，模拟相位面的系数是刻意选的：0.004*(y−128)² 在图像上下边缘已经让相位跨越了60多弧度，必然出现大量2π跳变，正好演示后续解包裹的必要性。第二，噪声必须加在强度域再做反正切，才能真实反映噪声对相位提取的影响；如果直接给 true_phase 加噪声再代入公式，结果会偏乐观，掩盖掉实际系统的毛刺。
跑完这段，wrapped_phase 的取值范围是(−π, π]，用 imshow 查看会看到密集的彩色条纹，这就是包裹相位图。modulation 是每个像素的条纹对比度，可以拿它生成掩码：设定一个阈值，把暗区、阴影、低对比度区域全部排除，后面给解包裹用。常见做法是把 mask = modulation > 阈值 写成函数，因为每个系统阈值不同，用Otsu自动分割也能凑合。
2.3 参数取舍：相移误差、条纹频率和噪声的三角关系
四步相移在理想假设下很干净，但实际系统里三个参数最容易让结果翻车。
第一是相移量精度。机械位移台或数字投影的相移步长如果偏离90°超过2°到3°，四步法会引入周期性误差，误差频率是条纹频率的两倍，表现为相位图上的波浪纹理。这个误差在解包裹阶段不会消失，反而被全局最小二乘“摊平”成更大范围的波纹。如果系统误差是线性且可重复的，换Hariharan五步相移能在很大程度上抵消；如果是随机抖动，只能靠缩短曝光与位移间隔、加固机械安装来压。
第二是条纹频率。条纹越密，对表面细节的分辨率越高，但相邻像素的相位差也越大。四步相移要求一个条纹周期内至少4个像素采样，实际操作建议保持6到8个像素每周期。超过这个密度，局部梯度会逼近π/像素，后端的解包裹算法再强也无能为力，因为信息已经被欠采样抹掉了。
第三是信噪比。反正切运算本质上是除法，会把强度噪声传递到相位域，尤其在I1−I3接近0的地方，也就是相位接近±π/2的区域，噪声被放大得最明显。我一般会在相移前对四张条纹图做σ=1到2像素的高斯平滑，代价是损失一点边缘锐度，但得到的相位图毛刺会少很多。曝光时间尽量让条纹图中间灰度落在100到180之间，别让投影仪gamma或相机饱和来添乱。
3. 最小二乘相位解包裹程序：包裹相位到连续相位的DCT求解
3.1 解包裹问题：2π跳变与相位展开的本质
包裹相位看起来“有规律地断开”，但每一处跳变都不确定是真实变化还是2π折叠。比如一个斜坡表面，真实相位从0平缓涨到20π，包裹结果却是每涨2π就跳回原点，形成锯齿。如果只有单个像素跳变，手动加2π就能恢复；但真实测量图里跳变无处不在，逐点处理根本不现实。
早期的主流做法是路径跟踪法：从一个起点出发，沿某条扫描路径展开，遇到相邻像素差值大于π就补偿2π。这个方法速度快，却有一个致命缺点：噪声或局部坏点会造成错误补偿，而且错误会沿着路径一路传播到最后，形成一条条像拉链一样的直线误差。更麻烦的是，物体表面有突起、孔洞、遮挡时，路径被切断，展开结果直接分家，不同区域的相位各自漂移。
最小二乘解包裹的思路完全不同：不找路径，而是把整个相位场当成一个全局优化问题，找一个连续相位场，使它的梯度在最小二乘意义下最接近包裹相位梯度。噪声被当作整体误差分摊到全场，不会因为某个坏点判断失误导致整条路径崩掉。这就是它成为工业界标配的原因：牺牲了一点速度，换来对噪声和遮挡的强得多。
3.2 最小二乘解包裹的数学模型：离散泊松方程
把包裹相位记为φ_w。相邻像素的相位差经过包裹算子W处理，得到“缠绕梯度”：
Δx(i, j) = W(φ_w(i, j+1) − φ_w(i, j)) Δy(i, j) = W(φ_w(i+1, j) − φ_w(i, j))
其中 W(t) = atan2(sin t, cos t)，结果落在(−π, π]。这个包裹算子很关键：原始相位差可能接近3.9 rad，但经过W处理后得到一个在(−π, π]内的等价值，它把由于2π折叠造成的“假梯度”重新映射成真实微小梯度。
然后建立优化目标：
min_φ Σ [φ(i, j+1) − φ(i, j) − Δx(i, j)]² + [φ(i+1, j) − φ(i, j) − Δy(i, j)]²
对每个像素的φ求偏导并令其为零，整理后得到离散泊松方程：
φ(i+1, j) + φ(i−1, j) + φ(i, j+1) + φ(i, j−1) − 4φ(i, j) = ρ(i, j)
其中ρ就是Δx和Δy的散度：
ρ(i, j) = Δx(i, j) − Δx(i, j−1) + Δy(i, j) − Δy(i−1, j)
方程左边是标准的五点拉普拉斯算子，右边是自己能算出来的已知量。直接解这个线性系统要面对 H×W 个未知数，矩阵稀疏，但用普通求逆依然慢。更快的方法是用离散余弦变换在频域求解：在 Neumann 边界条件下，DCT 基函数恰好是 Laplace 算子的特征函数，把泊松方程变换到频域后，每个频率分量只做一次除法即可。
3.3 DCT快速求解：完整可运行代码与参数说明
下面是我项目里一直在用的版本，加了详细注释方便裁剪。依赖 numpy 和 scipy，Python 3.8 以上都能跑：
import numpy as np from scipy.fft import dct, idct def least_squares_unwrap(wrapped_phase: np.ndarray, mask: np.ndarray | None = None) -> np.ndarray: """ 最小二乘相位解包裹（DCT 快速解法） 参数 ---------- wrapped_phase : 2D array，包裹相位，取值范围 (-π, π] mask : 2D bool array，有效区域为 True；None 表示全图有效 返回 ------- unwrapped : 2D array，解包裹后的连续相位 """ h, w = wrapped_phase.shape if mask is None: mask = np.ones((h, w), dtype=bool) # 1) 计算相邻像素的包裹相位差 dx = np.zeros((h, w), dtype=np.float64) dy = np.zeros((h, w), dtype=np.float64) dx[:, :-1] = np.angle(np.exp(1j * (wrapped_phase[:, 1:] - wrapped_phase[:, :-1]))) dy[:-1, :] = np.angle(np.exp(1j * (wrapped_phase[1:, :] - wrapped_phase[:-1, :]))) # 2) 无效区域的梯度置零 dx *= mask dy *= mask # 3) 计算散度 rho rho = np.zeros((h, w), dtype=np.float64) rho[:, :-1] += dx[:, :-1] rho[:, 1:] -= dx[:, :-1] rho[:-1, :] += dy[:-1, :] rho[1:, :] -= dy[:-1, :] # 4) 二维 DCT 正变换 rho_dct = dct(dct(rho, axis=0, norm='ortho'), axis=1, norm='ortho') # 5) 频域求解泊松方程 rows = np.arange(h, dtype=np.float64).reshape(-1, 1) cols = np.arange(w, dtype=np.float64).reshape(1, -1) denom = 2.0 * (np.cos(np.pi * rows / h) + np.cos(np.pi * cols / w) - 2.0) denom[0, 0] = 1.0 # 避免除零 phi_dct = rho_dct / denom phi_dct[0, 0] = 0.0 # 直流分量置零 # 6) 逆 DCT 回到空间域 unwrapped = idct(idct(phi_dct, axis=0, norm='ortho'), axis=1, norm='ortho') return unwrapped参数说明里有三个地方值得单独讲。
denom 矩阵是二维泊松方程在 DCT 域的特征值，对应公式 2cos(πi/H) + 2cos(πj/W) − 4。行列索引从0开始，(0,0) 是直流分量，特征值为0，方程在直流分量上本来就无解，所以单独设成1做除法，最后再把结果强制置0。
norm='ortho' 是归一化选项，保证正变换和逆变换互为逆运算。如果漏写，结果会被放大一个与尺寸相关的倍数；写错成 'norm=None' 会出现整幅图偏亮或偏暗的假象。
边界条件的处理藏在梯度计算里：dx[:, :-1] = np.angle(...) 把最后一列的梯度留成0，表示图像右边界外没有梯度，对应 Neumann 边界条件，即边界外法向梯度为零。上下左右四个边界都遵循这个约定，这是DCT解法与傅里叶解法最根本的区别，也是它能正确处理矩形状图的原因。
调用方式很简单，把第2章的 wrapped_phase 直接扔进去：
# 生成一个全有效区域的 mask；更稳妥是用 modulation 阈值生成 mask = np.ones((256, 256), dtype=bool) unwrapped = least_squares_unwrap(wrapped_phase, mask) # 减去一个参考点，让相位从 0 开始，方便后续标定 unwrapped = unwrapped - unwrapped[100, 100]mask 想偷懒就直接传 None，效果等同于全图有效。但真实场景里我强烈建议传 mask，否则阴影区的随机相位会被当成有效信号，拖累整个解包裹结果。性能方面，1024×1024 的相位图在这段代码上大约耗时 0.2 秒，比迭代加权最小二乘快一个数量级，足够用在线测量流水线。
4. 相位解包裹避坑指南：噪声、掩码与边界常见问题的排查
4.1 解包裹结果出现周期波浪纹理：相移误差的标志性症状
现象：解包裹后的相位场整体趋势是对的，但表面覆盖着一层周期性的波浪纹，像水面涟漪，周期和条纹周期呈倍数关系。单看包裹相位图时这层波纹已经存在，解包裹不会消除它，反而因为全局最小二乘把它“平均”到更大的区域，看起来更像噪声云。
原因：头号嫌疑是四步相移的相移量不精确。假设实际步长是90°+ε，四步法引入的相位误差近似为(ε²/2)·sin(2φ)，误差频率正好是条纹频率的两倍。第二位嫌疑是条纹响应非线性，比如投影仪gamma、相机传感器饱和，也会引入谐波，症状几乎一样。
解决：先用五步Hariharan相移做对照实验。如果波浪纹明显减弱，就说明是相移误差，可以改用五步法或对位移台做标定。如果波浪纹没变化，就要查gamma：对投影仪做gamma查找表校准，或对采集图像做幂次校正。直接对相位图做高通滤波能掩盖症状，但这是治标不治本，我一般不推荐。
4.2 mask边界处大台阶：散度计算与掩码传播的坑
现象：mask把物体圈起来后，解包裹结果在mask边界内外出现明显台阶，物体边缘和背景不是平滑过渡，而是硬生生跳了一段。有些方向的台阶特别明显，换个mask阈值台阶位置还会变。
原因：mask传入后，无效区域的梯度乘了0，但散度计算里仍隐含“边界内外梯度一致”的假设。DCT解法采用的是Neumann边界条件，它默认图像边界外法向梯度为0；当mask内部有空洞或形状不规则时，等价于在算法里额外塞进很多“内部边界”，这些边界的法向梯度被错误置0，导致边界两侧相位不连续。
解决：最实用的做法是对mask做形态学腐蚀，把物体边界往里收缩2到3个像素，让算法避开最不可靠的边缘区域。想要更精细就改加权最小二乘：有效区域权重设1，无效区域设非常小的值比如1e−6，然后迭代求解。注意mask里不要留孤立小洞，先用闭运算填掉，否则会在解包裹结果里形成小圆斑。
4.3 解不唯一导致常数偏移：直流分量与参考平面对齐
现象：同一组数据，把mask改大或改小后再跑，解包裹相位整体数值差了一个常数；在不同机器上跑同一个npy文件，结果也不一样。
原因：离散泊松方程解出来的相位场，任意加一个常数仍然满足方程。DCT解法在(0,0)分量强制把解置0，等价于让整个相位场的均值（严格说是最小二乘意义下的零均值）为0。所以只要mask或数据范围变了，这个“零均值”基准就跟着变。
解决：工程上不要盯着绝对相位值看，而是对比相位差。标准做法是放一个参考平面：分别采集参考平面和物体的条纹图，各跑一遍四步相移加解包裹，然后算 Δφ = φ_obj − φ_ref，两边的常数偏置在减法中自然抵消。如果只有一次测量，也可以手动指定一个已知平坦的点做零位，比如把物体台面的平均相位减掉。
4.4 梯度超过π/像素：最小二乘救不回来的混叠问题
现象：物体某处特别陡，比如90°台阶的侧壁，解包裹相位在那一段出现大块错误区域，而且错误会污染周围几十个像素，形成放射状伪影。
原因：无论哪种解包裹算法，第一步都要从包裹相位求相邻差。如果真实相邻相位差超过π，包裹算子W会把它错误映射到(−π, π]内的另一个值，比如3.9 rad会被认成−2.38 rad。这个错误直接进入散度ρ，泊松方程给出的解自然对不上。这不是算法缺陷，是采样不足。
解决：只有两条路。一是降低条纹频率，让最陡处每像素相位变化小于π，最好留到π/2以下，需要重新设计投影条纹。二是改用多频外差或时间相位展开：先用低频条纹做粗测，再用高频条纹细化，粗测结果给细测提供“这是第几个周期”的先验知识，从根上消除歧义。我建议在项目初期就预留多频接口，别等出现陡坡了再返工。
5. 从包裹相位到高度图：残差验证与进阶方向
5.1 残差热图验证解包裹质量
解包裹成功与否不能只看结果是否“平滑”。最客观的验证是把解包裹结果重新包裹，再与输入对比：
re_wrapped = np.angle(np.exp(1j * unwrapped)) residual = np.angle(np.exp(1j * (re_wrapped - wrapped_phase))) # 画残差热图，颜色越接近 0 越好残差热图里，如果只有原来2π跳变残留的窄条，说明解包裹是自洽的。如果出现大面积弥散的非零斑块，说明梯度计算或mask处理有问题，回到第4章排查。
5.2 参考平面标定：从相位差换算高度
条纹投影里，解包裹相位本身不是高度。平行光路情况下，高度与相对相位近似成正比：z(x, y) = k·(φ_obj − φ_plan) = k·Δφ。k由系统几何参数决定，常见做法是放一个已知高度的标准块，比如5mm，量出它的Δφ，k = 已知高度 / Δφ。经验更足时，用标定板的几个高度做最小二乘拟合，把k和光路的常数偏置一起估计出来。
5.3 进阶路线：从最小二乘到多频外差
如果你面对的表面大部分平缓、偶尔有台阶，上面这套最小二乘解包裹完全够用。如果台阶多、遮挡多、或者要绝对相位，那就接多频外差。做法是对三个频率的条纹各跑一遍四步相移，得到三个包裹相位图；频率两两相减生成新的“低频”包裹相位，再继续相减，直到获得一个全场无歧义的相位，然后逐级回推。代码层面，四步相移函数和最小二乘函数都能复用，只是调用流程多套一层循环。
我个人的习惯是先把模拟数据上的残差调到几乎全零，再上真实相机。这样设备调试时出现的任何异常都能直接怀疑光学或机械环节，而不是让算法背锅。希望帮到你。
本文还有配套的精品资源，点击获取
