首页 / 资讯中心 / 文章详情

互补格雷码与相移码结合的结构光相位解包裹实现

互补格雷码与相移码结合的结构光相位解包裹实现 ★ FEATURED ARTICLE
距离我第一次在结构光系统里把互补格雷码和相移码也就是多步相移正弦条纹结合起来解包裹相位已经过去好几年了。到现在凡是需要稳定恢复绝对相位的场景我还是首选这套组合原理不绕、代码不难、鲁棒性有惊喜。这篇把包裹相位出现的原因、互补格雷码的差分解码逻辑、还有一整套 Python 实现完整写出来适合正在做结构光三维重建、投影编码方案验证或者论文复现时被相位 unwrap 卡住的人直接参考。先说清楚这套方法到底解决了什么问题。相移法算出来的相位天生是“包”着的范围只有 (-π, π]每过一个条纹周期就跳变一次因此无法直接映射到三维坐标。要拿到连续、单调的绝对相位就得知道每一个像素落在第几个条纹周期里这一步就是相位展开。互补格雷码干的正是这件事它负责把“第几个周期”这个整数序号稳定地算出来再叠加到相移得到的包裹相位上最终得到绝对相位。整条链路里互补格雷码的差分解码是核心亮点本文会从原理到代码逐层拆开讲。我不会只贴一堆代码。这里面参数的取舍、为什么会踩坑、怎么排查比代码本身更值钱。建议你准备一个能跑 numpy 和 matplotlib 的环境跟着后面的完整脚本过一遍再回到真实投影仪和相机上去验证。1. 原理拆解包裹相位为何需要互补格雷码来辅助解缠1.1 相移法算出来的相位为什么是“包”着的四步相移是最经典的相位提取方式。投影仪依次投射四幅正弦条纹相移量分别是 0、π/2、π、3π/2相机同步采集到四幅光强图I1 R·(A B·cos φ) I2 R·(A B·cos(φ π/2)) I3 R·(A B·cos(φ π)) I4 R·(A B·cos(φ 3π/2))其中 R 是物体表面反射率A 是环境背景光B 是条纹调制幅度φ 是我们想求的相位。把两组差分信号做比值然后取反正切φ atan2(I4 - I2, I1 - I3)这样算出来的相位值只落在 (-π, π] 区间里。可以把 φ 理解为“真实相位对 2π 取余数后的结果”所以它在图像上呈现出一圈一圈的锯齿状每一条锯齿对应一个条纹周期。数学上这叫包裹相位英文 wrapped phase而包裹相位和真实相位之间的关系是Φ φ 2π·k这里 k 就是周期序号对应这个像素处于第 k 条条纹周期内。只要 k 定错了相位就会整体偏移 2π反映在三维重建结果上就是一条条明显的“台阶”断裂。相位展开的核心任务说白了就是把这个整数 k 求准。1.2 互补格雷码到底补在哪儿差分解码原理求 k 最直接的办法是给每个周期编一个号再用额外的条纹图把这个号“拍”出来。格雷码是一种相邻码字只差一位的二进制码。假设有 K 位格雷码就能编码 2^K 个周期投影 K 幅黑白条纹图之后每个像素可以得到一个 K 位码字把它解码成十进制就得到序号 k。传统格雷码的痛点在于二值化。因为解码时要在每幅格雷码图上判断当前像素是黑还是白最常见的做法是找个全局阈值比如大于 128 判 1、小于 128 判 0。这个做法在物体表面反射率均匀、环境光稳定的情况下问题不大但真实场景里物体表面有深色、浅色、金属高光、局部阴影同一个码值在不同位置的绝对亮度可能差好几倍全局阈值必然误判。一旦某一位判反k 可能错出十万八千里解包出来的相位直接跳变。互补格雷码针对这个痛点做了一个非常聪明的改动每个码位不只投一幅纯格雷码图还额外投一幅亮度分布完全相反的反码图。也就是说同一像素在正码图里如果亮那么它在反码图里一定暗。解码时不再找阈值而是逐像素比较正反两幅图的亮度正码比反码亮判 1正码比反码暗判 0。用公式表达就是位值 sign( I正 - I反 )这个“差”的好处一眼就能看出来。理想情况下正码图亮度为 P A R·C反码图亮度为 Q A R·(1-C)C 的取值是 0 或 1。两式相减得到 P - Q R·(2C - 1)。环境背景光 A 被直接消掉了反射率 R 虽然有乘性影响但 R 恒为正数不会改变差值的正负符号。所以不管物体表面是黑是白不管环境光强还是弱只要这个位置还在工作范围内符号判断的结论都一致。这就是互补格雷码比传统格雷码更稳的底层原因。代价也很明确投影幅数增加了。K 位格雷码原本投 K 幅图互补方案要投 K 幅正码加 K 幅反码总共 2K 幅。对于一个六位格雷码系统就是多投六幅图采集时间变长。所以它更适合静态或者准静态场景对于高速运动物体可能需要换成多频外差这种更省投影幅数的方案。2. 关键参数设计条纹周期、相移步数、格雷码位数怎么搭配2.1 周期数、分辨率与格雷码位数的计算公式参数设计的起点是投影仪横向分辨率。假设投影仪宽度为 W 像素我们把正弦条纹沿横向铺开让每个条纹周期占 PERIOD_PX 个像素那么全幅一共就有T W / PERIOD_PX个周期。格雷码要覆盖所有周期码位数 K 需要满足K ceil( log2(T) )如果 T 刚好等于 2^K那是最理想的解码出来的码值 0 到 T-1 全部有效。如果 T 不是 2 的幂实际编码时只使用前缀的 T 个码字解码后要把大于等于 T 的码值掩膜掉。比如投影仪宽度 1280周期宽度 32 像素T 40K 6能编码 64 个码字但实际只有 0 到 39 有效40 到 63 不能出现出现就说明解码异常或者该像素被遮挡。周期宽度怎么选这里有个工程经验。PERIOD_PX 太小比如小于 4正弦条纹在投影仪和相机分辨率下采样点太少相位求解误差明显增大PERIOD_PX 太大周期太粗相位展开后一个周期内能分辨的细节变少三维点云的空间分辨率被限制。我一般取 10 到 40 像素之间具体看物体尺寸和系统的工作距离。本文的演示代码取 PROJ_W 1280、PERIOD_PX 20这样 T 64K 6正好是 2 的幂省去掩膜问题。2.2 相移步数选择三步、四步还是五步相移步数的选择影响相位求解精度和投影幅数。三步相移最少只需要三幅图抗运动伪影最好但公式里依赖调制幅度 B如果投影仪或物体表面导致 B 不稳定误差会直接进入相位。四步相移多投一幅但使用 (I4-I2)/(I1-I3) 这种差分组合可以抵消恒定背景光和部分线性系统误差工程上最常用。五步相移进一步使用最小二乘拟合对投影仪 gamma 非线性更不敏感但投影幅数多采集时间变长。我的建议很直接第一套系统用四步公式简单好调试也容易推广到三步或五步。代码里的相移量可以写成循环想改成 N 步只动一个参数不需要重构逻辑。表里面列出了常规对比方便做方案评审时直接抄。相移步数投影幅数优点缺点适用场景3 步3 幅投影最少速度快对调制幅度和噪声较敏感动态测量、素材受限4 步4 幅差分抵消背景鲁棒性好比三步多一幅静态/准静态场景通用首选5 步5 幅抗 gamma 非线性最好投影幅数最多高精度测量、离线标定2.3 光栅方向与坐标约定光栅方向决定了相位展开后映射到哪个投影仪坐标。如果正弦条纹沿横向变化也就是亮度只随 x 变化那么相位展开得到的绝对相位就对应投影仪的水平坐标反之条纹沿纵向变化则对应垂直坐标。实际三维重建系统通常需要在水平、垂直两个方向分别投影和解包才能建立物体表面点到投影仪像素的完整坐标映射。本文先以水平条纹为例代码里我把生成函数的 vertical 参数预留好了改成纵向条纹时只需要把相位计算公式里的 x 换成 y其余逻辑完全一致。写代码之前一定要把坐标方向理清楚不然最后做相位到坐标映射时会发现解出来的坐标轴反了。3. 完整实现从条纹生成到绝对相位求解3.1 图像生成正弦条纹、格雷码条纹与互补条纹先用 Python 把需要投影的条纹图生成出来。正弦条纹的核心是余弦函数相位随 x 线性变化再加上相移量 offset。格雷码条纹则要先构造 K 位格雷码序列然后把某个位置对应的码位值横向铺开成一个宽度为 PERIOD_PX 的色块。互补条纹很简单就是正码图取反色。import numpy as np import matplotlib.pyplot as plt from scipy.ndimage import map_coordinates, median_filter PROJ_W 1280 PROJ_H 720 PERIOD_PX 20 # 单个条纹周期占 20 像素 T PROJ_W // PERIOD_PX # 周期数 64 K int(np.log2(T)) # 格雷码位数 6 PH_STEPS 4 # 四步相移 def generate_gray_codes(k): 生成 2^k 个 k 位格雷码返回十进制列表 return [i ^ (i 1) for i in range(2 ** k)] def sinusoidal_fringe(phase_shift, verticalFalse): 生成正弦条纹图phase_shift 为相移量 if not vertical: x np.arange(PROJ_W) phase 2 * np.pi * x / PERIOD_PX phase_shift img 127.5 127.5 * np.cos(phase) img np.tile(img, (PROJ_H, 1)) else: y np.arange(PROJ_H) phase 2 * np.pi * y / PERIOD_PX phase_shift img 127.5 127.5 * np.cos(phase) img np.tile(img[:, None], (1, PROJ_W)) return np.clip(img, 0, 255).astype(np.uint8) def gray_fringe(k, bit_idx, invertedFalse): 生成第 bit_idx 位的格雷码条纹图invertedTrue 时生成反码 codes generate_gray_codes(k) img np.zeros((PROJ_H, PROJ_W), dtypenp.uint8) for x in range(PROJ_W): idx min(x // PERIOD_PX, T - 1) bit (codes[idx] (k - 1 - bit_idx)) 1 if inverted: bit 1 - bit img[:, x] 255 if bit else 0 return img这里有个细节值得注意正弦条纹的数值范围是 0 到 255但我在生成时故意把直流分量放在 127.5调制幅度也是 127.5而不是直接乘一个 0 到 255 的随机值。这样做的目的是给后面模拟采集留足动态范围避免负值被截断。格雷码条纹直接使用 0 和 255 两个极值以便和正弦条纹分开处理。投影之前通常还会做 gamma 校正这一步在真实系统里做仿真流程里先忽略。3.2 模拟“投影-采集”过程加入物体高度相位没有真实硬件时可以用一个二维高斯凸起来模拟物体高度场。物体高度会让投影条纹发生横向偏移偏移量的大小和相位偏移成正比。只要对原始条纹图做一次水平位移变换就能模拟出“物体被投影后相机看到”的效果。这个模拟虽然简单但足够验证解包裹算法的正确性。def simulate_object_phase(): 生成模拟物体相位场形状为高斯凸起单位rad y, x np.mgrid[0:PROJ_H, 0:PROJ_W] h 80.0 * np.exp(-((x - 640) ** 2 (y - 360) ** 2) / (2 * 200 ** 2)) return 2 * np.pi * h / PERIOD_PX def warp_by_phase(img, obj_phase): 按相位场对图像做水平位移模拟物体对条纹的调制 h, w img.shape shift obj_phase / (2 * np.pi) * PERIOD_PX x np.arange(w) coords np.tile(x, (h, 1)) shift coords np.clip(coords, 0, w - 1) rows np.tile(np.arange(h)[:, None], (1, w)) return map_coordinates(img, [rows, coords], order1, modenearest) obj_phase simulate_object_phase() sin_ims [] for i in range(PH_STEPS): fringe sinusoidal_fringe(i * 2 * np.pi / PH_STEPS) captured warp_by_phase(fringe, obj_phase) captured np.clip(captured, 0, 255).astype(np.uint8) sin_ims.append(captured) gray_ims [] inv_gray_ims [] for b in range(K): g gray_fringe(K, b, invertedFalse) gi gray_fringe(K, b, invertedTrue) gray_ims.append(warp_by_phase(g, obj_phase).astype(np.uint8)) inv_gray_ims.append(warp_by_phase(gi, obj_phase).astype(np.uint8))warp 函数是这段代码里最容易出错的地方。map_coordinates 的输入坐标系必须匹配图像数组的行列顺序很多新手在这里把 x 和 y 换反结果图像变得一团糟。我的习惯是先在注释里写清楚“第一维是行第二维是列”再用 rows 和 coords 两个数组显式构造坐标索引肉眼检查一遍再跑。另外 modenearest 是为了处理物体边缘位移到图像外的问题代价是边缘一圈会出现拉伸假象真实采集时边缘本来就应该用掩膜剔除。3.3 包裹相位与互补格雷码解码四步相移求包裹相位直接套公式。注意 atan2 的返回值范围是 (-π, π]如果后续需要映射到 [0, 2π)可以加一条条件判断但展开绝对相位时不影响。互补格雷码解码的关键步骤是差分二值化正码亮度大于反码判 1否则判 0。把 K 个位值逐位移位组合成一个整数再通过 gray_to_binary 转换成周期索引。def gray_to_binary_vector(g): 把格雷码整数向量转换成二进制十进制索引 inv np.zeros_like(g) tmp g.copy() while True: inv ^ tmp tmp tmp 1 if np.all(tmp 0): break return inv I1, I2, I3, I4 [arr.astype(np.float64) for arr in sin_ims] wrapped_phase np.arctan2(I4 - I2, I1 - I3) decoded_bits np.zeros((PROJ_H, PROJ_W), dtypenp.int32) for b in range(K): diff gray_ims[b].astype(np.float64) - inv_gray_ims[b].astype(np.float64) bit (diff 0).astype(np.int32) decoded_bits (decoded_bits 1) | bit period_index gray_to_binary_vector(decoded_bits) period_index np.clip(period_index, 0, T - 1)为什么直接看“正码比反码亮”而不是单独对正码设阈值前面原理部分已经解释过这里我再强调一次实操意义真实采集时相机拍到的每一幅图都会受到环境光照影响背景光在差分里相互抵消反射率 R 作为正系数不会翻转符号。所以即使你拿到一整套灰度级差异很大的图像只要正反码对齐良好diff 的符号分布依然稳定。这是整个算法鲁棒性的命根子。有一件事必须提醒diff 的大小本身没有直接意义只有符号有用。如果物体表面有强烈镜面高光正反图都接近 255差值会变得很小噪声容易把符号翻转。遇到这种情况可以加一条“低置信度标记”当 abs(I正 - I反) 小于某个阈值时认为该像素为不可靠点后续用邻域插值或区域生长做修补。3.4 相位展开、物体相位恢复与误差验证有了周期索引和包裹相位绝对相位就是两者相加Φ φ 2π·k因为模拟中物体高度引起的相位偏移就是 obj_phase理论上可以拿“解包出来的相位”减去“原始条纹的基线相位”再和真值 obj_phase 比较从而验证整条链路有没有解错。base_phase 2 * np.pi * np.arange(PROJ_W)[None, :] / PERIOD_PX unwrapped_phase wrapped_phase 2 * np.pi * period_index reconstructed_obj unwrapped_phase - base_phase valid (obj_phase 1.0) # 只看物体区域 err reconstructed_obj[valid] - obj_phase[valid] rmse np.sqrt(np.mean(err ** 2)) print(f物体相位恢复 RMSE: {rmse:.4f} rad)我实际跑这个仿真时物体中心区域的 RMSE 通常在 0.02 到 0.05 rad 之间边缘因为 map_coordinates 插值的关系会偏大一些。如果你跑出来的 RMSE 很大先不要怀疑算法优先检查几个点周期索引是不是整体错位了 1 个周期格雷码解码出来的码字是不是和真实条纹周期对不上以及 warp 模拟物体相位时正负方向有没有取反。这里建议把 period_index 单独隔离开来看一眼伪彩图正常情况下它应该是一个按周期递增的阶梯状平滑图像任何不该出现的“雪花点”都能通过它快速定位。可视化阶段可以画三张图包裹相位图、解码周期索引图、绝对相位图。包裹相位有大量锯齿跳变周期索引是阶梯状整数图绝对相位是连续渐变图。三张图放在同一行对比如果算法正确它们之间的递进关系会非常直观。4. 工程落地中的坑与排查技巧4.1 我踩过的几个坑第一坑是码字边界和相位周期边界没有对齐。格雷码条纹的每个色块宽度是周期宽度但相移法求出的包裹相位在周期边界处会从 -π 跳到 π两个系统的边界如果偏差了一两个像素解码出的周期索引和真实相位之间就会产生随机的 2π 错位。这个问题在仿真里不显眼因为图像是理想生成的但在真实系统中非常常见。解决思路有两个层面一是在投影仪和相机标定时做好几何对齐二是用相位信息反过来补偿码字边界也就是在像素坐标的亚像素层面对周期索引做插值修正。第二坑是投影仪 gamma 非线性。很多廉价投影仪的光强响应不是线性的输入 128 可能实际输出亮度偏离理想值很多。正弦条纹被非线性扭曲后会出现高次谐波四步相移解出来的包裹相位带上周期性波纹误差。互补格雷码的差分解码对 gamma 不敏感因为它只看符号但正弦相位的波纹误差会直接变成三维重建的波浪面。应对办法是做一个全局 gamma 标定投影一组均匀灰度图测量相机响应建立查找表生成条纹前把灰度值预畸变回去。第三坑是黑色高光区域。物体表面如果是黑色的反射率 R 接近零正反码差值 P - Q 会小到和传感器噪声量级相当符号不稳定。这种情况下解码出来的码字会出现随机跳变。我目前的策略是同时保留一个“可靠度图”把 diff 绝对值低于阈值的像素直接标记为无效再后续处理里用空间邻域插值补齐。不要试图用更复杂的滤波器硬解标记无效反而更干净。第四坑藏在仿真和真实的差距里。仿真时物体相位是平滑的高斯凸起解包裹很难出错真实系统里物体表面可能有遮挡、台阶、高光格雷码条纹在遮挡边缘会产生半影区相机拍到的不是干净的 0 或 255而是过渡灰度。这时候差分解码的符号判断可能仍然有效但位值已经不可靠。投影仪离焦会让黑白边界模糊解决办法是适当缩小投影仪光圈或者选择景深更大的投影系统。4.2 常见问题速查表整理一张排查表遇到异常直接按表格对照能省不少调试时间。现象可能原因排查和解决思路解包后的绝对相位整体跳变 2π周期索引错位一个周期检查格雷码解码顺序正反码方向是否反了相位图像出现随机“雪花点”黑色低反射区或高光导致差分符号翻转增加可靠度阈值标记无效点后邻域插补物体边缘出现条纹撕裂边缘相位超出采集范围扩大视场或对边缘做掩膜剔除正弦相位上有周期性波纹投影仪 gamma 非线性做 gamma 预畸变或改用五步相移格雷码码字边界与相位边界不对齐投影仪相机几何误差标定投影仪坐标做亚像素边界修正重建出的物体高度反向条纹方向或相位符号取反检查 warp 方向和 atan2 参数顺序4.3 调优思路与扩展方向这套方案已经能支撑大多数静态场景的三维重建但它不是终点。我后期做过两个方向的扩展效果都不错。第一个方向是在周期索引的边界处结合空间连续性做修正。互补格雷码解码后偶尔还会残留少量合理范围内的误码可以利用相邻像素相位连续这一先验检测绝对相位在 x 方向上相邻像素的跳变量级。如果相邻像素绝对相位差超过了 π 的若干倍说明中间很可能有码字错误再用邻域中值替换。这种“时间编码 空间校验”的组合比单纯依赖某一边有效得多。第二个方向是把周期索引从整数提升到亚像素。格雷码给出的周期序号是整数相移相位在周期内部是连续的两者拼在一起后周期边界处会出现细微的不连续。改进办法是先用相移结果拟合出一个亚像素精度的码字边界再对周期索引进行线性插值。这一步能让重建点云在周期边界上的“棱线”效应明显减弱尤其在做高精度面型测量时值得做。如果场景对投影幅数极度敏感比如要拍运动物体可以换成多频外差方案用一个低频率的相位偏移去解另一个高频相位的包裹。但多频外差的噪声敏感度通常比格雷码高我更倾向于静态精度优先的场合继续用互补格雷码。最后分享一点实操心得这套代码我从仿真跑到真实硬件中间踩了无数回坑。说实话真正让互补格雷码方案发挥价值的不是那多投的几幅反码图而是它逼着你把“亮度”和“码值”彻底解耦。做结构光系统时我们最怕的不是噪声而是把噪声错误地当成信号。差分解码天然提供了一层“置信度”差值大码值可靠差值小就要警惕。如果你正准备从零搭建结构光三维测量系统我强烈建议先把本文的仿真链路完整跑通加上遮挡、反射率不均、噪声这些干扰观察一下周期索引在什么条件下会崩。这一步做到位之后再碰硬件你面对现实成像环境时会有完全不一样的排查手感和信心。
阅读完成 · 觉得有帮助?
咨询建站