关于两个对应的非刚体变换的2D点集的匹配，使用Thin Plate Spline,求助

dianyancao

http://en.wikipedia.org/wiki/Thin_plate_spline

假设x为二维向量，那么这里是需要对向量x的x轴分量与y轴分量分别计算Cix和Ciy,而生成一个(2+1)*K的矩阵C吗？
      [C1x  ...  Cix  ...  CKx]
C=[C1y  ...  Ciy  ...  Cky]    ?
      [1  ...  1  ...  1]


这里对||·||中的(2+1)*K矩阵，是如何计算范数？

呵呵，本人只有初中文化，很多东西还没有学会，打扰之处还望大家见谅。
怀疑上述链接的矩阵行和列和通常相比是转置的，可以从以下链接看出：
http://en.wikipedia.org/wiki/Affine_transformation

wayne

1# dianyancao
你给的链接已经做了很明确的解释，如果你英文足够好的话。
||x-w|| 表示求欧氏距离，即坐标分量相减的平方和，再开方。
phi(x) 表示核函数，如高斯核

dianyancao

感谢wayne的解答，请问上述矩阵的范数是矩阵元 p-范数(p=2)吗？

dianyancao

不知道我以下的理解是否有误，



ftps(zn)等价于f(xn,yn) = (g(xn), g(yn)),其中

g(xn) = m11·xn + m12·yn + m13 + Ki:Cx

g(yn) = m21·xn + m22·yn + m23 + Ki:Cy



zn= [ xn  yn  1]T

d是一个全局仿射矩阵，提供的参数为:

[m11  m12  m13]

d = [m21  m22  m23]

[0    0    1]

Cx和Cy分别是与向量z相关的x坐标测度和y坐标测度。

Cx = [C1x    C2x    C3x   …...   Cnx]T

Cy = [C1y    C2y    C3y   …...   Cny]T



Ki:为径向基核函数的提供矩阵的第i行矢量

Phi(rij)= rij2·log rij

rij = |zi - zj|

[  0      Phi(r12)  …   Phi(r1n)]

K = [Phi(r21)     0    …    Phi(r2n)]

[  …      …    0 ..     …  ]

[Phi(rn1)   Phi(rn2)  …     0   ]



按上述的链接中指出



等价于



由于不知道该矩阵采用的范数，所以用第二钟方法验证，

但是在计算二阶偏导数(der f2)/(der x2)的时候遇到困难，由于f的值是一个向量，

而x是该向量的x轴坐标，得到的结果：

lim

delta(x) -> 0

(m11·(xn+ delta(x)) + m12·yn + m13 + Ki:Cx ，m21·(xn+ delta(x)) + m22·yn + m23 + Ki:Cy) - (m11·xn + m12·yn + m13 + Ki:Cx ，m21·xn + m22·yn + m23 + Ki:Cy)

(der f)/(der x) =









              

                                                      delta(x)

            

lim

delta(x) -> 0

=

(m11·delta(x), m21·delta(x))

                            delta(x)



  = (m11, m21)



同理得到：

(der f)/(der y) = (m12,m22)

不知道以上结果是否正确，这些一阶偏导数均为常矢量，和x，y无关，

接下来如何求其二阶偏导数 (der ((der f)/(der x))/der x)

dianyancao

还是发图片吧...

dianyancao

上面计算偏导数时的Ki:Cx是不正确的，借助与Matlab，得到：
Etps=symsum(8*symsum((CX(i)*(XX(i)-XX(n))*(YY(i)-YY(n)))/(XX(i)^2+YY(i)^2+XX(n)^2+YY(n)^2-2*XX(i)*XX(n)-2*YY(i)*YY(n)),i,1,M)^2+8*symsum((CY(i)*(XX(i)-XX(n))*(YY(i)-YY(n)))/(XX(i)^2+YY(i)^2+XX(n)^2+YY(n)^2-2*XX(i)*XX(n)-2*YY(i)*YY(n)),i,1,M)^2+symsum(CX(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CX(i)+(2*CX(i)*(XX(i)-XX(n))*(2*XX(i)-2*XX(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CX(i)*(XX(i)-XX(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2+symsum(CY(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CY(i)+(2*CY(i)*(XX(i)-XX(n))*(2*XX(i)-2*XX(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CY(i)*(XX(i)-XX(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2+symsum(CX(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CX(i)+(2*CX(i)*(YY(i)-YY(n))*(2*YY(i)-2*YY(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CX(i)*(YY(i)-YY(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2+symsum(CY(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CY(i)+(2*CY(i)*(YY(i)-YY(n))*(2*YY(i)-2*YY(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CY(i)*(YY(i)-YY(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2,n,1,i-1)+symsum(8*symsum((CX(i)*(XX(i)-XX(n))*(YY(i)-YY(n)))/(XX(i)^2+YY(i)^2+XX(n)^2+YY(n)^2-2*XX(i)*XX(n)-2*YY(i)*YY(n)),i,1,M)^2+8*symsum((CY(i)*(XX(i)-XX(n))*(YY(i)-YY(n)))/(XX(i)^2+YY(i)^2+XX(n)^2+YY(n)^2-2*XX(i)*XX(n)-2*YY(i)*YY(n)),i,1,M)^2+symsum(CX(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CX(i)+(2*CX(i)*(XX(i)-XX(n))*(2*XX(i)-2*XX(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CX(i)*(XX(i)-XX(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2+symsum(CY(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CY(i)+(2*CY(i)*(XX(i)-XX(n))*(2*XX(i)-2*XX(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CY(i)*(XX(i)-XX(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2+symsum(CX(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CX(i)+(2*CX(i)*(YY(i)-YY(n))*(2*YY(i)-2*YY(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CX(i)*(YY(i)-YY(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2+symsum(CY(i)+log((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)*CY(i)+(2*CY(i)*(YY(i)-YY(n))*(2*YY(i)-2*YY(n)))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)-(3*CY(i)*(YY(i)-YY(n))^2*((XX(i)-XX(n))^2/2+(YY(i)-YY(n))^2/2))/((XX(i)-XX(n))^2+(YY(i)-YY(n))^2)^2,i,1,M)^2,n,i+1,M)

dianyancao

以下是对每个对应点求和，不知道能不能代替反常积分？
i != n

wayne

楼主研究这个好像有什么图谋，能说说吗
不介意我这样问吧。呵呵

dianyancao

额..是有的，可以用来对需要变形的形状进行插值，提高形状匹配的准确率，比如识别验证码之类的(干坏事，呵呵)，但是也不局限于此，还可以用于医学图像配准，地质建模等

dianyancao

杯具啊，单求一个全局最优的仿射矩阵就难倒我了
以下是Mathematica的代码，貌似求那个最高次数为2次的多项式很复杂

In[2]:= m = {{m11, m21, 0}, {m12, m22, 0}, {m13, m23, 1}}

Out[2]= {{m11, m21, 0}, {m12, m22, 0}, {m13, m23, 1}}

In[3]:= a = Table[{ax[i], ay[i], 1}, {i, 1, 4}]

Out[3]= {{ax[1], ay[1], 1}, {ax[2], ay[2], 1}, {ax[3], ay[3], 1}, {ax[4], ay[4], 1}}

In[4]:= b = Table[{bx[i], by[i], 1}, {i, 1, 4}]

Out[4]= {{bx[1], by[1], 1}, {bx[2], by[2], 1}, {bx[3], by[3], 1}, {bx[4], by[4], 1}}


In[5]:= s = 0

Out[5]= 0

In[6]:= Norm[a.m - b, 2]^2

Out[6]= Out[6]



In[7]:= For[i = 1, i < 5, i++,
s += (m13 + m11 ax[i] + m12 ay[i] - bx[i])^2 + (m23 + m21 ax[i] + m22 ay[i] -
       by[i])^2]

In[8]:= s

Out[8]= (m13 + m11 ax[1] + m12 ay[1] - bx[1])^2 + (m13 + m11 ax[2] + m12 ay[2] -
   bx[2])^2 + (m13 + m11 ax[3] + m12 ay[3] - bx[3])^2 + (m13 + m11 ax[4] +
   m12 ay[4] - bx[4])^2 + (m23 + m21 ax[1] + m22 ay[1] - by[1])^2 + (m23 +
   m21 ax[2] + m22 ay[2] - by[2])^2 + (m23 + m21 ax[3] + m22 ay[3] -
   by[3])^2 + (m23 + m21 ax[4] + m22 ay[4] - by[4])^2





In[9]:= Minimize[s, {m11, m21, m12, m22, m13, m23}]

During evaluation of In[15]:=

No more memory available.
Mathematica kernel has shut down.
Try quitting other applications and then retry.