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

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

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 | |
---|

20 | function EM = tps_eval(dsites,ctrs) |
---|

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

22 | EM = 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) |
---|

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

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

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

33 | end |
---|

34 | |
---|

35 | % calculate tps |
---|

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

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

38 | |
---|

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

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