1 | function [U,dU,delta]=tpsu(X,Y) |
---|
2 | % TPSU Compute the U matrix of a thin-plate spline transformation |
---|
3 | % U=TPSU(X,Y) returns the matrix |
---|
4 | % |
---|
5 | % [ U(|X(:,1) - Y(:,1)|) ... U(|X(:,1) - Y(:,N)|) ] |
---|
6 | % [ ] |
---|
7 | % [ U(|X(:,M) - Y(:,1)|) ... U(|X(:,M) - Y(:,N)|) ] |
---|
8 | % |
---|
9 | % where X is a 2xM matrix and Y a 2xN matrix of points and U(r) is |
---|
10 | % the opposite -r^2 log(r^2) of the radial basis function of the |
---|
11 | % thin plate spline specified by X and Y. |
---|
12 | % |
---|
13 | % [U,dU]=tpsu(x,y) returns the derivatives of the columns of U with |
---|
14 | % respect to the parameters Y. The derivatives are arranged in a |
---|
15 | % Mx2xN array, one layer per column of U. |
---|
16 | % |
---|
17 | % See TPS. |
---|
18 | |
---|
19 | if exist('tpsumx') |
---|
20 | U = tpsumx(X,Y) ; |
---|
21 | else |
---|
22 | M=size(X,2) ; |
---|
23 | N=size(Y,2) ; |
---|
24 | |
---|
25 | % Faster than repmat, but still fairly slow |
---|
26 | r2 = ... |
---|
27 | (X( ones(N,1), :)' - Y( ones(1,M), :)).^2 + ... |
---|
28 | (X( 1+ones(N,1), :)' - Y(1+ones(1,M), :)).^2 ; |
---|
29 | U = - rb(r2) ; |
---|
30 | end |
---|
31 | |
---|
32 | if nargout > 1 |
---|
33 | M=size(X,2) ; |
---|
34 | N=size(Y,2) ; |
---|
35 | |
---|
36 | dx = X( ones(N,1), :)' - Y( ones(1,M), :) ; |
---|
37 | dy = X(1+ones(N,1), :)' - Y(1+ones(1,M), :) ; |
---|
38 | r2 = (dx.^2 + dy.^2) ; |
---|
39 | r = sqrt(r2) ; |
---|
40 | coeff = drb(r)./(r+eps) ; |
---|
41 | dU = reshape( [coeff .* dx ; coeff .* dy], M, 2, N) ; |
---|
42 | end |
---|
43 | |
---|
44 | % The radial basis function |
---|
45 | function y = rb(r2) |
---|
46 | y = zeros(size(r2)) ; |
---|
47 | sel = find(r2 ~= 0) ; |
---|
48 | y(sel) = - r2(sel) .* log(r2(sel)) ; |
---|
49 | |
---|
50 | % The derivative of the radial basis function |
---|
51 | function y = drb(r) |
---|
52 | y = zeros(size(r)) ; |
---|
53 | sel = find(r ~= 0) ; |
---|
54 | y(sel) = - 4 * r(sel) .* log(r(sel)) - 2 * r(sel) ; |
---|
55 | |
---|
56 | |
---|