1 | function [x1,x2]=witps(xp1,xp2,Y,Yp) |
---|
2 | % WITPS Inverse thin-plate spline warping |
---|
3 | % [X1,X2]=wtps(XP1,XP2,Y,Yp) computes the inverse thin-plate spline |
---|
4 | % (TPS) warp of the points XP1,XP2. |
---|
5 | % |
---|
6 | % REMARK. The inverse of a thin-plate spline in general is NOT a |
---|
7 | % thin-plate spline and some splines do not have an inverse. This |
---|
8 | % function uses Gauss-Newton to compute a set of points (X1,X2) such |
---|
9 | % that [XP1,XP2]=WTPS(X1,X2,Y,Yp). |
---|
10 | % |
---|
11 | % See also WTPS. |
---|
12 | |
---|
13 | % Initial guess by inverting the control points |
---|
14 | [x1,x2] = wtps(tps(xp1,xp2,Yp),Y) ; |
---|
15 | |
---|
16 | X = [x1(:)';x2(:)'] ; |
---|
17 | Xp = [xp1(:)',;xp2(:)'] ; |
---|
18 | |
---|
19 | % Gauss-Newton |
---|
20 | K = size(Y,2) ; |
---|
21 | N = size(X,2) ; |
---|
22 | U = tpsu(Y,Y) ; |
---|
23 | L = [[ones(1,K); Y], zeros(3) ; U, ones(K,1), Y'] ; |
---|
24 | invL = inv(L) ; |
---|
25 | A = [Yp, zeros(2,3)] * invL ; |
---|
26 | |
---|
27 | for t=1:5 |
---|
28 | [U,dU] = tpsu(Y,X); |
---|
29 | W = A * [repmat([0 0;1 0;0 1],1,N); reshape(dU, K, 2*N)] ; |
---|
30 | err = Xp - A * [ ones(1,N) ; X(1,:) ; X(2,:) ; U ] ; |
---|
31 | |
---|
32 | W = reshape(W,4,N) ; |
---|
33 | dets = W(1,:).*W(4,:) - W(3,:).*W(2,:) ; |
---|
34 | dX = [ ( W(4,:).*err(1,:) - W(3,:).*err(2,:) ) ./ dets ; ... |
---|
35 | (- W(2,:).*err(1,:) + W(1,:).*err(2,:) ) ./ dets ] ; |
---|
36 | X = X + dX ; |
---|
37 | end |
---|
38 | |
---|
39 | [M,N] = size(xp1) ; |
---|
40 | x1 = reshape(X(1,:),M,N) ; |
---|
41 | x2 = reshape(X(2,:),M,N) ; |
---|
42 | |
---|