source: trunk/src/tps_coeff.m @ 504

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

corrections in the use of get_field
test_tps introduced, to test thin plate shell functions
tps-eval_dxy corrected (bug on the calculation of the derivatives fixed)

File size: 1.8 KB
RevLine 
[356]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.
[246]10%------------------------------------------------------------------------
[434]11% [U_smooth,U_tps]=tps_coeff(ctrs,U,Smoothing)
[246]12%------------------------------------------------------------------------
[356]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:
[434]18% ctrs: Nxs matrix  representing the postions of the N centers, sources of the tps (s=space dimension)
19% U: Nx1 column vector representing the values of the considered scalar measured at the centres ctrs
20% Smoothing: smoothing parameter: the result is smoother for larger Smoothing.
[356]21
22
[434]23function [U_smooth,U_tps]=tps_coeff(ctrs,U,Smoothing)
[356]24%------------------------------------------------------------------------
[434]25%Smoothing smoothing parameter
[356]26% X=reshape(X,[],1);
27% Y=reshape(Y,[],1);
[362]28N=size(ctrs,1);
[356]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
[246]32EM = tps_eval(ctrs,ctrs);
[434]33SmoothingMat=Smoothing*eye(N,N);%  Smoothing=1/(2*omega) , omega given by fasshauer;
34SmoothingMat=[SmoothingMat zeros(N,3)];
[246]35PM=[ones(N,1) ctrs];
[434]36IM=[EM+SmoothingMat; [PM' zeros(3,3)]];
[356]37U_tps=(IM\U);
[246]38U_smooth=EM *U_tps;
Note: See TracBrowser for help on using the repository browser.