前阵子有个研究生来问我，MATLAB做结构光三维重建到底该从哪儿入手。很多新手一上来就翻论文，三频四步相移法、多频外差、包裹相位展开这些术语看得头大，真正能跑的代码却拼不出一套。其实这套方法远没有想象中那么神秘：投影仪往被测物体上依次投三组不同频率的四步相移正弦条纹，相机同步拍照，再从12幅图像里把每个像素的绝对相位解出来，最后利用标定好的相位-高度关系映射到三维坐标。整个过程用MATLAB实现，代码量一百多行，关键是把每一步的数学原理对应到矩阵运算上。这篇文章不绕弯子，直接按"原理→选型→代码→避坑→进阶"的顺序带你把整条链路跑通，适合正在做结构光课题、或者想用MATLAB快速验证算法的同学参考。
1. 先搞懂三频四步相移法到底在解决什么问题
1.1 结构光的本质是"投影-拍摄-反推"
结构光三维重建说白了就是一句话：给被测物体表面打上已知图案，再用相机观察图案被物体形状扭曲成了什么样，反推出表面深度。投影仪投正弦条纹，相机拍回来的条纹在物体表面发生了形变，这个形变量里就藏着高度信息。
从数学模型上看，物体表面某一点的高度h，最终会转化成该点的相位偏移Δφ。换句话说，如果我们能算出每个像素的相位值，再通过预先标定的"相位-高度"映射关系，就能直接得到高度。这就是主动条纹投影类方法的基本逻辑，而三频四步相移法就是目前实验室和工业界最常用的相位提取与展开方案之一。
1.2 四步相移如何逐像素计算相位
相移法解决的是"怎么从光强图里算出相位"的问题。投影仪投一个正弦条纹，相机的接收光强可以写成：
I_n(x,y) = A(x,y) + B(x,y)·cos(φ(x,y) + δ_n)
其中A是背景光强，B是条纹调制幅度，φ就是我们要解的相位，δ_n是每次投影的相移量。所谓四步相移，就是依次投影δ = 0、π/2、π、3π/2四张图，得到四张光强图。
把这四张图按公式一推，奇迹就发生了：
- I_4 - I_2 = 2B·sin(φ)
- I_1 - I_3 = 2B·cos(φ)
所以φ = atan2(I_4 - I_2, I_1 - I_3)，背景A和调制B全部被消掉了。这正是相移法的精妙之处，它不需要知道物体表面是亮是暗，只要条纹清晰，就能逐像素算出相位。
1.3 三频外差为什么能消除2π歧义
四步相移算出来的是包裹相位，范围只在[-π, π)之间。真实相位如果超过了一个周期，就会被折叠回来，造成2π歧义，这就是所谓"相位展开"问题。
空间相位展开算法（比如MATLAB里的unwrap）沿着像素路径做展开，遇到光滑表面没问题，一旦遇到物体表面有台阶、断裂、遮挡，展开路径就会出错，误差像多米诺骨牌一样传下去。
三频外差是时间相位展开算法中的一种。它投影三种不同频率的条纹，分别求出三种包裹相位，再通过差频运算构造出一个频率为1的绝对相位。频率为1意味着整幅图像只有一个条纹周期，相位从0到2π单调递增，不再存在歧义。然后用这个绝对相位作为基准，逐级把高频相位展开。因为每个像素独立运算，完全不依赖邻域信息，天然免疫表面不连续的问题。
1.4 和其他方案放在一起看
做结构光的朋友应该都听过格雷码+相移，或者光栅投影。
格雷码+相移的思路是用格雷码给每个条纹周期编号，相移负责周期内部细分，两者结合也能得到绝对相位，但要对物体投影一二十甚至更多幅图。
三频四步只需要12幅图，效率明显更高，而且三个频率投影顺序对精度的影响有明确的物理解释。
从解算稳定性来看，多频外差属于时间编码，每个像素的展开结果只依赖自身，不像空间展开那样会"跨像素传染"错误。这也是为什么三频四步在动态测量和高精度工业检测中如此受欢迎。
2. 环境准备、系统选型与三频参数设计
2.1 MATLAB环境与硬件准备
本文所有代码基于MATLAB R2025a编写，实际上R2020之后的版本都能直接运行，不需要额外工具箱。图像处理工具箱里有imshow、imagesc这些可视化函数，即使没有，换成plot也能跑，核心计算全是矩阵运算。
如果你是新装MATLAB，建议检查一下许可证路径是否配置正确，有些人装了附加功能却无法调用，多半是许可证没有关联到MathWorks账户，这是个常见问题。
硬件方面，实验室搭建结构光系统通常需要一台工业相机、一台DLP投影仪和同步触发装置。DLP投影仪的优势是灰度线性度和刷新率高，工业相机负责采集条纹。如果是做算法验证，不急于上真实硬件，先用下面的模拟数据把算法链路跑通，再逐步替换成真实图像，调试成本和心理压力都会小很多。
2.2 三频条纹参数到底怎么定
三频外差的核心约束是：三个频率的差频之间还得差出1。假设频率为f1、f2、f3（周期数），那么：
- f12 = f1 - f2
- f23 = f2 - f3
- f123 = f12 - f23 = 1
这个f123必须等于1，才能让整幅图像只有一个周期，从而消除相位歧义。
| 组合方案 | f1 | f2 | f3 | f12 | f23 | f123 | 特点 |
|---|---|---|---|---|---|---|---|
| 方案A | 70 | 64 | 59 | 6 | 5 | 1 | 差频小，展开最稳定，入门首选 |
| 方案B | 64 | 56 | 49 | 8 | 7 | 1 | 差频中等，抗噪和精度比较均衡 |
| 方案C | 100 | 85 | 71 | 15 | 14 | 1 | 最高频高，细节精度好，但噪声敏感 |
从表格能看出，f3越大，条纹越密，相位对高度越敏感，测量精度越高；但噪声也会被放大，解包裹的错误概率同步上升。我平时在实验室习惯先用方案A把整个流程调通，再根据实际精度需求增大f3。
2.3 相移步数为什么选四步而不是三步或五步
三步相移只需要3幅图，速度最快，但对相移误差和投影仪Gamma畸变比较敏感。
五步相移对谐波的抑制能力更强，因为它的误差项在频域中被推到了更高频，代价是多投两幅图，计算量也略大。
四步处在中间：投影数量适中，公式形态最简洁，对常见的线性相移误差有天然的抑制效果。工业上如果扫描速度允许，很多人宁可多投一步用五步求稳；但如果是在做动态测量、需要控制总投影幅数，四步是业界最常见的平衡选择。
3. 完整MATLAB实现：从条纹生成到绝对相位解算
3.1 第一步：模拟物体与生成三频四步条纹
为了让代码可以直接运行，我用一个高斯形状模拟三维物体，生成三种频率、每个频率四步相移的条纹图像。物体高度造成一个与频率成正比的相位偏移，这样保证三个频率的绝对相位严格成比例，满足三频外差的展开条件。
%% 参数设置 H = 480; % 图像高度 W = 640; % 图像宽度 freqs = [70, 64, 59]; % 三频条纹周期数 steps = 4; % 四步相移 obj_A = 12; % 物体最大高度 (mm) obj_sigma = 80; % 物体高斯宽度 (像素) D = 300; % 等效相位-高度常数 (mm) noise_sd = 0.01; % 高斯噪声标准差 rng(2025); % 固定随机种子，保证结果可复现 %% 生成模拟物体高度 [x, y] = meshgrid(1:W, 1:H); cx = W / 2; cy = H / 2; h_true = obj_A * exp(-((x - cx).^2 + (y - cy).^2) / (2 * obj_sigma^2)); %% 生成三频四步条纹，模拟投影仪投影 + 相机采集 I_stack = zeros(H, W, length(freqs) * steps); for fi = 1:length(freqs) f = freqs(fi); for n = 1:steps phase = 2 * pi * f * ((x - 1) / W + h_true / D) + (n - 1) * pi / 2; I = 0.5 + 0.5 * cos(phase); I = uint8(I * 255); % 量化为8bit灰度 I = double(I) / 255 + noise_sd * randn(H, W); % 加噪声 I_stack(:, :, (fi - 1) * steps + n) = I; end end我把物体高度用高斯函数来模拟，中心最高12毫米，四周平滑过渡到0，看上去像一座小山丘。这个模型的好处是平滑且数学形式固定，后面恢复出高度后可以直接和真值做像素级对比。条纹生成时用的是归一化坐标(x-1)/W，这样第一列的初始相位为0，后续全局常数校准会省很多事。
加噪声之前先把图像量化到uint8，是为了模拟真实相机8bit采集的动态范围。很多教材代码直接生成浮点条纹，这在模拟里当然能跑通，但一旦切换到真实图像，量化误差、噪声、暗电流这些因素全都会冒出来，提前在模拟里加入反而能帮你建立更准确的直觉。
3.2 第二步：四步相移提取三个包裹相位
有了三组四步条纹，解包裹相位的公式直接套用，每个频率独立计算，得到三个包裹相位图。注意atan2返回的范围是[-π, π)，这是后面所有展开运算的前提。
%% 四步相移提取包裹相位 phi_wrap = zeros(H, W, length(freqs)); for fi = 1:length(freqs) idx = (fi - 1) * steps; I0 = I_stack(:, :, idx + 1); I1 = I_stack(:, :, idx + 2); I2 = I_stack(:, :, idx + 3); I3 = I_stack(:, :, idx + 4); phi_wrap(:, :, fi) = atan2(I3 - I1, I0 - I2); end phi1 = phi_wrap(:, :, 1); phi2 = phi_wrap(:, :, 2); phi3 = phi_wrap(:, :, 3);这里容易踩的坑是：有人会写成atan2(I0-I2, I3-I1)，顺序反了导致相位符号反号。符号在解包裹时可能带来整体翻转，最终高度恢复就会变成凹陷变凸起。建议拿到包裹相位后，先画出来看一眼条纹走向，再做下一步，不要闷头往下算。
3.3 第三步：三频外差逐级展开绝对相位
三频外差的完整展开分两步：先通过两次差频得到频率为1的相位基准，再做全局常数校准；然后用这个基准逐级展开中频、高频。
%% 三频外差：差频得到等效频率 f12 = freqs(1) - freqs(2); f23 = freqs(2) - freqs(3); f123 = f12 - f23; % 必然等于1 phi12 = mod(phi1 - phi2, 2 * pi); % 等效频率 6 phi23 = mod(phi2 - phi3, 2 * pi); % 等效频率 5 phi123 = mod(phi12 - phi23, 2 * pi); % 等效频率 1 %% 全局常数校准：用unwrap展开基准相位并归零 phi123_unwrap = unwrap(phi123, [], 2); k_offset = round((0 - phi123_unwrap(cy, 1)) / (2 * pi)); Phi123 = phi123_unwrap + 2 * pi * k_offset; %% 用基准相位逐级展开 k23 = round((f23 / f123 * Phi123 - phi23) / (2 * pi)); Phi23 = phi23 + 2 * pi * k23; k3 = round((freqs(3) / f23 * Phi23 - phi3) / (2 * pi)); Phi3 = phi3 + 2 * pi * k3;很多初学者看不懂为什么要做"全局常数校准"，这里必须说透。unwrap得到的是相对相位，它可以是2π的任意整数倍偏移，而后续展开公式里clock周期的round操作依赖这个常数。如果常数偏了0.3、0.7这种值，round可能取错整数，导致整幅图的展开相位全部差2π，却还看不出哪里断了。校准的方法很直接：在参考平面位置，我们事先知道理论绝对相位应该是0，那就把这里的unwrap结果强行拉回到0，其他像素跟着平移。
理论上，对于完全无噪声的数据，不校准也能得到正确展开。但真实场景中噪声会让unwrap后的基准相位产生一个整体偏移，不校准就很容易翻车。这个细节是我在真实系统里踩过坑之后才理解的。
3.4 第四步：从绝对相位恢复高度并评估误差
最高频绝对相位Φ3里包含了飞机的载频（由x坐标引入）和物体高度（由h引入）。载频是我们自己设计的，可以直接减去，剩下的就是高度。
%% 从绝对相位恢复高度 h_est = D * (Phi3 / (2 * pi * freqs(3)) - (x - 1) / W); %% 误差评估 err = h_est - h_true; rmse = sqrt(mean(err(:).^2)); fprintf('高度重建RMSE = %.4f mm\n', rmse); %% 可视化 figure; subplot(1, 3, 1); imagesc(h_true); axis image; colorbar; title('真实高度'); subplot(1, 3, 2); imagesc(h_est); axis image; colorbar; title('重建高度'); subplot(1, 3, 3); imagesc(err); axis image; colorbar; title(sprintf('误差 (RMSE=%.4f)', rmse)); colormap jet;把真实高度、重建高度和误差放到一起看，如果RMSE在微米到毫米量级，说明相位解算链路是通的。我跑这个模型时，10毫米左右的高斯物体，误差通常在0.01毫米量级，主要来自于量化噪声和随机噪声。你可以试着把noise_sd调到0.05甚至0.1，观察RMSE如何上升，这会让你直观感受到噪声对三频外差算法的影响有多直接。
4. 从绝对相位到三维点云：标定模型与坐标映射
4.1 相位-高度映射的标定逻辑
前面的模拟代码用了一个解析表达式把相位换算成高度，但真实系统里相位和高度之间不是这么简单的线性关系。由于投影仪、相机之间存在夹角，镜头有畸变，实际关系通常是二次甚至更高次的非线性映射。
实验室最常见的做法是：把一块标准平面放在几个已知高度位置，比如用位移平台每次移动1毫米，在每个高度投影三频四步条纹，解出该高度的绝对相位分布。建立一个多组"高度-相位"对应点，然后用二次多项式拟合：
h = a0 + a1·Φ + a2·Φ²
有了系数a0、a1、a2，再拍摄任意物体时，直接代入相位Φ就能得到高度。如果系统畸变较大，可以改为逐像素拟合，每个像素单独拟合一组系数，计算量略大但精度提升明显。
MATLAB里的polyfit和polyval就能完成这项工作，几十行代码就能搞定标定流程。具体做法是在每个已知高度h_k下采集一组条纹，解出该平面的绝对相位图Φ_k(x,y)，然后对每个像素收集一系列(h_k, Φ_k)，用polyfit(x, y, 2)拟合二次曲线。
4.2 三维坐标生成与点云可视化
高度图已经算出来了，三维点云的X、Y坐标还需要相机内外参数映射。这里我用一个虚拟针孔相机模型来演示，相机的等效焦距fx、fy，主点坐标cx、cy，这些参数在真实系统里由标定得到。
%% 虚拟相机内参映射到三维点云 Z0 = 500; % 相机到参考平面距离 (mm) fx = 800; fy = 800; % 等效焦距 (像素) cx = W / 2; cy = H / 2; Z = Z0 - h_est; % 物体表面点的深度 X = (x - cx) .* Z / fx; % 相机坐标系下的X坐标 Y = (y - cy) .* Z / fy; % 相机坐标系下的Y坐标 figure; scatter3(X(:), Y(:), Z(:), 1, h_est(:), '.'); axis equal; colormap jet; colorbar; xlabel('X (mm)'); ylabel('Y (mm)'); zlabel('Z (mm)'); title('三维重建点云');点云的密度取决于图像分辨率，480×640的图像能生成约30万个点，足够看出物体形状了。需要导出PLY或者XYZ文件时，直接concatenate X、Y、Z矩阵并写文本文件即可。很多人会忽略轴的比例，导致点云看起来被压扁或拉长，记得设axis equal。
4.3 单目结构光方案的局限：为什么不能拍就完事
这套方案只有一个相机加一个投影仪，本质上是"单目结构光"。它能测出物体表面的深度Z，但前提是标定精度足够高，而且被测表面必须在相机和投影仪都能同时看到的区域。
如果物体表面有深凹槽，投过去的条纹会被自身的凸起挡住，相机拍不到，该区域就成了阴影盲区。这是单目结构光天生的硬伤，不是算法能解决的。要覆盖更完整的表面，通常需要多个相机视角，或者旋转物体多次测量。做实验时被阴影区折磨过的人，对这个限制一定深有同感。
5. 换到真实环境后最容易翻车的四个细节
5.1 投影仪Gamma畸变与二次谐波
投影仪输出亮度和输入灰度通常不是线性关系，很多投影仪的实际输出近似为输入的2.2次方。一个理想正弦条纹经过幂函数变换后，波形会被"压扁"，产生明显的二次谐波分量。
在四步相移中，二次谐波会让相位误差出现2倍频震荡，表现为重建表面上有规律的水波纹。这个误差在实验室里非常常见，很多人折腾半天以为是标定问题，其实是Gamma搞的鬼。
对策之一是对投影图像做预校正：预先给投影仪测一条"灰度-亮度"响应曲线，生成一个反函数查找表，投影前先反算一遍灰度。例如标定出Gamma系数为2.2，就把要投影的灰度做1/2.2次幂变换：
gamma = 2.2; lut = uint8(((0:255)./255).^(1./gamma).*255); I_prewarp = lut(I_proj + 1);在实际系统中，这个预校正步骤几乎能消除大部分低频相位波纹。如果你发现重建出来的平面始终有周期性的起伏，优先怀疑Gamma，而不是急着加滤波。
5.2 环境光、饱和与低调制度区域
四步相移理论上可以消掉背景A，但前提是背景A在每个像素上恒定，且没有像素饱和。如果环境光过强，相机的动态范围被压缩，亮区域的条纹对比度B很低，甚至直接饱和成纯白，这时I_3-I_1变成零，atan2的输入全是噪声。
解决办法有三个层面：一是物理层面，尽量在暗室或遮光环境下操作；二是曝光层面，先投影全黑和全白图，根据两图的灰度分布把曝光时间调到中间灰度不饱和；三是算法层面，利用调制度B做掩膜，把B过低的无效像素置为NaN。
调制度的估算方法很简单：B = sqrt((I1-I3)^2 + (I4-I2)^2)/2，如果它低于全局最大值的10%，几乎可以肯定该像素条纹质量不合格。
5.3 物体边缘遮挡造成的阴影与跳变
投影仪和相机之间有夹角，物体轮廓边缘在物理上必然有一侧是"相机能看到但条纹照不到"的阴影区。这些区域没有有效条纹信息，相移解出的相位是噪声，展开后会出现大量不连续的跳变点。
处理阴影没有万能药，但有三招可以组合使用：第一，用上述调制度掩膜把阴影区直接删掉；第二，只对物体表面一定梯度范围内的像素做后续三维重建；第三，如果边缘轮廓特别重要，考虑增加第二台投影仪或者转动物体多次测量，让每个表面区域至少在某个视角下光照完整。
我第一次用三频四步测一个齿轮时，齿轮齿根的阴影区跳变点密密麻麻，三维点云看起来像长了刺。加了调制度掩膜之后，所有刺全部消失，重建表面立刻干净了。这招的效果比任何滤波都立竿见影。
5.4 空间相位展开的适用边界
MATLAB自带的unwrap属于空间展开算法，沿着某一行路径展开包裹相位。它的前提是相邻像素的真实相位差小于π。表面连续时展开得很好，一旦相位跳变超过π，unwrap就会出错，而且错误会沿路径向后传播。
三频四步的优势就在于它是时间展开，每个像素独立运算，理论上不存在跨像素错误传播。但要注意，时间展开要求每个像素都能给出一系列有效的包裹相位，如果某个频率的条纹在某个像素上被噪声污染，对应的展开仍会出错。
因此，实际项目中通常用"时间展开为主，空间展开为辅"的组合：先由三频四步得到绝对相位，再用空间滤波去除少量孤立坏点，而不是单独依赖任何一种策略。
6. 还能往哪些方向延伸
6.1 用查表替代多项式拟合
二次多项式拟合在测量范围不大时够用，但系统畸变较大时，多项式的残差会很扎眼。更稳健的做法是先标定一大组"高度-相位"数据，然后直接建一张查找表，用线性插值查相位对应的精确高度。
MATLAB的interp1或者interp2都能实现这个功能。标定越密集，查表精度越高，成本只是多存几个数组。很多工业结构光系统最终都是查表法，因为它在保证精度的同时，代码逻辑极其简单，几乎没有模型假设偏差。
6.2 引入BP神经网络做非线性映射
如果相位和高度之间存在非常复杂的非线性关系，多项式拟合力不从心，可以考虑用BP神经网络直接学习这个映射。
输入是若干频率的绝对相位，输出是高度值。训练数据可以用标定平台在不同高度采集的相位图构建，用MATLAB的神经网络工具箱可以很快训出一个小型拟合网络。神经网络的优势是不需要预设模型阶数，天然适合处理相机畸变、Gamma残留等综合非线性误差叠加出来的复杂映射。
但注意，模型换相机、换投影仪甚至换摆放位置后，原来的神经网络很可能失效，需要重新采集数据训练。它适合系统长期固定、追求极致精度的情况。
6.3 用层次分析法优化投影与采集参数
三频四步在前面的代码里用了一组频率、一步曝光、固定投影亮度。如果你的应用对精度、速度、成本都敏感，可以在调参阶段用层次分析法来排优先级。
比如把条纹频率、曝光时间、投影亮度作为候选方案，把相位噪声、重建精度、处理速度作为评价指标，通过成对比较矩阵算出各方案的综合得分。这套方法不一定能给出全局最优解，但能帮你从一堆经验参数里理出主次关系。工程上，我倾向于先固定相机曝光，再用层次分析法在频率和亮度之间权衡，能让调试过程少很多纠结。
6.4 可替换的结构光方案盘点
如果做完三频四步，想探索其他结构光方案，主流方向有这些：
| 方案 | 投影数量 | 主要优势 | 主要限制 |
|---|---|---|---|
| 格雷码+相移 | 10-20幅 | 绝对编码，鲁棒性极高 | 投影数量多，速度受限 |
| 三频四步 | 12幅 | 精度高，速度适中 | 对Gamma畸变敏感 |
| 二值条纹离焦投影 | 8-12幅 | 适用于DLP高速投影 | 对焦精度要求高 |
| 散斑相关 | 1-2幅 | 动态测量，无需扫条纹 | 精度较低，依赖纹理 |
如果要做高速动态三维测量，二值条纹离焦投影加DLP是主流路线，但它的核心思想仍然是相移法，只是把正弦条纹换成了二值条纹，靠投影仪离焦模糊生成正弦效果。理解了四步相移的公式，再去理解二值条纹离焦方案会容易得多。
最后说一点实践层面的建议。我搭真实系统时，调试顺序基本固定：先看条纹图质量，再检查包裹相位是否光滑，然后确认展开相位的连续性，最后才看重建误差。如果直接盯着最终点云找原因，往往会被各种误差叠加搞到崩溃。先把模拟代码在MATLAB里跑通，把每个中间量都可视化一次，再去碰真实硬件，这套流程能帮你省下大量排错时间。三频四步的原理不难，代码也不复杂，真正有门槛的是你对待中间结果的耐心程度。
