[37] | 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 | |
---|