# source:trunk/src/tps_eval.m@660

Last change on this file since 660 was 581, checked in by sommeria, 8 years ago

clean the transform field functions

File size: 1.6 KB
Line
1%'tps_eval': calculate the thin plate spline (tps) interpolation at a set of points
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.
10%  use : U_interp=EM*U_tps
11%
12%INPUT:
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,
15%
16% related functions:
17% tps_coeff, tps_eval_dxy
18% tps_coeff_field, set_subdomains, filter_tps, calc_field
19
20function EM = tps_eval(dsites,ctrs)
21[M,s] = size(dsites); [N,s] = size(ctrs);
22EM = zeros(M,N);
23
24% calculate distance matrix: accumulate sum of squares of coordinate differences
25% The ndgrid command produces two MxN matrices:
26%   Dsite, consisting of N identical columns (each containing
27%       the d-th coordinate of the M data sites)
28%   Ctrs, consisting of M identical rows (each containing
29%       the d-th coordinate of the N centers)
30for d=1:s
31 [Dsites,Ctrs] = ndgrid(dsites(:,d),ctrs(:,d));
32 EM = EM + (Dsites-Ctrs).^2;%EM=square of distance matrices
33end
34
35% calculate tps
36np=find(EM~=0);
37EM(np) = EM(np).*log(EM(np))/2;%= tps formula r^2 log(r) (EM=r^2)
38