用MATLAB实现傅里叶变换轮廓提取与三维重建：从频域滤波到相位解包裹

1. 从“傅里叶变换轮廓”这个说法聊起：它到底在解决什么问题

我最早接触到“傅里叶变换轮廓”这个词，是在一次图像处理大作业的题目里。当时第一反应是：傅里叶变换和轮廓提取有什么关系？轮廓不应该是梯度、边缘检测算子（比如Sobel、Canny）干的事吗？等到自己动手在MATLAB里跑了一遍，才理解这个词背后的两层含义——第一层是用傅里叶变换提取图像中物体的二维边缘轮廓，第二层是光学测量领域里一个专门的三维轮廓重建技术，全称叫傅里叶变换轮廓术（Fourier Transform Profilometry，简称FTP）。

这篇文章不打算只讲理论公式，重点是把两条路线都在MATLAB里走通：先用频域方法把一张普通照片的边缘轮廓“抠”出来，再沿着同一套思路，用条纹投影加傅里叶变换的方法恢复物体表面的三维高度轮廓。无论你是正在做图像处理大作业、想给毕设加一个亮点，还是纯粹想搞懂傅里叶变换在图像里到底“看见了什么”，这篇文章都适用。我会尽量把每一步的MATLAB代码、参数含义和调试心得都写清楚，不需要你有多深的理论基础，能跟着敲完代码就算赢。

先说一个容易劝退新手的点：傅里叶变换的数学形式看起来吓人，但用MATLAB处理图像时，你只需要记住三个关键操作——fft2做二维变换、fftshift把零频挪到中心、ifft2把频域结果变回空间域。剩下的所有细节，都是围绕这三个操作在做文章。下面我就从最基础的频谱可视化开始，一步步把“轮廓”这件事讲透。

2. 先把图像的频谱看懂：为什么轮廓藏在“高频”里

2.1 用MATLAB把频谱“拍”出来

在动手提取轮廓之前，我强烈建议你先做一件事：把图像的频谱打印出来看一眼。这一步没有任何技术含量，但是能帮你建立最关键的直觉——傅里叶变换到底把图像变成了什么。

假设你有一张灰度图lena.png（用自带图片或其他图片都可以），下面这段代码可以把频谱显示成一个直观的灰度图：

img = imread('lena.png'); if size(img, 3) == 3 img = rgb2gray(img); end img = im2double(img); % 转成double，方便后续运算 F = fft2(img); % 二维傅里叶变换 F_shifted = fftshift(F); % 把零频分量移到频谱中心 magnitude = log(abs(F_shifted) + 1); % 加1再取log，压缩动态范围 figure; subplot(1, 2, 1); imshow(img, []); title('原始图像'); subplot(1, 2, 2); imshow(magnitude, []); title('频谱幅度（对数显示）');

你可能已经见过类似的频谱图：中间一个十字亮带，四周零星分布着一些亮点和暗纹。中间最亮的点就是直流分量，它代表图像的“平均亮度”——整张图最基础的灰度基调，和轮廓没有任何关系。从中心往外延伸的亮线，对应图像里方向性很强的纹理结构；而四周那些散落的亮点，才是我们真正关心的边缘和细节信息。

这里有一个关键点：log(abs(F) + 1)里的+1不是随便加的。原始傅里叶变换的数值跨度极大，直流分量可能上百万，高频分量可能只有零点几，直接显示会变成一片黑或一片白。取对数能把动态范围压缩到人眼可分辨的范围。你在任何教程里看到的那种漂亮的频谱图，几乎都是做了对数变换的。

2.2 边缘在频域里到底是什么形态

理解了频谱长什么样，接下来就要回答核心问题：图像中的轮廓（也就是边缘）在频域里对应什么？

我自己的理解是这样的：一幅图像可以看成无数个不同频率、不同方向的正弦光栅叠加而成。平坦区域亮度变化缓慢，对应低频；边缘区域灰度在几个像素内发生剧烈跳变，对应高频。所以一个物体的轮廓，本质上是大量高频分量沿着边缘法线方向叠加的结果。

这里有个非常实用的推论：如果我想提取轮廓，只要把低频成分抑制掉，只保留高频成分，再逆变换回空间域，那么平坦区域会变得灰灰的，边缘位置则特别亮。这就是频域高通滤波提取边缘的基本原理，也是“傅里叶变换轮廓”最朴素的一种实现。

为了验证这个直觉，你可以做一个简单的实验：构造一张纯粹的棋盘格图像，然后看它的频谱。棋盘格的边缘是水平和垂直方向的，所以频谱里会在水平和垂直方向出现一堆离散的亮点。如果旋转棋盘格，频谱里的亮点也会跟着旋转。这个实验我建议亲手做一次，做过之后你对“频域方向对应空间方向”这句话的体会会深很多。

3. 动手实现：用频域高通滤波把轮廓“抠”出来

3.1 滤波器的选择：理想高通、Butterworth还是高斯

既然思路明确了，接下来就是设计滤波器。在MATLAB里做频域滤波，本质上就是构造一个和图像同样大小的矩阵，这个矩阵在低频位置取值接近0，在高频位置取值接近1，然后把它和频谱逐点相乘。设计这个“取值矩阵”的方式有很多种，我在实际项目里主要用过以下三类：

理想高通滤波器的MATLAB实现最简单，但我不建议你在真实图片上用它，因为频域里的一次“硬切”会在空间域产生振铃效应——轮廓周围出现一圈一圈的明暗纹路，跟鬼影一样。更稳的做法是Butterworth滤波器，它的过渡带是渐变的，阶数n控制过渡的陡峭程度。n越大，过渡越陡，但也越容易出现轻微振铃，所以一般取2到4比较合适。

下面是一段完整的、可以直接运行的MATLAB脚本，用Butterworth高通滤波器提取图像的边缘轮廓：

img = imread('lena.png'); if size(img, 3) == 3 img = rgb2gray(img); end img = im2double(img); [M, N] = size(img); F = fft2(img); F_shifted = fftshift(F); % 构造频率坐标网格，范围是[-1, 1] u = linspace(-1, 1, N); v = linspace(-1, 1, M); [U, V] = meshgrid(u, v); D = sqrt(U.^2 + V.^2); % 每个点到频谱中心的距离 D0 = 0.1; % 截止频率，可调参数 n = 2; % Butterworth阶数 H = 1 ./ (1 + (D0 ./ (D + eps)).^(2*n)); % 高通滤波器 % 频域滤波 G_shifted = F_shifted .* H; % 逆变换回空间域 g = real(ifft2(ifftshift(G_shifted))); % 归一化到[0,1]方便显示 g = (g - min(g(:))) / (max(g(:)) - min(g(:))); figure; subplot(1, 3, 1); imshow(img, []); title('原始图像'); subplot(1, 3, 2); imshow(magnitude, []); title('频谱'); subplot(1, 3, 3); imshow(g, []); title('高通滤波提取的轮廓');

执行完这段代码，你会看到一个“勾勒”了人物主要轮廓、但背景相对较暗的结果。这正是我们想要的效果：平坦的区域被抑制，边缘被保留并放大。

3.2 从灰度轮廓到干净的二值轮廓

imshow(g, [])里显示的是一幅灰度图像，边缘位置的亮度高，但还不够“干净”。如果要做进一步分析（比如计算轮廓长度、识别物体形状），通常还需要把灰度结果转成二值图，也就是把边缘像素变成白色、其他像素变成黑色。

这一步看起来简单，实际上有一个大坑。如果你直接对g用全局阈值（比如im2bw(g, 0.5)），结果往往会很糟糕——因为不同区域的边缘亮度差异很大，部分弱边缘会被误杀，而一些纹理丰富的地方又会变成一大片白斑。我一般会先用形态学操作做后处理，流程如下：

% 自适应阈值二值化 bw = imbinarize(g, 'adaptive'); % 去除小面积噪声区域 bw = bwareaopen(bw, 20); % 形态学细化，让轮廓变成单像素宽度 skel = bwskel(bw); figure; imshowpair(bw, skel, 'montage'); title('左：二值化结果 右：细化后的骨架轮廓');

关于bwskel这个函数，我想多说两句。它是MATLAB R2019b之后推出的骨骼化函数，用来把二值区域细化为单像素宽的骨架。对轮廓提取来说，这个操作特别有价值——因为频域高通滤波恢复出来的边缘往往有“重影”现象（边缘位置有两条接近的线，甚至更宽），直接测宽度或长度都会不准确。bwskel能把这些粗线条压成一条线，方便后续做几何测量。

还有一点必须提醒：如果你的MATLAB版本比较老，imbinarize可能不可用，可以用graythresh加im2bw的替代方案。具体来说就是level = graythresh(g); bw = im2bw(g, level);。效果差不太多，但imbinarize的自适应模式在光照不均匀的图片上明显更稳。

3.3 参数选择的实操心得

这部分是纯经验，代码里看不出，但往往决定了你的轮廓质量。

第一个参数是截止频率D0。它的物理意义是：把频谱中心半径D0范围内的低频成分全部压下去。D0太小，滤波保留的低频过多，背景会残留大量灰度渐变，轮廓不突出；D0太大，又会把边缘本身的中频信息也削掉，轮廓变细甚至断裂。我第一次做的时候，D0取了0.3，结果人物面部轮廓几乎消失了，只留下发丝级别的细节。后来慢慢试，发现0.05到0.15之间通常比较合适，你可以从D0 = 0.1开始尝试，再根据结果每次增减0.02。

第二个容易忽略的点是在构造函数时给分母加一个eps。因为频谱中心点的距离D恰好为0，如果不加这一项，0.1 / 0会直接得到Inf的提示甚至NaN，整个滤波矩阵就废了。加eps是数值计算里的常规操作，目的是避免除零，代价是中心点极小范围的精度损失，对结果没有任何可感知的影响。

第三个经验是：如果原图里有规律性纹理（比如衣服上的格纹、墙上的砖缝），高频滤波后这些纹理的边缘也会被当成轮廓提取出来，形成大量“假边缘”。这时候不要急着调滤波器，先用imgaussfilt(img, 2)对原图做一个轻微高斯平滑，把细碎纹理抹掉再进频域处理，效果往往立竿见影。这个操作相当于在频域滤波之前，先在空间域做了一次预处理。

4. 进阶玩法：傅里叶变换轮廓术（FTP）的三维重建

前面讲的是用傅里叶变换提取二维图像中的平面轮廓，但“傅里叶变换轮廓”这个词在光学测量领域还有另一个完全不同的含义——通过分析一个被物体调制过的条纹图案，恢复物体的三维表面轮廓。这就是傅里叶变换轮廓术。

4.1 FTP的基本原理：条纹是“载波”，物体是“包络”

这个技术最有意思的地方在于，它把光学测量问题转化成了信号解调问题，思路和收音机调频广播有异曲同工之妙。

你想象一下：拿一台投影仪，往被测物体表面投一组竖直方向的正弦条纹，条纹的亮度按I(x, y) = a(x, y) + b(x, y) * cos(2π f0 x + φ(x, y))分布。当物体表面有高低起伏时，条纹看起来就弯曲了——凸起的地方条纹向左弯，凹陷的地方向右弯。这里的φ(x, y)就是物体高度引入的相位调制，它恰恰就是我们要恢复的东西。

关键是这个φ被“编码”在载波频率f0上，就像广播信号被编码在特定频道上一样。如果我们对这张条纹图像做傅里叶变换，频谱上会出现三个明显的东西：中心是零频和背景，左右两侧各有一个“旁瓣”，这两个旁瓣就是把φ带着走的“频道”。为了提取φ，做法就是对图像做傅里叶变换，然后用滤波器把其中一个旁瓣单独抠出来，搬回原点（相当于解调），再做一次逆傅里叶变换，得到的复数信号的相位角就是tan^(-1)(Im/Re)，即包裹相位。

下面这个流程是FTP的经典三件套，每一行代码都有明确含义：

% stripe是采集到的条纹图像，已归一化到[0,1] [M, N] = size(stripe); F = fft2(stripe); F_shifted = fftshift(F); % 构造带通滤波器，只保留右侧旁瓣 % 旁瓣中心大致在 (M/2, N/2 + f0*N 对应位置)，f0是条纹频率 mask = zeros(M, N); center_x = round(N/2 + f0 * N); % 右旁瓣中心 radius = round(f0 * N * 0.5); % 滤波半径 % 这里用高斯带通更稳 [X, Y] = meshgrid(1:N, 1:M); mask = exp(-((X - center_x).^2 + (Y - M/2).^2) / (2*radius^2)); % 提取旁瓣并搬到原点 filtered = F_shifted .* mask; filtered = circshift(filtered, [0, - (center_x - N/2)]); % 逆变换得到复数信号 complex_signal = ifft2(ifftshift(filtered)); phase_wrapped = angle(complex_signal); % 包裹相位，范围 [-pi, pi]

拿到phase_wrapped之后还没有结束，因为angle函数返回的相位被限制在[-π, π]之间，而真实的物体高度导致的相位变化往往远超这个范围。这就好比你开车记里程表，走完一圈表盘又归零——包裹相位就是这种“归零”的相位，需要“解包裹”才能得到连续的真实相位。

4.2 相位解包裹和高度映射

MATLAB自带一个一维解包裹函数unwrap(p)，它按“相邻点相位差超过π就加减2π”的规则把跳变修正掉。二维场景下可以逐行做unwrap，再逐列做unwrap，虽然算法朴素，但多数仿真实验够用：

phase_unwrapped = unwrap(phase_wrapped, [], 2); % 按行解包裹 phase_unwrapped = unwrap(phase_unwrapped, [], 1); % 再按列解包裹

解包裹之后得到的是连续相位场Δφ(x, y)，它和物体高度h(x, y)是线性关系：h = (λ / (2πd)) * Δφ，其中λ是条纹波长，d是投影仪和相机之间的几何参数。更严谨的做法是标定，但对大作业或课程设计来说，用一个已知高度的标准平面先测一次，反推出比例系数k = h_known / Δφ_measured，之后待测物体高度直接乘k就行。这个方法简单粗暴，但足够你交一份像模像样的实验报告。

我在第一次跑FTP仿真时踩过一个特别典型的坑：旁瓣滤波半径选得太小，导致解调出来的相位在物体边缘剧烈抖动。原因是旁瓣在频域里并不是一个干净的圆点，它周围有因为物体遮挡产生的连续频谱泄露。你如果发现重建出来的高度图边缘有“波浪纹”，八九不离十是带通滤波器半径太小。解决办法是把radius放大到0.7 * f0 * N左右，或者改用二维汉宁窗做平滑截断。

4.3 为什么说FTP比相移法更适合入门

做三维轮廓测量，除了FTP还有相移法。相移法精度更高，但要至少采集三张不同相移的条纹图，而且要求环境稳定；FTP最大的优势是单张图就能恢复出相位信息，非常适合同一个物体在动态变化场合的测量。MATLAB做FTP仿真还有一个额外好处：所有中间结果——频谱、滤波后的旁瓣、包裹相位、解包裹相位——都能直接用imshow一步一图地显示出来，每一步都看得见摸得着，这对理解整个算法链路帮助极大。

5. 调试中的高频“坑”：我在MATLAB里实际踩过的6个问题

这部分我打算集中写问题排查。原因很简单：算法的原理看一遍就懂，但代码第一次跑不出来、或者跑出来的结果很奇怪，才是真正劝退人的地方。以下问题我全部用MATLAB实际验证过，按出现频率从高到低排列。

5.1 结果图像全黑或全白

这个问题的几乎都是数据类型惹的祸。傅里叶变换处理的是double类型数据，值域可能远超过[0,1]，如果你直接用imshow(g)显示，MATLAB会把数据范围当作[0,1]，于是一堆负值或大于1的像素直接被截断成黑点或白点。解决方案永远是先归一化：g = (g - min(g(:))) / (max(g(:)) - min(g(:)))。我见过太多人漏掉这一步，对着全黑图怀疑人生。

5.2 边缘有严重的明暗“涟漪”

振铃效应。多半是你用了理想高通滤波器，或者Butterworth阶数太高（超过5阶）。在频域里突然把某些频率分量砍到零，相当于在空间域卷积了一个带旁瓣的核，于是边缘两侧出现周期性明暗条纹。解决办法优先换高斯高通，或者把Butterworth阶数降到1到2之间，并接受过渡带变宽的事实。

5.3 提取出的轮廓不是一条线，而是“两条线”

尤其是白色物体在黑色背景上，轮廓边缘会同时出现内侧和外侧两条响应，看起来像物体被描了双线。原因在于频域高通滤波同时保留了边缘两侧的正负跳变，它们叠加后互相独立。处理方法就是前面提到的bwskel骨骼化：先二值化，再细化。如果骨骼化后还是有分叉，加一步skeleton_endpoints相关处理，或者用bwmorph(skel, 'spur', 3)去掉小毛刺。

5.4 频域滤波后图像坐标“翻转”了

这类问题多半和fftshift/ifftshift的配对有关。规则只有一条：变换之前用了fftshift，变换之后必须用ifftshift，不能搞混。ifftshift不是fftshift的逆吗？对二维矩阵来说它们大部分时候等价，但当矩阵维度是奇数时两者结果不同。所以安全起见，写代码时就固定用ifftshift，别每次凭感觉写。

5.5 大尺寸图片处理时内存爆炸

这是每个用MATLAB处理图像的人迟早会遇到的事。一张4K图像大概是800万像素，fft2和后续的矩阵运算会产生若干个double类型的中间变量，一个变量60多MB，同时存在五六个就逼近500MB内存。处理办法有两个：一是把图片缩到合适大小再处理，img = imresize(img, [512, 512])；二是分段处理，对条纹图做FTP时在Y方向逐行处理，也能显著降低内存占用。

5.6 频域滤波后整体图像变暗

这不是bug，而是物理规律。大幅度抑制低频后，图像的平均亮度几乎必然下降。高通滤波的初衷就是只保留变化剧烈的部分，平坦区域像素趋近于零是正常结果。如果你需要“轮廓叠加在原图上”的效果，可以先把高通结果二值化，再把白色边缘点映射回原图变成红色或绿色。我常做的是用overlay = imoverlay(img, bw, [1, 0, 0]);这样的思路，视觉效果很直观。

6. 我对傅里叶变换轮廓提取的几点实践心得

如果只让你记住一句话，那就是：傅里叶变换不是用来“找轮廓”的标准工具，而是用来“理解轮廓”的强有力工具。标准边缘检测算子（Sobel、Canny）在空间域直接算梯度，速度快、效果好；傅里叶方法的好处在于思路清晰、可解释性强、能从全局频域视角理解图像结构，而且它是很多高级图像处理技术（比如FTP三维测量）的基础。所以我的建议是，做课程设计或大作业时，用频域方法多画几张图、多写一点原理分析，拿分效果远比直接调一个edge()函数好得多。

另外我还想强调一个习惯：每次修改参数之后，一定要把原始图像、频谱图、滤波结果图并排放一起对比。我见过不少同学调参数完全靠猜，看到结果变丑就随便改一个数字，最后陷入“调来调去也不知道自己在干嘛”的循环。正确做法是每次只改一个参数，认真对比三张图的变化，记录下参数和效果之间的因果关系。这个过程看起来慢，但积累下来，你对频域滤波的直觉就会非常准。

最后分享一个小技巧：写MATLAB脚本时，把所有关键参数集中放在文件头部，用注释标注每个参数的含义和常用区间。这样不管是自己后面回来调试，还是给老师、同学演示，都能快速定位问题。代码本身不值钱，值钱的是你调试过程中的判断力——而这恰恰是从一次次“运行失败”和“结果不对”里长出来的。