第一次拿到 GelSight 传感器拍回来的原始图像我的第一反应是这玩意儿看着就像一张高光塑料膜的照片红色、绿色、蓝色反光混在一起完全不像能直接用来做触觉反馈的数据。但真正懂行的人知道这一张图里其实编码着接触表面的微观形貌——凝胶被压下去多少、哪个方向变形、表面是凸还是凹全都藏在颜色变化里。问题是怎么把这些颜色变成一张真正的 height-map高度图/表面形貌图答案绕不开两个东西一个是泊松方程另一个是 DST-II 快速求解。这篇内容写给两类人一类是做机器人触觉感知、想把 GelSight 数据变成可量化 surface 模型的工程师另一类是做三维重建、光度立体、shape-from-shading 时被“梯度积分不一致”折磨过的同学。我会从数学原理讲到可直接复用的 Python 实现再把我踩过的真实数据坑一个个摆出来。1. 为什么重建 height-map 本质上是解泊松方程1.1 GelSight 到底给了我们什么从三色光照到表层法向量GelSight 的核心设计很有意思一块弹性凝胶表面镀了反射涂层底下放着彩色相机四周用红、绿、蓝三种颜色的 LED 从不同角度照射。当传感器接触物体时凝胶表面发生形变局部坡度会改变反射光进入相机的方向——颜色就变成了表面梯度的编码。红色通道对应一个方向的倾斜绿色对应另一个方向蓝色再对应一个方向。之后再通过标定好的查找表或者光度立体算法把 RGB 图像转换成每个像素的表面法向量 n (nx, ny, nz)。也就是说GelSight 输出的并不是高度本身而是梯度的信息。这有点像我们看一座山肉眼能看出“这里陡、那里缓”但很难直接量出“山顶到底比山脚高多少米”。要得到高度场必须把无数个局部斜率“拼”回一个连续曲面。这里面最关键的一步就是把法向量先转成梯度场。公式很直接p -nx / nz, q -ny / nz其中 p 表示高度在 x 方向的偏导数 ∂h/∂xq 表示 y 方向的偏导数 ∂h/∂y。这个负号经常有人搞反我建议你自己推一遍如果表面高度在某个方向上升法向量会往相反方向偏所以梯度是法向量水平分量的负比例。1.2 直接沿路径积分会把误差累积成裂纹很多人拿到 p 和 q 之后第一反应是随便找个起点沿行方向做累加积分h(x1, y) ≈ h(x, y) p(x, y)这种朴素做法在无噪声的理想数据上能跑但真实 GelSight 数据几乎必崩。原因很简单p 和 q 是从光度信息估计出来的天然带噪声还夹杂着标定误差和光照不均带来的系统偏差。一个理想的梯度场必须满足可积性条件也就是混合偏导数相等 ∂p/∂y ∂q/∂x。但实测梯度场的旋度根本不是零。一旦旋度不为零沿不同路径积分就会得到不同的结果。你从左边往右积和从右边往左积回到同一个点时高度差一截画出来就是一条条“裂纹”或者“台阶”。这就是经典的 shape-from-shading 里的积分不一致问题。我在刚做这个项目时试过先对梯度做中值滤波再积分结果只是把裂纹变模糊了并没有消除。后来才明白这个问题不能靠“修梯度”解决而是要把问题重定义为一个全局优化问题。1.3 泊松形式把不可积梯度场拉回可积空间既然找不到一个完美满足观测梯度的曲面那就找一个“误差最小”的曲面。定义目标函数E(h) Σ ( ∂h/∂x - p )² ( ∂h/∂y - q )²这个式子的意思是我们要找一个高度场 h让它求出来的梯度尽量接近观测到的 p 和 q。对所有像素求和差异以平方计。用变分法对 h 求导并令导数为零就得到著名的泊松方程Δh ∂p/∂x ∂q/∂y等号左边是 h 的拉普拉斯算子二阶差分右边是观测梯度场的散度。这一步的意义非常深刻我们不再逐点积分而是把“全场面上的拼接误差”压成一个整体最小二乘问题。只要解出这个线性方程就得到了一致性最好的高度场。到这里问题从“怎么积分”变成了“怎么快速求解一个大规模的稀疏线性系统”。对一个 512×512 的图来说未知量有 26 万个直接构造矩阵求解的内存和耗时都不可接受所以必须要用变换域方法。2. DST-II 凭什么能快拉普拉斯算子的谱对角化2.1 一维差分矩阵和正弦基的第一个巧合先看一维情况。假设有一个长度为 N 的列向量 u离散二阶差分算子在内部点可以写成矩阵形式 A它的对角线是 -2相邻对角线是 1。这个矩阵长得很规矩是三对角对称矩阵。学过线性代数的都知道对称矩阵一定可以对角化而且它的特征向量有非常漂亮的形式v_k(i) sin(π (k1) (i 0.5) / N)也就是说这个矩阵的特征向量就是离散正弦序列。为什么会这样可以从差分方程直接验证把 v_k 代入二阶差分利用三角恒等式 sin(AB) sin(A-B) 2cos(B) sin(A)可以得到A v_k ( 2cos(π (k1)/N) - 2 ) v_k这里产生了一个奇妙的巧合差分矩阵的本征分解恰好就是离散正弦变换DST本身。这意味着如果我们把输入信号变换到 DST 域矩阵 A 就变成了一个对角矩阵求解线性系统变成一次逐元素的除法然后再逆变换回来。整个过程从解一个大规模稀疏矩阵变成了三次 O(N log N) 的快速变换。2.2 DST-I / DST-II 的边界语义差别很多资料会笼统地说“解 Dirichlet 边界条件的泊松方程用 DST”但实际使用时你会遇到一个困惑scipy.fft.dst 有 type1、type2 之分到底用哪个它们对应的是两种略有不同的边界采样约定。以长度为 N 的序列为例DST-Itype1基函数是 sin(π (k1) (i1) / (N1))它对应的是“网格两端点都固定为零”的经典 Dirichlet 边界。MATLAB 老代码和许多图像处理教材里的 dst 默认是这个。DST-IItype2基函数是 sin(π (k1) (i0.5) / N)它隐含的边界是“网格左右两侧半像素处值为零”也就是半格点奇延拓。scipy 的 dst 默认是 type2。两类变换在对角化拉普拉斯算子时都能用特征值公式略有差别。如果用 DST-II长度 N 对应的一维特征值是λ_k 2cos(π (k1)/N) - 2如果写成 DST-I通常把内点数量记为 N则特征值是 2cos(π (k1)/(N1)) - 2。我见过很多人把这两个公式混着用最后解出来的图边界出现奇怪波纹却不知道原因十有八九就是库的 type 编号和特征值公式没对上。所以我的建议是不要死记编号要跑合成测试。构造一个已知解析解的曲面把泊松右端项算出来用你的实现还原如果误差在 1e-10 级别说明公式对了如果边界有明显偏差先查特征值分母是 N 还是 N1。2.3 二维扩展与特征值公式二维情况实际上是一维的直积推广。离散拉普拉斯算子可以写成 A_x ⊗ I I ⊗ A_y也就是先在 x 方向做二阶差分再在 y 方向做二阶差分。由于 DST 对这种张量积结构天然兼容二维 DST 可以拆成“先沿每一列做变换再沿每一行做变换”。在 DST-II 的约定下对于 m×n 的图像二维特征值就是两个一维特征值相加λ_{i,j} 2cos(π i / m) 2cos(π j / n) - 4注意这里 i 从 1 到 mj 从 1 到 n。这个公式中 i、j 都是“1-indexed”也就是说最小频率不是 0。这正好避免了传统 FFT 解法里零频率带来的除零问题——Dirichlet 边界下根本没有常数特征向量绝对高度本来就是不可恢复的。整个快速求解流程可以写成一行伪代码h IDST2( DST2(f) / λ )其中 f 是泊松方程右端项也就是散度场 ∂p/∂x ∂q/∂y。所有计算都在谱域完成复杂度 O(m n log(m n))。2.4 和 FFT / DCT 的对比边界条件决定选型既然有 FFT为什么不用 FFT因为 FFT 隐含的是周期性边界条件图像的左边会“卷绕”到右边。对于自然图像或者 GelSight 接触区域左右边缘的高度并不连续用 FFT 会在边界造成严重的接缝伪影。DCT 也很常见它对应 Neumann 边界条件也就是边界处梯度为零适合处理“曲面在边缘自然延展”的假设。很多做光度立体的经典代码用的就是 DCT 版本Frankot-Chellappa 的后续改进里也常出现 DCT。DST 对应的是 Dirichlet 边界条件边界处高度固定为零。对 GelSight 来说这个假设其实很自然接触区域外面就是未变形的凝胶高度就是零。所以我个人优先选择 DST 系列。如果实际数据的接触区域没有完全贴到图像边缘DST 和 DCT 的结果差别非常小主要影响集中在最外圈。三种变换的选择本质上是问你图像边界以外的世界是什么样的是循环重复FFT、镜面延续DCT还是一无所有DST想清楚这个选型就不会错。3. 手写一个可直接复用的快速泊松求解器3.1 法向量到梯度场的转换细节先处理输入。假设你已经通过光度立体得到了三个通道的法向量分量 nx、ny、nz都是和原图同尺寸的数组。第一步是求梯度场import numpy as np # nx, ny, nz: 法向量分量float32/float64nz 理论上应该 0 def normals_to_grad(nx, ny, nz, eps1e-8): valid nz eps p np.zeros_like(nx) q np.zeros_like(ny) p[valid] -nx[valid] / (nz[valid] eps) q[valid] -ny[valid] / (nz[valid] eps) return p, q这里把非法像素nz 接近 0的梯度直接置 0对应“未接触区域没有形变”的物理事实。如果你发现真实数据的接触区域外法向量噪声很大建议先做一个 mask 再置 0后面第四节会详细说。3.2 核心函数两行 DST 一次点除接下来是主角。下面的函数接收泊松方程右端项 f散度场返回重建的高度场 h。我用的是 scipy.fft 的 dst正变换用 type2逆变换用 type3并且都开启 normortho保证正交归一化。from scipy.fft import dst def fast_poisson_solver(f): m, n f.shape # 正向二维 DST-II F dst(f, type2, normortho, axis0) F dst(F, type2, normortho, axis1) # DST-II 约定下的拉普拉斯特征值谱 kx np.pi * np.arange(1, m 1)[:, None] / m ky np.pi * np.arange(1, n 1)[None, :] / n denom 2 * np.cos(kx) 2 * np.cos(ky) - 4 # 谱域除法 H F / denom # 逆向二维 DST-III h dst(H, type3, normortho, axis0) h dst(h, type3, normortho, axis1) return h就这么简单。正变换两行、逆变换两行、中间一次除法。整个流程没有任何迭代也没有构造任何稀疏矩阵。构造右端项 f 时要注意散度要用中心差分而不是 np.gradient 的默认边界处理。我一般这样写def divergence(px, py): dpx np.zeros_like(px) dpy np.zeros_like(py) dpx[1:-1, :] (px[2:, :] - px[:-2, :]) / 2 dpy[:, 1:-1] (py[:, 2:] - py[:, :-2]) / 2 return dpx dpy把 p、q 传进 divergence 后得到 f再把 f 交给 fast_poisson_solver就得到重建的高度场。整个 pipeline 只有十几行代码。3.3 合成曲面验证误差要到 1e-10 级别任何变换域求解器都要做同一件事用一个已知解析解的曲面验证实现是否正确。我用一个“余弦帽”曲面内部余弦形状外部为零边界处值和一阶导数都连续非常适合 Dirichlet 边界。def cosine_cap(n, radius, cx, cy): yy, xx np.mgrid[0:n, 0:n] r np.sqrt((xx - cx) ** 2 (yy - cy) ** 2) z np.where(r radius, 0.5 * (1 np.cos(np.pi * r / radius)), 0.0) return z验证流程z_true cosine_cap(128, 50, 63.5, 63.5) # 用中心差分近似解析梯度 p np.zeros_like(z_true) q np.zeros_like(z_true) p[1:-1, :] (z_true[2:, :] - z_true[:-2, :]) / 2 q[:, 1:-1] (z_true[:, 2:] - z_true[:, :-2]) / 2 f divergence(p, q) z_rec fast_poisson_solver(f) print(np.max(np.abs(z_rec - z_true))) # 输出通常是 1e-12 这个量级如果输出在 1e-12 附近说明你的 DST 方向、特征值分母、逆变换 type 全都对上了。如果出来的误差是 1e-3 甚至更大先检查特征值公式里分母是 m 还是 m1再检查逆变换是不是用 type3。3.4 从梯度场到 height-map 的完整流水线把上面的模块串起来一个完整的 GelSight height-map 重建函数长这样def gelsight_heightmap(nx, ny, nz): p, q normals_to_grad(nx, ny, nz) f divergence(p, q) h fast_poisson_solver(f) # 去掉常数偏移让最小值归零方便可视化 h - h.min() return h生成 h 之后可以用 matplotlib 的 imshow 加 colormap 显示或者保存成 16-bit PNG 做后续处理。注意这里的高度单位是“像素网格单位”并不是实际毫米要换算成物理尺寸需要标定这个我在下一节讲。4. 拿到真实 GelSight 数据后最容易翻车的四个地方4.1 全局倾斜伪影梯度均值不等于零合成数据永远是完美的真实 GelSight 数据第一个坑就是重建出来的表面经常是一个“锅盖”形状——中间鼓起、四周下沉或者整体往一个方向倾斜。原因在于光度立体标定不完美LED 亮度不对称导致估计出的法向量场有一个全局偏置p 和 q 的整体均值不为零。这个偏置在频域里表现为极低频分量泊松方程会忠实把它还原成大尺度的倾斜面。解决办法是在重建完成后做一个平面拟合把高度场投影到去除平面def remove_plane(h): m, n h.shape yy, xx np.mgrid[0:m, 0:n] A np.stack([xx.ravel(), yy.ravel(), np.ones_like(xx.ravel())], axis1) coeff, _, _, _ np.linalg.lstsq(A, h.ravel(), rcondNone) plane (coeff[0] * xx coeff[1] * yy coeff[2]).reshape(m, n) return h - plane我实测过这个操作能去掉绝大多数锅盖伪影。如果是更复杂的光照不均可能还需要用高通滤波器或者先对梯度场做加权均值零化。核心思想是形貌信息集中在局部梯度里全局倾斜是标定误差的低频污染。4.2 非矩形接触区域mask 里的泊松是个坑GelSight 的视场通常是圆形或者圆角矩形而 DST 要求整个矩阵都参与变换。如果你直接把 mask 外的法向量置 0那么 mask 边缘会出现从有效梯度到零的突变。这个突变反映到散度场里就像一个巨大的“电荷”堆积重建出来的高度场在接触区域外会有一圈明显的“围墙”。我曾经在这个问题上卡了两天。后来总结出三个层次的应对方案对 mask 做膨胀几像素然后再置零外部梯度让过渡区平滑一些在 mask 边缘加一个软过渡带比如用高斯模糊后的 mask 做加权把边界附近的梯度逐渐压到零如果精度要求很高那就别再依赖纯 DST 了改用带 mask 的迭代式泊松求解这类方法能显式指定哪些像素参与计算。DST 的快速性牺牲的是边界灵活性它只能处理矩形域。对大多数 GelSight 接触检测场景来说方案 1 和 2 足够而如果你要做精密计量建议在 DST 结果基础上再用共轭梯度法精修几轮。4.3 高频噪声放大除以零附近要加正则项泊松方程求解在频域里就是一个除法F / denom。denom 在高频区域约等于 -4低频谱分量则接近 0。这意味着低频信号会被放大得厉害这是它能把表面形状恢复出来的原因但同时也意味着非常微弱的低频噪声会变成大尺度波纹。更隐蔽的问题是高频噪声GelSight 图像里每个像素的法向量估计都带噪声这些噪声经过散度计算后已经比较高频除以 denom 后会形成看起来像“砂纸”一样的颗粒纹理。如果你觉得重建出来的表面太毛糙可以给 denom 加一个小的阻尼项eps 1e-4 H F / (denom eps)注意 eps 不能太大否则会把真实形貌的低频能量也一并压掉。实际项目中我会先跑一遍无正则版本看看噪声水平再决定 eps 的量级。也可以做成自适应的denom 绝对值大的高频位置几乎不受影响denom 接近零的低频位置加一点约束。4.4 绝对高度和单位height-map 是相对尺度最后是个经常被忽略的物理问题泊松方程只能确定曲面形状不能确定绝对高度。二维 Dirichlet 边界下不存在零频分量重建出来的 h 总是一个“相对值”——你把它整体平移 100 个单位梯度完全不变误差函数也不变。所以看到重建结果后不要问“这个像素的高度绝对值是多少毫米”而要问“这个凸起相对背景高多少”。要做绝对标定必须找已知深度的参考物比如一个半径确定的钢珠压出已知曲率半径的凹痕然后用重建结果反推像素高度到微米的换算系数。GelSight 社区把这一步通常叫“深度标定”每个传感器都需要单独做因为凝胶厚度和弹性模量会影响响应曲线。5. 从快速泊松再往外走一步5.1 与 Frankot-Chellappa 的关系做光度立体的朋友应该听过 Frankot-ChellappaFC算法它在频域把梯度场投影到可积函数空间。FC 的核心目标函数和泊松重建几乎完全一样最小化重建梯度与观测梯度差的平方和。实际上在规则网格和特定边界条件下FC 的闭式解和泊松方程解就是同一个东西。区别主要在两点一是 FC 传统上用 FFT 基隐含周期边界图像边缘会有明显卷绕二是在频域处理中FC 可以很方便地加入权重或者用不同的积分基函数。如果你的工程代码里已经有 FC 的实现想切到 DST 版本只需要把特征的核从 exp(-iux) 换成 sin 基然后除以对应的拉普拉斯特征值本质就是我上面写的 fast_poisson_solver。所以这篇文章给的求解器本质上是一个边界条件更合理的 FC 变体。5.2 什么时候该放弃 DST去用迭代法DST 不是万能的。当出现以下情况时建议直接切换到预条件共轭梯度PCG或者多重网格法接触区域的 mask 很复杂不能近似成矩形法向量的有效区域只有图像的一小部分大部分像素是无效背景需要对泊松问题加非标准的权重例如对高梯度区域置信度更高需要实时处理但实时线程上没有现成的 DST 库。作为参考在我自己的机器上对 512×512 的图像DST 快速求解只需要大约 15 毫秒而 scipy 的稀疏 Cholesky 直接解法大约需要 1 秒到 2 秒PCG 要好一些但也要几十毫秒到上百毫秒具体看迭代次数。DST 在规则网格上的优势是碾压级的前提是边界条件适配。5.3 重建结果能做什么有了干净的 height-map后面可以做的事情就多了可以直接做表面纹理分类比如区分丝绸、砂纸、牛仔布可以计算曲率图用于检测接触区域的凸起和凹陷特征可以提取接触区域的面积、质心、接触轮廓反馈给机械臂做力控还能做滑移检测——当物体在凝胶表面滑动时重建高度图会出现明显的纹理运动场这个信号比单纯看原始 RGB 图像更稳定。我自己做的下游是抓取稳定性判断在机械手抓取过程中实时重建 height-map计算接触区域的曲率变化一旦检测到接触轮廓发生非连续跳变就判定物体有滑动趋势立刻加大夹持力。这个闭环系统里快速泊松求解是真正卡实时性的关键模块DST 版本让我把整个重建周期从几百毫秒压到了几十毫秒以内。最后分享一个我自己调试时的小技巧无论你用 DST 还是 DCT拿到代码第一件事不是接真机 GelSight 数据而是用一个解析已知的余弦帽跑通打印误差。如果误差不是 1e-10 级别多半是特征值公式里的分母写成了 N1 而实际应是 N或者逆变换用错了 type。这类算子级错误从重建图像的轮廓上很难一眼看出来但合成测试一秒就能定位。等你把合成测试跑通再接真实数据你就知道后面那些锅盖、围墙、砂纸纹理全都是数据本身的坑跟求解器没关系了。