source: trunk/src/tps_coeff.m @ 356

Last change on this file since 356 was 356, checked in by sommeria, 12 years ago

civ updated with new functions for opening files, consistently with uvmat
Bugs to be expected (use previous version then)

File size: 1.7 KB
Line 
1%'tps_coeff': calculate the thin plate spline (tps) coefficients
2% (ref fasshauer@iit.edu MATH 590 ? Chapter 19 32)
3% this interpolation/smoothing minimises a linear combination of the squared curvature
4%  and squared difference form the initial data.
5% This function calculates the weight coefficients U_tps of the N sites where
6% data are known. Interpolated data are then obtained as the matrix product
7% EM*U_tps where the matrix EM is obtained by the function tps_eval.
8% The spatial derivatives are obtained as EMDX*U_tps and EMDY*U_tps, where
9% EMDX and EMDY are obtained from the function tps_eval_dxy.
10%------------------------------------------------------------------------
11% [U_smooth,U_tps]=tps_coeff(ctrs,U,rho)
12%------------------------------------------------------------------------
13% OUPUT:
14% U_smooth: values of the quantity U at the N centres after smoothing
15% U_tps: tps weights of the centres
16
17%INPUT:
18% ctrs: Nxs matrix  representing the postions of the M centers, sources of the tps (s=space dimension)
19% U: Nx1 column vector representing the initial values of the considered scalar at the centres ctrs
20% rho: smoothing parameter: the result is smoother for larger rho.
21
22
23function [U_smooth,U_tps]=tps_coeff(ctrs,U,rho)
24%------------------------------------------------------------------------
25%rho smoothing parameter
26% X=reshape(X,[],1);
27% Y=reshape(Y,[],1);
28% N=numel(X);
29% rhs = reshape(U,[],1);
30U = [U; zeros(3,1)];
31% ctrs = [X Y];% coordinates of measurement sites, radial base functions are located at the measurement sites
32EM = tps_eval(ctrs,ctrs);
33RhoMat=rho*eye(N,N);%  rho=1/(2*omega) , omega given by fasshauer;
34RhoMat=[RhoMat zeros(N,3)];
35PM=[ones(N,1) ctrs];
36IM=[EM+RhoMat; [PM' zeros(3,3)]];
37U_tps=(IM\U);
38U_smooth=EM *U_tps;
Note: See TracBrowser for help on using the repository browser.