做触觉传感器和三维形貌重建的人多半绕不开 GelSight Sensor 和它输出的 height-map。我最早接触这套传感器时天真地以为拿到表面法向图之后高度场无非是沿着像素路径积一下分拉一条跨平台式的 scanline 就能搞定。直到真的在真实数据上跑了一轮才发现噪声、边界和累积误差会把结果毁得没法看。后来我把问题改写成泊松方程配合 DST-II 算子做快速求解计算量从“跑一次要等到怀疑人生”降到“实时预览毫无压力”。这篇文章就讲讲这套方案的思路、实现和踩坑记录适合正在做触觉重建、机器人抓取感知或者想了解快速数值求解如何落到图像处理管线的朋友。1. 先搞清楚我们要解什么问题1.1 GelSight 怎么把“触觉”变成“图像”GelSight 的物理结构并不神秘一块软凝胶覆盖在透明的支撑板上凝胶表面喷涂了反光涂层。当物体压上去凝胶表面跟着物体形状发生微米级的形变像一面镜子。传感器内部用红、绿、蓝三组 LED 从不同方向照射这面“镜子”顶上再放一个相机拍下反射光。形变越大、局部倾斜越陡反射光进入相机的角度就变化越明显RGB 三个通道的亮度就会呈现规律性的变化。这背后的光学逻辑是局部表面法向决定了反光方向反光方向决定了相机能不能在某个 LED 的照射下看到高光。把三个通道的亮度分别理解成“来自三个方向的光照强度”就能通过比色法反推出每个像素处表面的倾斜分量。换句话说GelSight 直接测量的是一个梯度场而不是高度本身。它输出的原始信号经过标定后可以转成每个像素处的法向量或者是表面的斜率分量 gx、gy。这个阶段最重要的标记是我们要重建的 height-map 不是直接从图像里读出来的它像一个“积分结果”。传感器给的是导数的证据我们要从证据里还原原函数。于是问题天然落到数值计算上。1.2 height-map 的直接积分法为什么不行我第一次上手时用的是最直观的路径积分从图像左上角出发gx 是 x 方向的斜率gy 是 y 方向的斜率于是逐像素累加h[i][j1] h[i][j] gx[i][j] h[i1][j] h[i][j] gy[i][j]听起来没有任何问题可真实数据一跑就露馅了。第一GelSight 图像里的梯度场带有大量噪声每个像素的 gx、gy 都有小误差累加时误差被一路放大积分久了画面里全是条纹状的漂移。第二路径选择会直接影响结果从左上角积分得到的右下角高度和从右下角反着积分回来的数值完全对不上因为噪声让梯度场不再是一个“可积”的保守场。这就像你在一个起伏不平的地形图上沿着不同路线量海拔每条路线量出来的终点高度都不一样你会怀疑到底是地形有问题还是尺子有问题。更本质的原因是路径积分只利用了局部信息。每一个点的高度只依赖于它同一条路径上的前一个点它不知道侧向的梯度值和周围一大片区域传递过来的约束。在存在噪声时这种“一条路走到黑”的算法几乎是灾难。所以行业里很少直接用路径积分来做 GelSight 的形貌重建除非只是看一个粗略的剖面趋势。1.3 泊松方程把局部梯度场“拼”回全局高度场既然路径积分不靠谱那就换成全局优化思路。我们希望找到一个高度场 h使它的梯度 ∇h 在最小二乘意义下尽量接近传感器测出来的梯度场 g。也就是min ∫ || ∇h - g ||² dx dy对这个目标函数做变分等价于求解一个泊松方程∇²h div(g)这里 ∇² 是拉普拉斯算子div 是散度。右边 div(g) 完全由实测梯度场决定对 gx 求 x 方向差分加上 gy 求 y 方向差分。左边是标准二阶微分算子。整个方程的含义非常直观高度场的“弯曲程度”应该等于实测梯度场的“膨胀程度”而不再是逐点累加。这个形式的好处是它在全局范围内分配误差不会出现路径积分那样的沿路径漂移。即使某个局部梯度值带了噪声它只会影响该点附近的平滑度而不会顺着一条线一路污染到图像边缘。用图像处理的术语来说泊松方程像一个全局低通滤波器把不可积的梯度场“投射”到可积函数空间里投影之后剩下的残差就是噪声。所以剩下的核心问题只有一个如何快速、稳定地求解 ∇²h div(g)尤其是当图像分辨率来到百万像素级别的时候。2. 为什么选 DST-II 算子而不是硬解矩阵2.1 从拉普拉斯算子到线性方程组离散化之后拉普拉斯算子变成经典的卷积模板。二维情况用五点模板∇²h[i][j] ≈ (h[i1][j] h[i-1][j] h[i][j1] h[i][j-1] - 4h[i][j]) / Δ²把整个 height-map 摊成一个一维向量这个离散拉普拉斯就变成一个巨大的稀疏矩阵 A。泊松方程变成A h bb 就是散度项 div(g)。如果不考虑效率直接调一个稀疏线性求解器也能跑。但问题是当图像尺寸是 1024 x 1024 时未知量超过一百万即使 A 是稀疏矩阵通用求解器也要吃大量内存和时间。在实际项目里我经常要在一台工控机上做实时重建这种“通用但昂贵”的算法肯定没法用。换个角度想A 是五对角块矩阵它对应的离散泊松方程是一个结构极其规则的线性系统。对规则系统数学上早就有比通用求解器快几个数量级的方法——换到频率域去解。2.2 特征分解思维换到“频率域”去解方程回忆一下信号处理里解微分方程的常规操作对时间信号做傅里叶变换微分算子就变成乘以 jω于是微分方程变成一个代数方程。微分算子在频域里是对角的这个性质简直是为快速求解量身定做的。离散泊松方程也一样。如果我们能找到一组基使得离散拉普拉斯矩阵 A 在这组基下变成对角矩阵那么解方程就变成三步1. 把 b 变换到频域b_hat transform(b) 2. 在频域做逐点除法h_hat b_hat / lambda 3. 逆变换回空间域h inverse_transform(h_hat)整个流程里没有迭代没有矩阵求逆只有三次变换加一次逐元素除法。如果变换本身是 O(N log N) 的快速算法那整体求解速度比通用稀疏求解器快几个数量级尤其在高分辨率下优势越发明显。这其实就是“用特征分解理解快速求解”的核心思路。A 的维度再大只要它能被一组固定基同时对角化我们就可以绕开显式的矩阵存储和分解。离散余弦变换DCT和离散正弦变换DST之所以在图像编码里随处可见本质上也是因为它们能对角化某些特殊边界的二阶微分矩阵。2.3 DST-II 到底做了什么边界条件如何对号入座具体到 GelSight 重建我为什么不用更常见的 DCT-II而选了 DST-II这里有一个特别容易被忽略的细节边界条件。在标准的二维泊松重建里如果高度场四周都是自由边界也就是法向导数在边缘处为零最自然的选择是 DCT-II。DCT 的基函数是余弦函数它自带“边界外推平坦”的 Neumann 特征。但 GelSight 重建 height-map 时很多场景下我们并不是在整个二维平面上解方程而是沿着扫描方向一条剖面、一条剖面地重建。比如机器人按压时只关心接触中心区域的高度轮廓或者我们用滚筒式扫描接触痕迹每列像素都是一条独立的轮廓线。对每条剖面线边界条件通常是一端固定、一端自由。固定端可以取接触区域的起点自由端是远离接触的位置。这种混合边界的离散二阶微分矩阵特征函数正好是 DST-II 的基函数也就是sin(π * (n 0.5) * (k 1) / N), n, k 0, 1, ..., N-1这个公式不要被吓到它可以理解为在像素中心位置采样、但边界上满足一硬一软条件的正弦基。如果你实际做过数值验证就会发现把离散拉普拉斯矩阵的特征向量画出来它们和 DST-II 的基函数逐点重合。而 Gradient 和散度在 staggered grid 上的定义方式也天然让 DST-II 成为对界最准的那个算子。二维整面重建当然也可以用 DCT-II 一类的方法但如果要处理的是大量一维剖面DST-II 是最顺手的方案。用 DST-II 的好处是边界条件和物理模型完全一致不用人为在图像外侧补一堆假数据来凑 Neumann 边界。3. 实现流程从 GelSight 图像到 height-map3.1 图像采集与梯度场提取在实际系统里GelSight 的原始输出是一张彩色图像。以常见的 RGB 三光源配置为例红、绿、蓝三个通道各自对应一个方向的照射响应。拿到图像后的第一步不是直接算梯度而是先做预处理。我习惯先把三个通道拆开。这个操作在不同视觉库里有不同叫法在 OpenCV 里是split在类 Halcon 的更底层操作里可以理解为一个通道切片channel slice算子把三通道图按波段切成三个单通道图。拆完通道后要做一次亮度归一化把白平衡偏差去掉。凝胶表面反射涂层如果不是完全均匀图像里会出现低频明暗不均这一步能明显缓解后续梯度场的系统性倾斜。接下来用一阶微分算子提取梯度。最常用的是 Sobel 算子它的原理是对图像做带权重的差分同时起平滑作用。X 方向的 Sobel 核提取水平梯度Y 方向的 Sobel 核提取垂直梯度。对 GelSight 这种表面纹理较弱的图像Sobel 十分合适因为它对噪声有一定抑制不会像单纯的前向差分那样把像素噪声全部放大。import cv2 import numpy as np img cv2.imread(gelsight_sample.png) b, g, r cv2.split(img) # 用三个通道的亮度关系换算出表面梯度分量 # 具体换算矩阵由 LED 方向标定得到这里用近似系数示意 gx 0.7 * (r.astype(np.float32) - g.astype(np.float32)) gy 0.7 * (g.astype(np.float32) - b.astype(np.float32))这段代码是一个高度简化的示意。严格的做法是先在标定块上采集已知平面和若干已知曲率目标的光照响应建立“通道亮度到局部倾斜角度”的查找表或线性映射矩阵。做一次完整标定之后GelSight 的测量精度可以到微米级。3.2 构建散度项与边界条件有了 gx 和 gy 之后下一步是把它们合成泊松方程的右端项 b div(g)。离散散度用的是差分格式。为了和 DST-II 的 staggered grid 配合我通常把梯度定义在半像素位置散度则用相邻梯度之差来近似div(g)[i][j] (gx[i][j] - gx[i-1][j]) (gy[i][j] - gy[i-1][j])边界处则根据具体模型做处理。如果是逐列剖面重建每一列的起点高度直接固定为 0这是 Dirichlet 条件剖面末端认为法向倾角为 0这是 Neumann 条件。这种一边固定、一边自由的组合正好和 DST-II 的特征函数匹配。如果做整面二维重建边界处理会更繁琐。我的习惯是先把图像边缘外扩 4 个像素外扩方式用镜像对称然后再算散度。这样能显著减少重建出的高度图边缘出现的“碗状”扭曲也就是边界伪影。3.3 DST-II 快速求解步骤与代码级示意求解阶段用 DST-II 变换。Python 生态里没有像 FFT 那样开箱即用的标准 DST-II 封装但可以用 FFTW 的 RODFT10 变换类型或者自己在 NumPy 里手写一个基于 FFT 的快速实现。一维剖面重建的核心逻辑如下def dst2_1d(x, axis0): # DST-II 可以借助 FFT 实现这里为可读性直接用矩阵乘法示意 N x.shape[axis] n np.arange(N) k np.arange(N).reshape(-1, 1) basis np.sin(np.pi * (n 0.5) * (k 1) / N) return basis x # 假设 gx 是某一行剖面的 x 方向梯度 N gx.shape[0] b np.zeros(N) b[1:-1] gx[1:] - gx[:-1] # 内部散度 b[0] gx[0] # 起点固定 h[0] 0散度单独处理 # 变换到频域 b_hat dst2_1d(b) # 频率域逐点除以特征值 # 特征值与离散拉普拉斯矩阵和边界条件绑定实际建议数值验证后使用 k np.arange(N) lam 2 * np.cos(np.pi * (k 1) / N) - 2 lam[0] 1.0 # 避免除零直流分量在固定起点场景下无意义 h_hat b_hat / lam # 逆变换DST-III 或 DST-II 的逆注意归一化系数 h np.zeros(N) # 实际工程中建议用 FFTW 的 RODFT10 和 RODFT01 配对 # 这里的矩阵版本省略了归一化作示意用实际工程项目里我不会建议手写这个变换。成熟的方案是FFTW 的 RODFT10 做正变换RODFT01 做逆变换FFTW 里这两个变换类型就是标准 DST-II 和它的逆。如果你用的是基于 FFTW 封装的 Python 库例如pyfftw可以直接调用。配合多线程1080p 图像的逐列剖面重建能在几十毫秒内完成。二维整面重建也可以照搬这个思路先在 x 方向做一次 DST-II再在 y 方向做一次 DST-II散度项用二维差分构造最后在频域做一次二维逐元素除法。边界条件需要和变换类型严格匹配实际操作时我用一维剖面版本做了一个 small test再扩展到二维。3.4 可视化与灰度拉伸重建出的 height-map 是浮点数据数值范围往往非常小。比如一个深度 50 微米的微小压痕它的高度值在 -2.5e-5 到 2.5e-5 之间浮动直接换算成 8 位图就是一片漆黑。这时候必须要做灰度拉伸。在 Halcon 里我经常用灰度值拉伸算子配合min_max_gray统计图像中的极值再把整个深度范围线性映射到 0 到 255。OpenCV 的话可以用normalize配合 percentile 截断把 0.5% 到 99.5% 的数值映射到全灰度范围防止个别极端噪声点拖垮整体对比度。除了可视化灰度拉伸还能辅助后续的特征提取。比如做缺陷检测时经过拉伸后微小划痕变成清晰的亮暗纹理再用阈值分割或者边缘检测就方便多了。4. 实操避坑与常见问题4.1 边界伪影与条带噪声第一个坑是边界伪影。用 DCT/DST 类变换求解泊松方程时边界条件如果和变换的特征结构不匹配会看到重建结果在图像边缘出现明显的“拱起”或“下陷”。尤其是二维重建时如果四周都当自由边界处理但实际接触区域外 GelSight 的梯度场并不为零重建出的高度图边缘就会像碗边一样翘起来。我常用的解决办法是外部扩展。在进入散度计算前先把原始梯度图向外侧镜像扩展 8 到 16 个像素求解完成后再把扩展部分裁切掉。这样外扩区域承担了边界效应目标区域的边缘就干净很多。条带噪声则通常来自光源不均匀或者反射涂层缺陷处理思路是先对 gx、gy 做沿扫描方向的中值滤波再用灰度拉伸增强对比。4.2 标定偏差与光源不均GelSight 的梯度场是从三通道亮度换算过来的换算矩阵的准确度直接决定 height-map 的精度。我在早期试验时发现重建出的球面轮廓总是带一个莫名其妙的“锥形偏置”排查了很久才发现是红色 LED 的亮度衰减比蓝色 LED 快导致换算系数在高亮度区域出现了非线性偏移。解决方法是做多点标定用一堆已知曲率半径的标准球或标准斜块在不同亮度区间采集样本拟合分段线性或多项式映射。另外光源老化会导致标定参数漂移所以项目上线后最好定期重新标定。如果环境光干扰严重可以考虑在传感器外面加遮光罩或者在算法里加一个暗帧扣除盖上不透光盖板拍一张全黑图后续每帧都减去这张暗帧能有效消除固定模式噪声。4.3 高度场漂移与倾斜修正即使用了泊松方程重建出的高度场有时还是会带整体倾斜项。原因是传感器本身的相机光轴和凝胶表面法向不完全平行导致梯度场里混入了一个恒定偏置。解决方式是在重建后做一个平面拟合把拟合得到的倾斜平面从高度场里减去def remove_tilt(h): H, W h.shape y, x np.mgrid[0:H, 0:W] A np.stack([x.ravel(), y.ravel(), np.ones_like(x.ravel())], axis1) coef, _, _, _ np.linalg.lstsq(A, h.ravel(), rcondNone) tilt (coef[0] * x coef[1] * y coef[2]).reshape(H, W) return h - tilt这个操作在离线分析和实时处理中都很常见。如果高度场还伴随低频波浪一般是因为凝胶表面自身不是绝对平面这时候可以用高通滤波或者多帧平均背景相减来处理。4.4 性能与硬件实现算子优化方向整套 GelSight 重建流程实际是由一连串算子组成的拆分通道、灰度归一化、Sobel 梯度提取、散度计算、DST-II 正变换、频域除法、DST-II 逆变换、灰度拉伸。在桌面处理器上跑1080p 分辨率下面向实时毫无压力但一旦要把这套能力部署到机器人夹爪、嵌入式视觉模块里硬件压力马上就来了。大量的算子顺序执行每一遍都涉及对整幅图像的内存读写瓶颈通常不在算力而在内存带宽。所以优化第一优先级是算子融合把 Sobel 梯度提取和散度计算合并成一次卷积省一趟内存读写把拆分通道和白平衡合并到相机原始数据到 RGB 的 ISP 流程里频域除法和逆变换的前半部分可以在同一个循环里完成。我做过一个嵌入式平台上的 DST-II 实现发现手写蝶形运算时主要性能杀手是三角函数查表带来的 cache miss。把角度预计算成 LUT 后整体算子快了近 30%。如果再往底层走可以用定点整数近似替代浮点运算不过要注意动态范围建议先仿真验证误差。DST-II 本身是一个高度规整的移位加乘结构很适合在 GPU 或 DSP 上用 SIMD 指令加速。这也是现在算子开发工作里很常见的优化方向从数学形态不变的算法等价变换开始再针对目标硬件做手脚。问题可能原因排查方向重建结果边缘翘起边界条件与变换不匹配扩展边界再做散度或改用混合边界模型高度图整体倾斜相机光轴与凝胶面不垂直重建后平面拟合去倾斜或重新标定出现横竖条带纹路光源不均、反射涂层磨损暗帧扣除、梯度图滤波、定期重新标定局部突变噪声反射涂层划伤或污渍中值滤波、限制单点误差权重求解非常慢通用求解器直接搬用换成 DST-II 快速解法注意变换归一化高度值整体偏小/偏大标定系数不准用已知高度标准块回归校验最后再分享一个小技巧。DST-II 的特征值在不同库里的排列顺序和归一化系数差异很大第一次写代码时不要直接套论文公式先用一个已知的小矩阵做数值验证随手构造一个小尺寸测试高度图算出它的梯度场再用你的 DST-II 求解器重建回来看误差是不是在浮点精度范围内。我每次换平台重写这个算法都会先跑一遍这个小测试能省掉后面排查一堆莫名其妙问题的功夫。这套方法后续还可以扩展比如用在柔性触觉阵列的非接触标定、压痕深度统计、或者机器人滑移检测里的微几何变化追踪上核心都是同一套“梯度场到高度场”的快速求解逻辑。