简介面向光学测量、结构光三维成像与相位重建领域开发者的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, j1) − φ_w(i, j)) Δy(i, j) W(φ_w(i1, j) − φ_w(i, j))其中 W(t) atan2(sin t, cos t)结果落在(−π, π]。这个包裹算子很关键原始相位差可能接近3.9 rad但经过W处理后得到一个在(−π, π]内的等价值它把由于2π折叠造成的“假梯度”重新映射成真实微小梯度。然后建立优化目标min_φ Σ [φ(i, j1) − φ(i, j) − Δx(i, j)]² [φ(i1, j) − φ(i, j) − Δy(i, j)]²对每个像素的φ求偏导并令其为零整理后得到离散泊松方程φ(i1, j) φ(i−1, j) φ(i, j1) φ(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 和 scipyPython 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有效区域为 TrueNone 表示全图有效 返回 ------- unwrapped : 2D array解包裹后的连续相位 h, w wrapped_phase.shape if mask is None: mask np.ones((h, w), dtypebool) # 1) 计算相邻像素的包裹相位差 dx np.zeros((h, w), dtypenp.float64) dy np.zeros((h, w), dtypenp.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), dtypenp.float64) rho[:, :-1] dx[:, :-1] rho[:, 1:] - dx[:, :-1] rho[:-1, :] dy[:-1, :] rho[1:, :] - dy[:-1, :] # 4) 二维 DCT 正变换 rho_dct dct(dct(rho, axis0, normortho), axis1, normortho) # 5) 频域求解泊松方程 rows np.arange(h, dtypenp.float64).reshape(-1, 1) cols np.arange(w, dtypenp.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, axis0, normortho), axis1, normortho) return unwrapped参数说明里有三个地方值得单独讲。denom 矩阵是二维泊松方程在 DCT 域的特征值对应公式 2cos(πi/H) 2cos(πj/W) − 4。行列索引从0开始(0,0) 是直流分量特征值为0方程在直流分量上本来就无解所以单独设成1做除法最后再把结果强制置0。normortho 是归一化选项保证正变换和逆变换互为逆运算。如果漏写结果会被放大一个与尺寸相关的倍数写错成 normNone 会出现整幅图偏亮或偏暗的假象。边界条件的处理藏在梯度计算里dx[:, :-1] np.angle(...) 把最后一列的梯度留成0表示图像右边界外没有梯度对应 Neumann 边界条件即边界外法向梯度为零。上下左右四个边界都遵循这个约定这是DCT解法与傅里叶解法最根本的区别也是它能正确处理矩形状图的原因。调用方式很简单把第2章的 wrapped_phase 直接扔进去# 生成一个全有效区域的 mask更稳妥是用 modulation 阈值生成 mask np.ones((256, 256), dtypebool) 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 进阶路线从最小二乘到多频外差如果你面对的表面大部分平缓、偶尔有台阶上面这套最小二乘解包裹完全够用。如果台阶多、遮挡多、或者要绝对相位那就接多频外差。做法是对三个频率的条纹各跑一遍四步相移得到三个包裹相位图频率两两相减生成新的“低频”包裹相位再继续相减直到获得一个全场无歧义的相位然后逐级回推。代码层面四步相移函数和最小二乘函数都能复用只是调用流程多套一层循环。我个人的习惯是先把模拟数据上的残差调到几乎全零再上真实相机。这样设备调试时出现的任何异常都能直接怀疑光学或机械环节而不是让算法背锅。希望帮到你。本文还有配套的精品资源点击获取