source: trunk/src/tps_coeff.m @ 643

Last change on this file since 643 was 586, checked in by sommeria, 11 years ago

tps_coeff_field introduced and several bugs corrected

File size: 1.9 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% for big data sets, a splitting in subdomains is needed, see functions
11% set_subdomains and tps_coeff_field.
12%
13%------------------------------------------------------------------------
14% [U_smooth,U_tps]=tps_coeff(ctrs,U,Smoothing)
15%------------------------------------------------------------------------
16% OUPUT:
17% U_smooth: values of the quantity U at the N centres after smoothing
18% U_tps: tps weights of the centres and columns of the linear
19
20%INPUT:
21% ctrs: NxNbDim matrix  representing the positions of the N centers, sources of the tps (NbDim=space dimension)
22% U: Nx1 column vector representing the values of the considered scalar measured at the centres ctrs
23% Smoothing: smoothing parameter: the result is smoother for larger Smoothing.
24%
25%related functions:
26% tps_eval, tps_eval_dxy
27% tps_coeff_field, set_subdomains, filter_tps, calc_field
28
29function [U_smooth,U_tps]=tps_coeff(ctrs,U,Smoothing)
30%------------------------------------------------------------------------
31warning off
32N=size(ctrs,1);% nbre of source centres
33NbDim=size(ctrs,2);% space dimension (2 or 3)
34U = [U; zeros(NbDim+1,1)];
35EM = tps_eval(ctrs,ctrs);
36SmoothingMat=Smoothing*eye(N,N);%  Smoothing=1/(2*omega) , omega given by fasshauer;
37SmoothingMat=[SmoothingMat zeros(N,NbDim+1)];
38PM=[ones(N,1) ctrs];
39IM=[EM+SmoothingMat; [PM' zeros(NbDim+1,NbDim+1)]];
40U_tps=(IM\U);
41U_smooth=EM *U_tps;
Note: See TracBrowser for help on using the repository browser.