[356] | 1 | %'tps_eval': calculate the thin plate spline (tps) interpolation at a set of points |
---|

[434] | 2 | % see tps_coeff.m for more information and test_tps.m for an example |
---|

[356] | 3 | %------------------------------------------------------------------------ |
---|

| 4 | % function EM = tps_eval(dsites,ctrs) |
---|

| 5 | %------------------------------------------------------------------------ |
---|

| 6 | % OUPUT: |
---|

| 7 | % EM: Mx(N+s) matrix representing the contributions at the M sites |
---|

| 8 | % from unit sources located at each of the N centers, + (s+1) columns |
---|

| 9 | % representing the contribution of the linear gradient part. |
---|

[434] | 10 | % use : U_interp=EM*U_tps |
---|

[356] | 11 | % |
---|

| 12 | %INPUT: |
---|

[434] | 13 | %dsites: Mxs matrix representing the postions of the M 'observation' sites, with s the space dimension |
---|

| 14 | %ctrs: Nxs matrix representing the postions of the N centers, sources of the tps, |
---|

[356] | 15 | % |
---|

| 16 | % related functions: |
---|

| 17 | % tps_coeff, tps_eval_dxy |
---|

[246] | 18 | |
---|

[356] | 19 | function EM = tps_eval(dsites,ctrs) |
---|

| 20 | [M,s] = size(dsites); [N,s] = size(ctrs); |
---|

| 21 | EM = zeros(M,N); |
---|

| 22 | |
---|

| 23 | % calculate distance matrix: accumulate sum of squares of coordinate differences |
---|

| 24 | % The ndgrid command produces two MxN matrices: |
---|

| 25 | % Dsite, consisting of N identical columns (each containing |
---|

| 26 | % the d-th coordinate of the M data sites) |
---|

| 27 | % Ctrs, consisting of M identical rows (each containing |
---|

| 28 | % the d-th coordinate of the N centers) |
---|

| 29 | for d=1:s |
---|

| 30 | [Dsites,Ctrs] = ndgrid(dsites(:,d),ctrs(:,d)); |
---|

| 31 | EM = EM + (Dsites-Ctrs).^2;%EM=square of distance matrices |
---|

| 32 | end |
---|

| 33 | |
---|

| 34 | % calculate tps |
---|

| 35 | np=find(EM~=0); |
---|

| 36 | EM(np) = EM(np).*log(EM(np))/2;%= tps formula r^2 log(r) (EM=r^2) |
---|

| 37 | |
---|

| 38 | % add linear gradient part: |
---|

| 39 | EM = [EM ones(M,1) dsites]; |
---|