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