Date: 2015-06-30 21:53

1. 基本数学描述

薄板样条(Thin Plate Spline)映射根据两幅相关图像中的对应控制点集来决定一个变形函数。它寻找通过所有给定点的饶度最小的光滑曲面。“薄板”这个名字的由来，就表示薄板样条是用来近似的仿真一块金属薄片在通过相同的控制点时的行为特征。

上图是使用薄板样条进行变换的例子，薄板样条插值可表示为多变量插值问题，对于d维空间中提取出的n个特征点对xi和yi,问题描述为在一个合适的Hilbert空间H上寻找连续变换f:Rd→Rd,并满足如下两个条件:

(1)使得既定泛函E(f):H→R最小化;

(2)满足插值条件(通过指定的几个点)。

(这里有一系列的数学推导，一维如何计算，二维如何计算，三维如何计算，四维如何计算，暂时代过 … …)

2. 二维图像配准中的应用

以二维图像配准为例,在2维2阶导数情况下,则f(x,y)在满足插值方程的同时使得能量函数达到最小，即有：

f=argminfEfE(f)=∫∫⎛⎝(∂2f∂x2)2+2(∂2f∂x∂y)2+(∂2f∂y2)⎞⎠dxdy

此即为Bookstein于1989年首次将薄板样条应用于二维图像配准中时所采用的目标函数.使得能量函数取得最优解的

f(x,y)

的表达式为：

f(x,y)=a1+axx+ayy+∑i=1nwiU(∥(xi,yi)−(x,y)∥)

其中

U(r)=r2logr

为径向基函数。

wi,i=1,2,⋯,n

以及

a1,ax,ay

为未知系数。

为了使

f(x,y)

的二阶导数平方可积，必须有：

∑i=1nwi=∑i=1nwixi=∑i=1nwiyi=0

故可得薄板样条函数的系数满足一个线性方程组：

[KPTPO][wa]=[vo]

其中Ki,j=U(∥(xi,yi)−(xj,yj)∥)，P的第i行为(1,xi,yi)，O为3×3的零矩阵，o为3×1零向量，w和v为由wi和vi得到的列向量，而列向量a中元素为a1,ax,ay。令

L=[KPTPO]

对上述线性方程组求解，可得：

w=L−1M

薄板样条的映射函数可以确定源图像(图像A)到目标图像(图像B)的映射变换关键系数(w1,w2,⋯,wn;a1,ax,ay)，然后将图像A中任意一点的坐标代入公式，可得到图像B中对应点的坐标。但是，薄板映射并不是一对一的，通常若干个点会映射到同一个点上。这就需要对整幅图像进行遍历，找到所有漏掉的点，通过插值得到这些点的像素值。插值的过程这里不再过多介绍。

3. 结合Matlab代码说明

主函数如下，基中，Xp和Yp分别指源图像中的特征点，Xs和Ys指参考图像中的特征点：

%% Algebra of Thin-plate splines

% Compute thin-plate spline mapping [W|a1 ax ay] using landmarks

[wL]=computeWl(Xp, Yp, NPs);

% Y = ( V| 0 0 0)' where V = [G] where G is landmark homologous (nx2) ;

% Y is col vector of length (n+3)

wY = [Xs(:) Ys(:); zeros(3,2)];

wW = inv(wL)*wY; % (W|a1 ax ay)' = inv(L)*Y

% Thin-plate spline mapping (Map all points in the plane)

% f(x,y) = a1 + ax * x + ay * y + SUM(wi * U(|Pi-(x,y)|)) for i = 1 to n

[Xw, Yw]=tpsMap(wW, imgH, imgW, Xp, Yp, NPs);

从中可以看出，首先是计算L矩阵，代码如下：

%% [L] = [[K P];[P' 0]]

% np - number of landmark points

% (xp, yp) - coordinate of landmark points

function [wL]=computeWl(xp, yp, np)

rXp = repmat(xp(:),1,np); % 1xNp to NpxNp

rYp = repmat(yp(:),1,np); % 1xNp to NpxNp

wR = sqrt((rXp-rXp').^2 + (rYp-rYp').^2); % compute r(i,j)

wK = radialBasis(wR); % compute [K] with elements U(r)=r^2 * log (r^2)

wP = [ones(np,1) xp(:) yp(:)]; % [P] = [1 xp' yp'] where (xp',yp') are n landmark points (nx2)

wL = [wK wP;wP' zeros(3,3)]; % [L] = [[K P];[P' 0]]

return

得以L矩阵后，下面一行代码用来对w求解：

wW = inv(wL)*wY;

最后利用下面一行代码实现源图像到目标图像的变换：

[Xw, Yw]=tpsMap(wW, imgH, imgW, Xp, Yp, NPs);

其中，tpsMap的代码如下：

%% Mapping: f(x,y) = a1 + ax * x + ay * y + SUM(wi * U(|Pi-(x,y)|)) for i = 1 to n

% np - number of landmark points

% (xp, yp) - coordinate of landmark points

function [Xw, Yw]=tpsMap(wW, imgH, imgW, xp, yp, np)

[X Y] = meshgrid(1:imgH,1:imgW); % HxW

X=X(:)'; % convert to 1D array by reading columnwise (NWs=H*W) Y=Y(:)'; % convert to 1D array (NWs)

NWs = length(X); % total number of points in the plane

% all points in plane

rX = repmat(X,np,1); % Np x NWs

rY = repmat(Y,np,1); % Np x NWs

% landmark points

rxp = repmat(xp(:),1,NWs); % 1xNp to Np x NWs

ryp = repmat(yp(:),1,NWs); % 1xNp to Np x NWs

% Mapping Algebra

wR = sqrt((rxp-rX).^2 + (ryp-rY).^2); % distance measure r(i,j)=|Pi-(x,y)|

wK = radialBasis(wR); % compute [K] with elements U(r)=r^2 * log (r)

wP = [ones(NWs,1) X(:) Y(:)]'; % [P] = [1 x' y']'

wL = [wK;wP]'; % [L] = [w1 w2 ... wn 1 x y]

Xw = wL*wW(:,1); % [Pw] = [L]*[W]

Yw = wL*wW(:,2); % [Pw] = [L]*[W]

return

薄板样条插值用于图像配准的原理是非常容易理解的。此外，代码中还有一些是用于插值填图像孔洞 ，这里就不介绍了

注意：这个代码中的X，Y和我们传统图像处理中的X，Y是反的，也就是说，X是行，Y是列，作者为什么这么写，不清楚了

其它细节可以参考如下论文：

F.L. Bookstein, “Principal Warps: Thin-Plate Splines and the Decomposition of Deformations,” IEEE TPAMI, 1989.

张春生, 基于CT图像的肺部轮廓非刚性配准方法[D], 哈尔滨工业大学 2012.

李菁菁, 基于控制点平滑的人脸变形算法及其在人脸动画中的应用[D], 湘潭大学 2008.