1 | %'tps_uvmat': calculate thin plate shell coefficients |
---|

2 | %------------------------------------------------------------------------ |
---|

3 | %fasshauer@iit.edu MATH 590 ? Chapter 19 32 |
---|

4 | % X,Y initial coordiantes |
---|

5 | % XI vector, YI column vector for the grid of interpolation points |
---|

6 | function [U_smooth,U_tps]=tps_coeff(X,Y,U,rho) |
---|

7 | %------------------------------------------------------------------------ |
---|

8 | %rho smoothing parameter |
---|

9 | % ep = 1; |
---|

10 | X=reshape(X,[],1); |
---|

11 | Y=reshape(Y,[],1); |
---|

12 | N=numel(X); |
---|

13 | rhs = reshape(U,[],1); |
---|

14 | rhs = [rhs; zeros(3,1)]; |
---|

15 | ctrs = [X Y];% coordinates of measurement sites, radial base functions are located at the measurement sites |
---|

16 | %ctrs = dsites;%radial base functions are located at the measurement sites |
---|

17 | EM = tps_eval(ctrs,ctrs); |
---|

18 | % DM_data = DistanceMatrix(ctrs,ctrs);%2D matrix of distances between spline centres (=initial points) ctrs |
---|

19 | % IM_sites = tps(1,DM_data);%values of thin plate at site points |
---|

20 | % PM=[ones(N,1) ctrs]; |
---|

21 | % EM = [IM_sites PM]; |
---|

22 | %IM = IM_sites + rho*eye(size(IM_sites));% rho=1/(2*omega) , omega given by fasshauer; |
---|

23 | |
---|

24 | %IM=[IM PM; [PM' zeros(3,3)]]; |
---|

25 | RhoMat=rho*eye(N,N);% rho=1/(2*omega) , omega given by fasshauer; |
---|

26 | RhoMat=[RhoMat zeros(N,3)]; |
---|

27 | PM=[ones(N,1) ctrs]; |
---|

28 | |
---|

29 | IM=[EM+RhoMat; [PM' zeros(3,3)]]; |
---|

30 | %fprintf('Condition number estimate: %e\n',condest(IM)) |
---|

31 | %DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs |
---|

32 | %EM = tps(ep,DM_eval);%values of thin plate |
---|

33 | %PM = [ones(size(epoints,1),1) epoints]; |
---|

34 | %EM = [EM PM]; |
---|

35 | U_tps=(IM\rhs); |
---|

36 | % PM = [ones(size(dsites,1),1) dsites]; |
---|

37 | |
---|

38 | U_smooth=EM *U_tps; |
---|