TPS薄板样条变换：从物理直觉到图像配准实战

做图像配准、人脸形变或者点云处理的朋友，大概率都撞见过“TPS”这三个字母。更巧的是，如果你同时做性能测试，会发现性能圈子里的TPS是Transactions Per Second（每秒事务数），一度让我在查资料时怀疑人生。这篇文章要聊的，是数学和图像处理领域的那个TPS——Thin Plate Spline，薄板样条变换。简单说，它解决的是这样一个问题：给你一组控制点在变形前和变形后的对应位置，能不能构造一个平滑的映射，把整个空间（或整张图）都自然地“揉”过去，既精确对齐已知点，又不让形变显得生硬。适合想入门图像配准、需要做数据增强变形、或者被一堆数学公式劝退过的朋友阅读。

我从最早碰TPS到现在，踩过不少坑，也踩出了一些心得。这篇文章不会只堆公式，而是从物理直觉讲到数学原理，再讲实际代码里那些没人明说的细节，争取让你看完能直接上手用。

1. TPS到底在干什么？一个物理直觉先立起来

1.1 想象一块会被“按”弯的金属薄板

TPS的全称是Thin Plate Spline，翻译过来就是“薄板样条”。这个名字其实已经把核心思想说得很直白了：想象一块又薄又平的金属板，你在上面选若干个点，用手把这些点往某些方向按下去或者提起来，金属板会怎么样？它会形成一个光滑的、连续的整体曲面，而且这个曲面的弯曲是“尽量省力”的——金属板会自然地选择弯曲能量最小的那种形态，不会无缘无故地剧烈褶皱。

这个物理类比，就是TPS的核心逻辑。在数学上，我们拿到的是若干控制点对：原始坐标((x_i, y_i))和目标坐标((x_i', y_i'))，或者更一般地，某个标量场在控制点处的取值。我们要找一个函数(f(x, y))，让它满足所有控制点的约束，同时让整个曲面的“弯曲能量”最小。它不追求用一条直线去拟合所有点，也不追求简单的全局缩放旋转，而是允许不同区域有不同尺度的形变，就像金属板在多个钉子约束下形成的自然曲面。这是TPS和传统线性变换最本质的区别：后者是刚性的、全局的；前者是柔性的、局部的。

这个直觉能帮你理解后面所有公式：为什么TPS能同时处理全局旋转和局部扭曲，为什么它在控制点之外的地方也能平滑过渡，为什么控制点太少会形变过度、控制点太密又会病态震荡——本质上，都是因为这块“虚拟金属板”在物理规律下做出的选择。

1.2 和仿射变换、透视变换放在一起看，差异就很清晰

很多刚接触图像几何变换的同学会困惑：已经有仿射变换（Affine）和透视变换（Perspective/Projective），为什么还要TPS？这三种变换的应用场景和自由度完全不同，放在一起对比就一目了然。

举一个直观的例子：人脸从面无表情到咧嘴笑，嘴角附近的皮肤被拉伸，眼角会产生细微的纹理变化，颧骨周围甚至会有轻微的鼓起和凹陷。这种形变是高度局部的——嘴部区域变化大，额头区域几乎不动。用仿射变换处理，整张脸只会被统一拉伸或旋转，效果就是一张“僵硬的贴图”；而TPS可以根据每个关键点的位移自由调节局部区域，额头保持不动、嘴角附近拉伸、脸颊轻微隆起，最后的输出才符合物理直觉。这也是为什么很多老式基于关键点的人脸变形Demo，背后都是TPS在做空间映射。

还有一点值得注意：仿射变换是TPS的特例。如果所有控制点的位移都刚好满足同一个全局仿射关系，那么TPS计算出来的非刚性项权重会趋近于零，最终退化为一个纯仿射变换。反过来，TPS总能自动在全局线性成分和局部扭曲成分之间找到一个平衡。所以很多时候，你可以把TPS当作一个“升级版的仿射变换”，它天然兼容全局和局部的混合形变。

2. 数学原理：从物理模型到可计算的方程

2.1 弯曲能量函数：什么才算“平滑”

理解了物理直觉之后，就该把“金属板自然弯曲”翻译成数学语言了。物理上，一块薄板的弯曲能量和它的曲率有关。在二维函数(f(x, y))的语境下，我们用所有二阶偏导数的平方和来度量曲面在某一点的弯曲程度，对整个平面做积分，就得到总弯曲能量：

[ E[f] = \iint \left[ \left( \frac{\partial^2 f}{\partial x^2} \right)^2 + 2 \left( \frac{\partial^2 f}{\partial x \partial y} \right)^2 + \left( \frac{\partial^2 f}{\partial y^2} \right)^2 \right] dx,dy ]

这个公式的直觉非常直接：一阶导数代表曲面“倾斜”了多少，倾斜本身并不消耗能量——一块斜着的平板依然是平的；二阶导数代表的才是“弯曲”——它是曲率变化的度量。如果曲面在这一点的二阶导数大，说明它弯得很急；处处二阶导数为零，就是标准的平面。TPS求解的目标，就是在满足控制点约束的众多函数里，找到让这个能量值最小的那一个。

这里的关键是“在所有满足约束的函数里找最光滑的那个”。生活里可以这样类比：你有一组散点，让你用一根软铁丝做一条穿过所有点的曲线。你可以把铁丝拧成一个极度曲折的形状，也可以用一个大圆弧舒缓地穿过所有点，两者都满足“穿过所有点”的约束，但后者明显更“自然”。TPS做的，就是在无穷多个候选函数中，自动选出那个最“舒缓”的。数学上可以证明，这个最优化问题的解是存在且唯一的，而且具有非常优美的解析形式——这正是下面要讲的径向基函数展开。

有意思的是，弯曲能量函数这一定义天然依赖于维度。上面这个公式是二维函数(f(x, y))的版本，如果是三维空间中的曲面，能量定义会涉及更多偏导数项；如果是处理一维曲线，能量会变成(\int (f'')^2 dx)。这也是为什么薄板样条在不同维度下有不同的径向基函数形式，千万别拿二维的公式硬套到三维配准上。

2.2 径向基函数：为什么偏偏是 (r^2 \log r^2)

有了弯曲能量最小化这个目标，下一步就是找到最优函数的表达式。这里直接给出结论（也是TPS最美妙的地方）：最小化上述弯曲能量的函数，可以写成一组径向基函数加上一个仿射项的形式：

[ f(x, y) = a_0 + a_1 x + a_2 y + \sum_{i=1}^{N} w_i , U(|\mathbf{P}_i - (x, y)|) ]

其中(U(r) = r^2 \log(r^2))，(\mathbf{P}_i)是第(i)个控制点，(w_i)是对应的权重系数，(a_0, a_1, a_2)是仿射系数。

为什么径向基函数偏偏是(r^2 \log r^2)而不是别的？这要从数学的变分法说起。求解“在约束下最小化泛函”的问题，等价于求解一个偏微分方程（欧拉-拉格朗日方程），而(U(r) = r^2 \log r^2)恰好是这个方程的基本解（格林函数）。换句话说，它天生就是“在一点施加单位冲击后，薄板在该处形成的最小能量形态”。既然每个控制点都像一个“钉子”，组合它们的形态就能表达出整块薄板在多个钉子约束下的最终形态。你完全可以把(U(r))想象成一个“基础鼓包”，每个控制点贡献一个强度为(w_i)的鼓包，而仿射项提供一个整体的倾斜台面，叠加在一起就是最终曲面。

这个基函数还有一个很实用的特性：(r=0)时函数值平滑地趋向于0，但在原点附近斜率变化很快，这意味着它能表达非常局部的剧烈形变；而在远离控制点的区域，它又不会无限发散，曲面可以平滑延展。如果是三维的TPS，基函数会变成(U(r)=r)，一维情况下是(U(r)=|r|^3)，这也再次说明“不同维度要用不同形式的径向基函数”这件事不能想当然。

至于为什么要加仿射项(a_0 + a_1 x + a_2 y)，原因也很直观：径向基函数只负责“局部扭曲”，但整块板整体平移、旋转、缩放这些全局运动，必须由仿射项来承担。如果没有仿射项，哪怕所有控制点都只是在做一个整体的平移旋转，TPS也可能被迫用大量径向基函数去硬凑这种全局运动，结果就是严重的数值不稳定和无法预测的边缘外推行为。两个部分各司其职，是TPS设计上的巧妙之处。

2.3 为什么还需要三个“平衡约束”

构造线性系统时，光有N个控制点约束还不够，学术上还会额外加上三个约束条件：

[ \sum_{i=1}^{N} w_i = 0, \quad \sum_{i=1}^{N} w_i x_i = 0, \quad \sum_{i=1}^{N} w_i y_i = 0 ]

这三个约束的物理意义，是让径向基函数部分和仿射部分“正交解耦”。直白点说，如果没有这些约束，径向基函数系数可能会偷偷吸收掉一部分全局仿射运动——比如整体平移完全可以由一组常数权重(w_i)模拟出来，但这样做毫无必要且数值极差。加上约束后，仿射运动只能由仿射项表达，径向基函数只负责真正的局部残差，矩阵条件数也会更好。这就像拟合数据时，先把均值减去再拟合，剩下的才交给波动部分，逻辑是一样的。

有了N个插值方程加3个约束方程，我就能联立求解N+3个未知数（N个权重(w_i)加3个仿射系数）。这个方程组通常写成块矩阵的形式，左上角是控制点之间的径向基函数矩阵，右上角和左下角是控制点坐标构成的„A“字型矩阵。具体形式下一节详细拆解。

3. 怎么求出这些系数？一个线性方程组的解法

3.1 搭建完整的线性系统

实际求解TPS时，我们通常按如下方式组矩阵。假设有N个控制点(\mathbf{P}i = (x_i, y_i))，对应的目标值(z_i = f(x_i, y_i))。先计算一个(N \times N)的矩阵(\mathbf{K})，其中(\mathbf{K}{ij} = U(|\mathbf{P}_i - \mathbf{P}_j|))，也就是控制点两两之间的距离对应的径向基函数值。另建一个(N \times 3)的矩阵(\mathbf{P})，每一行是([1, x_i, y_i])。然后拼出下面的分块线性系统：

[ \begin{bmatrix} \mathbf{K} & \mathbf{P} \ \mathbf{P}^{T} & \mathbf{0}_{3\times 3} \end{bmatrix} \begin{bmatrix} \mathbf{w} \ \mathbf{a} \end{bmatrix}

\begin{bmatrix} \mathbf{z} \ \mathbf{0}_{3\times 1} \end{bmatrix} ]

其中(\mathbf{w} = [w_1, \dots, w_N]^T)，(\mathbf{a} = [a_0, a_1, a_2]^T)，(\mathbf{z} = [z_1, \dots, z_N]^T)。右下角的零矩阵对应前面说的三个平衡约束。这个矩阵的大小是((N+3) \times (N+3))，解这个方程组就得到了所有系数。

以三个控制点为例会更直观。假设三个控制点分别是((0,0))、((1,0))、((0,1))，目标值分别是0、1、2。这时候为了计算简化，我们用(U(r) = r^2 \log r^2)，那么(\mathbf{K})矩阵的对角线全是0（因为(U(0) = 0)），非对角线元素要根据点间距离算出来。比如((1,0))和((0,0))距离为1，那么(\mathbf{K}_{12} = U(1) = 1^2 \cdot \log 1^2 = 0)。再看((0,0))和((0,1))同样距离为1，也是0；((1,0))和((0,1))距离为(\sqrt{2})，对应元素是((\sqrt{2})^2 \log((\sqrt{2})^2) = 2 \log 2 \approx 1.386)。然后把(\mathbf{P})矩阵和零矩阵拼起来，求解这个6×6的方程组。实际中很少手算，但看懂这个结构对理解后续原理很有帮助。

这里需要解释一个看起来反直觉的现象：为什么矩阵对角线上的距离为0，对应的基函数值也是0？难道控制点自身对自己的影响是零？这恰恰说明径向基函数u(r)=r^2 log r^2在原点处虽然为0，但它对相邻点的“牵扯力”是非零的。每个控制点的影响是通过它与其他控制点之间的距离关系共同决定的，自身的0只是这个特定基函数在原点处的取值。最终解出来的权重会综合考虑所有点之间的相互影响，形成一个自洽的形变场。

3.2 计算复杂度与数值稳定性：为什么需要预处理

TPS看起来很美，但直接求解有一个绕不开的问题：这是个稠密线性系统，求解的复杂度是(O(N^3))。控制点数量只有几十个时完全无压力，几百个还能接受，一旦超过几千个，矩阵求逆的计算量就会变得非常可观，内存占用也会以二次方速度增长。这也是为什么在实际做图像配准时，很少有人直接对几万个特征点做TPS，往往先对控制点做降采样、分块处理，或者用迭代方法求近似解。

数值稳定性是另一个更隐蔽的坑。(\mathbf{K})矩阵的条件数可能非常大，尤其当控制点分布紧密、坐标数值很大时，矩阵会接近奇异。一个极端的例子：如果两对控制点距离极近，那么径向基函数值几乎相同，矩阵的若干行几乎线性相关，求解结果就会出现剧烈震荡。我的习惯是求解前先把所有控制点坐标归一化到([0,1])或者([-1,1])区间，也就是对坐标做一次尺度缩放，让(\mathbf{K})矩阵的元素保持在一个合理的数量级。这个处理不会改变形变结构，只是在数值上让解更稳定。

还要注意，在图像配准中，每个控制点并不只有一个标量值，而是有两个分量——x方向的位移(dx)和y方向的位移(dy)。通常的实操方式是分别对(dx)和(dy)各训练一个TPS模型，共用同一组控制点坐标，但解两个线性方程组的系数。也就是说，一个形变场需要两套((w, a))系数，分别是x方向和y方向的映射函数。千万别试图用一个标量TPS搞定二维位移，那是概念性错误。

3.3 正则化参数lambda：过拟合和欠拟合之间的天平

经典TPS的插值要求是“严格穿过所有控制点”，这在控制点本身有噪声的场景里会出问题。比如用特征点匹配算法自动提取控制点时，哪怕大部分点匹配正确，也总有少数点的坐标有像素级误差。如果强行让曲面精确穿过这些带噪声的点，形变场就会在这些错误点附近产生局部突起——听起来很厉害，实际上毫无意义。

解决办法是引入正则化项，把目标函数改成弯曲能量加一个拟合误差惩罚项：

[ E_{\text{reg}}[f] = \sum_{i=1}^N \left( f(x_i, y_i) - z_i \right)^2 + \lambda \iint \left[ \left( f_{xx} \right)^2 + 2 \left( f_{xy} \right)^2 + \left( f_{yy} \right)^2 \right] dx,dy ]

当(\lambda = 0)时，退化为严格插值；(\lambda)越大，曲面越平滑，但控制点处的拟合误差也越大。这个(\lambda)就是“平滑度”和“保真度”的天平。选(\lambda)没有一劳永逸的公式，工程上我试过有效的方法有两类：一类是留一交叉验证，每次留出一个控制点计算误差，网格搜索最优(\lambda)；另一类是经验值起步，根据坐标归一化后控制点噪声水平，从(\lambda = 0.001)到(\lambda = 1)之间做几轮实验，肉眼观察形变场的平滑度再微调。

有一点特别容易踩坑：归一化坐标之后，同样的(\lambda)值对应的平滑度会完全不同。同样一组数据，坐标范围从0到1000缩放到0到1之后，原本合适的(\lambda = 0.01)可能会变成明显的过拟合。所以如果你在某个坐标系下标定好的(\lambda)，换了一套坐标范围后不能照搬，需要重新标定。

4. 实操细节与避坑指南

4.1 控制点怎么选：分布比数量更重要

控制点的选取，直接决定TPS形变的质量。我见过不少新手一上来就提取几千个特征点，然后直接灌进TPS，结果矩阵求解慢不说，形变场在控制点密集区域出现明显的“波动”。控制点数量不是越多越好，关键是分布要合理。

第一个原则是覆盖边界。TPS的径向基函数在远离控制点的区域外推时，行为由仿射项主导，但边界附近如果控制点太少，形变场就会出现无法控制的“飞边”现象——明明是处理一块区域，边缘却莫名其妙地翘起来。我处理图像配准时，不管内部特征点有多少，都会在图像的四角和四边中点补上若干控制点，让形变场的边界行为可控。

第二个原则是分布要均匀。如果某个区域控制点特别密集，另一个区域非常稀疏，TPS解算时会在稀疏区域产生不自然的过渡。这个问题在特征点匹配中尤其常见：纹理丰富的区域特征点扎堆，纹理平坦区域一个点都没有。补救方法是做一次均匀采样，或者用类似网格化的方法把控制点的密度拉平。第三个原则是控制点对的匹配精度要尽量高，尤其是异常值必须剔除。TPS本身对匹配误差很敏感，一个偏差了10个像素的错误对应点都可能让局部区域形变失控。我通常会在跑TPS之前，先用RANSAC或人工筛选把明显错误的匹配点对清掉，这一步对结果的改善远超后续任何参数调优。

4.2 坐标归一化：数值稳定性的第一步

坐标归一化这件事，在TPS实操里的重要性怎么强调都不过分。直接用原始像素坐标（比如0到4000灰度值）计算径向基函数(U(r) = r^2 \log r^2)时，矩阵元素会出现非常夸张的量级差别，比如距离为1时函数值为0，距离为4000时函数值可能是上亿级别，这会让线性系统的条件数变得极其糟糕，求解出来的权重可能出现百万级别的数值，形变场也会剧烈震荡。

我通常在求解前做这样几步归一化：把所有控制点坐标整体减去均值，让控制点中心位于原点；再除以一个统一的比例因子，把坐标范围压缩到接近([-1, 1])。这个比例因子通常用所有控制点到中心点的最大距离。归一化之后，径向基函数的值域、矩阵的条件数都会保持在合理范围，解出来的系数也不会有吓人的量级。更重要的是，形变场在数值上是“同一种语言”的，后续做可视化、插值采样都会更稳定。

还有一个细节：归一化之后，求出来的仿射系数和权重系数都是针对归一化坐标系的。实际应用时如果要对任意像素坐标做映射，不能直接把原始坐标代进去，而是要先对这个坐标做同样的归一化变换，计算完结果后再映射回原坐标空间。这个换算很容易写错，我在代码里通常会封装一个“TPSModel”类，把归一化的均值和尺度一并保存，映射时在内部自动完成换算，避免每次手工转换。

4.3 常见问题速查表

这份速查表是我在实际项目里反复踩坑后整理的。其中“远离控制点区域失控”这个问题最容易被新手忽略，很多人做完TPS只看控制点拟合得好不好，完全不管非控制点区域，结果整张形变图一渲染出来就露馅。还有一个建议：每次运行TPS后，最好把形变场可视化出来，比如画成网格图、箭矢图或颜色图，肉眼扫一遍就基本能看出问题在哪。

5. 实际应用：从图片配准到人脸变形

5.1 图像配准与医学影像对齐

TPS在图像配准里是经典中的经典。场景通常是这样的：有两幅图像，它们拍的是同一个东西，但拍摄角度、姿态、甚至拍摄时间都不一样——比如脑部CT和MRI的断层影像、眼底在不同时间拍摄的彩照。先在两幅图上各自找到一组对应结构点（例如血管分叉点、器官轮廓的标记点），这些点就是控制点对。然后让一个图像的空间坐标经过TPS映射到另一个图像的坐标体系，再做像素重采样，就能把两幅图对齐。

医学影像领域尤其喜欢TPS，因为它能处理器官的非刚性变形。人的器官在呼吸、心跳作用下本来就会轻微移动和变形，用刚体变换或仿射变换对齐远远不够，而TPS提供的平滑变形可以很好地逼近这种软组织形变。具体操作中还有一个细节：配准不一定是“源图到目标图”的单向映射，如果两幅图之间有较大的非线性差异，通常用双向映射并做对称化处理，或者直接拟合一个平均形变场，减少配准偏差。

5.2 人脸形变与数据增强

人脸相关的应用可能是普通用户感知最强的一块。很多年前的“人脸渐变”特效、表情迁移、五官微调，核心都有关键点驱动的TPS形变。流程是先检测人脸关键点（眉毛、眼睛、鼻子、嘴巴、脸轮廓），然后手动或自动修改部分关键点位置——比如把嘴角点往上移到“微笑”位置，保持其他点不动；用TPS计算从原始关键点到新关键点的形变场，最后应用到整张脸部图像上。由于TPS本身足够平滑，面部纹理不会被撕裂，形变后的脸看起来依然自然。

在深度学习领域，TPS还被广泛应用于数据增强。训练一个识别模型时，往往需要对训练图像做各种变形以增加多样性，比如对OCR文字做细微弯曲、对车牌做透视之外的扭曲。用TPS对图像坐标做随机平滑形变，比单纯用仿射变换能产生更丰富的训练样本。这里有一个实操建议：形变的幅度要控制好，lambda可以设成0（严格插值），但控制点的随机位移幅度不能太大，一般控制在图像尺寸的2%到5%之间，否则生成的数据会失真，模型训练反而受影响。

5.3 别和性能测试里的TPS搞混

查资料时发现一个有趣现象：性能测试圈也有一个TPS，全称Transactions Per Second，每秒事务数。这个TPS衡量的是系统每秒能处理的完整事务数目，是压测报告里的核心指标。我见过有人搜“TPS虚高”搜到性能调优文章，又搜到一堆“怎么给某个交易限制TPS”的限流方案，结果跟薄板样条变换完全不搭边。

这两个TPS就像同名同姓的两个陌生人，一个住在图像和几何的数学世界，一个住在服务端性能测评的世界。如果你是在做图像配准、金属板形变、空间插值这类事情，你需要的Thin Plate Spline；如果你看的是吞吐量、压测、QPS、限流，那你需要的其实是Transactions Per Second。就连关键词“tps虚高”，说的也是压测时事务数虚高（比如错误请求被计入了成功率、缓存命中导致重复计算等），和薄板样条没有半毛钱关系。写这篇文章时我特意提这个，是因为我当年就因为这个同名缩写，在收藏夹里堆了一堆完全无关的链接。

6. 一个快速上手的小实验：用Python跑通TPS

6.1 伪代码：控制点求解与形变场计算

理论说了一堆，不如跑个最小例子。核心代码的逻辑如下（用的是最基础的方法，方便理解原理）：

import numpy as np def u(r): return r**2 * np.log(r + 1e-10) # 防止log(0) def build_tps_matrix(src_pts): n = len(src_pts) K = np.zeros((n, n)) for i in range(n): for j in range(n): r = np.linalg.norm(src_pts[i] - src_pts[j]) K[i, j] = u(r) P = np.hstack([np.ones((n, 1)), src_pts]) # n x 3 M = np.zeros((n + 3, n + 3)) M[:n, :n] = K M[:n, n:] = P M[n:, :n] = P.T # M[n:, n:] 保持零矩阵 return M def solve_tps(src_pts, dst_pts): n = len(src_pts) M = build_tps_matrix(src_pts) # dst_pts 可以是标量位移值，或者分坐标处理 b = np.zeros(n + 3) b[:n] = dst_pts # 注意：这里正则化lambda可以通过给K加上lambda*I来实现 w_a = np.linalg.solve(M, b) return w_a[:n], w_a[n:] def tps_map(src_pts, w, a, xy): # xy: (m, 2)的坐标点 # 结果 = 仿射项 + 径向基项 m = len(xy) res = a[0] + xy @ a[1:3] for i in range(len(src_pts)): r = np.linalg.norm(xy - src_pts[i], axis=1) res += w[i] * u(r) return res

注意代码里的np.log(r + 1e-10)是必要的数值保护，因为(r=0)时(r^2 \log r^2)在数学上的极限是0，但直接计算会得到-0或nan。加上小量后结果更稳定。上面示例里没有做坐标归一化，实际使用时一定要先对src_pts做中心化和尺度缩放。

如果要求解带正则化的版本，只需要在M[:n, :n]上加一个(\lambda I)即可，这样曲面就不必严格穿过每个控制点。更精细的做法是求解时用numpy.linalg.lstsq而不是solve，对矩阵奇异更鲁棒，但速度慢一些。

6.2 对整张图像做映射时的重采样细节

拿到待映射的任意像素坐标点之后，最有用的操作通常是逆向映射。怎么理解逆向映射？假设我们要把源图像变成目标图像，与其从源图的像素出发，看它跑到目标图的哪个位置，不如反过来：先定义目标图像的网格坐标（也就是最终输出图像每个像素的坐标），然后求这个网格坐标经过TPS逆变换映射回源图坐标的对应位置，再到源图上采样像素值。这样做能避免输出图像里出现空洞或重叠，因为目标图像的每个像素都能在源图找到采样点。

但这里立刻引出两个问题。第一，TPS正向形变是源到目标的映射，怎么求逆？如果你只有正向的TPS系数，求逆本身又是一个非线性优化问题。实操中常用的技巧是：交换控制点对重新求解。也就是说，解TPS时控制点输入从“源点→目标点”变成“目标点→源点”，这样直接就得到逆向的TPS模型，无需额外求逆。在控制点映射关系是一一对应的前提下，这个逆向TPS通常足够好。更复杂的情况可以用迭代最近点类算法逼近逆映射，但一般配准场景中交换控制点对就是首选方案。

第二，采样时要用插值。源图的像素坐标通常不是整数，最常用的是双线性插值，简单、快速、效果好；如果对质量有极致要求，可以用双三次插值，但计算量更大。这一步本质上是把离散图像离散网格之间的映射转化为连续坐标上的采值，插值方法会直接影响配准后图像的锐利度，对精度要求高的医学影像，通常会用高阶插值或者基于样条的插值方法。

实际写代码时还有一个坑：图像坐标系的原点通常位于左上角，而很多TPS实现假设原点在左下角（数学坐标系），如果不做Y轴翻转，形变结果会上下颠倒，而且很难一眼看出来。建议在代码中明确约定坐标系的定义，并在测试阶段用简单的平移控制点验证一下方向。

6.3 选型建议：什么时候用TPS，什么时候选别的

TPS不是万能的，也不该是唯一选择。我的经验是，判断标准很简单：形变是否局部、控制点是否有足够可信的对应关系。如果整幅图就是全局旋转加缩放，仿射变换就够了，用TPS纯属杀鸡用牛刀，还可能因为控制点噪声引入不自然的局部扭曲；如果形变极度局部且需要刚性保持（比如骨骼动画的局部转动），TPS又不如局部薄板样条（LPS）或者RBF插值的变种。

还有一种常见需求是处理“大位移形变”，比如两张图像拍摄视角差异极大，控制点对之间的位移可以达到图像尺寸的一半以上。这时候普通TPS会在中间区域产生剧烈的过冲或褶皱，因为笛卡尔坐标系下的位移被直接用径向基函数插值，大位移不满足小变形假设。一个改进方案是使用组合变换（先做全局仿射对齐，再做TPS残差），或者用多尺度金字塔策略，先粗后细地逐步配准。很多成熟配准框架（比如SimpleITK里的B样条配准）不直接用TPS，而是用B样条基函数做形变场建模，就是因为B样条对局部控制更灵活、数值性质也更好。所以，如果你发现TPS在某些大形变场景中表现不佳，不要硬调参，可以考虑换B样条形变模型。

还有一类场景，TPS容易被人忽略的优势是可解释性和对控制点的精确拟合。医学影像标注里的解剖点往往有明确的生物学含义，用TPS生成形变场后，你可以非常直观地看到每个关键点在形变前后的位置变化，这对医生审核结果很重要。而深度学习方法虽然可以取得更准确的配准结果，但可解释性、可复现性和受限于训练数据的泛化能力，在医疗等监管严格领域反而是劣势。所以TPS直到今天仍是很多配准任务中的稳健首选，不是因为它花哨，而是因为它可靠、可解释、可复现。

7. 一点个人心得作为收尾

最后分享一个我自己的经验。第一次实践TPS时，我对着一堆控制点直接求解，结果输出图上出现巨大的震荡波纹，我一度以为算法实现有误，排查了好久才发现是坐标没有归一化，控制点里有两点的坐标几乎是相邻像素，导致矩阵接近奇异。从那以后，我每次写TPS代码的第一件事永远是“坐标归一化+检查控制点距离”，第二件事是“在边界补点”，第三件事才是调(\lambda)。顺序颠倒过来，往往会浪费大量时间在错误方向上。

另外一个小技巧是：跑完TPS后不要只看控制点的误差，把形变场画出来——我习惯画网格线或者箭矢图，立刻就能判断形变是否合理。很多“看起来没问题其实在边缘已经失控”的情况，只有可视化才能发现。你现在如果有类似的项目在手，不妨从一小块区域开始，先跑通一个最小例子，再逐步扩展到全图，这个过程能帮你省下无数排查的时间。

如果你后续想深入，可以往这些方向继续挖：TPS与迭代最近点结合做自动配准、局部薄板样条（LPS）和多层级TPS处理大形变、以及把TPS嵌入深度学习框架作为可微空间变换层使用。薄板样条本身不是什么高深玄学，理清了它的物理图像、数学形式和计算细节，很多形变相关的场景都能顺理成章地解决。希望这篇文章能帮你在自己的项目里少踩几个坑，顺利把形变效果做出来。