Changeset 356 for trunk/src/tps_coeff.m

Timestamp:

Jan 3, 2012, 12:58:52 AM (13 years ago)

Author:

sommeria

Message:

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

File:

• trunk/src/tps_coeff.m (modified) (1 diff)

trunk/src/tps_coeff.m

@@ -1 @@
-%'tps_uvmat': calculate thin plate shell coefficients
+%'tps_coeff': calculate the thin plate spline (tps) coefficients
+% (ref fasshauer@iit.edu MATH 590 ? Chapter 19 32)
+% this interpolation/smoothing minimises a linear combination of the squared curvature
+%  and squared difference form the initial data. 
+% This function calculates the weight coefficients U_tps of the N sites where
+% data are known. Interpolated data are then obtained as the matrix product
+% EM*U_tps where the matrix EM is obtained by the function tps_eval.
+% The spatial derivatives are obtained as EMDX*U_tps and EMDY*U_tps, where
+% EMDX and EMDY are obtained from the function tps_eval_dxy.
 %------------------------------------------------------------------------
-%fasshauer@iit.edu MATH 590 ? Chapter 19 32
-% X,Y initial coordiantes
-% XI vector, YI column vector for the grid of interpolation points
-function [U_smooth,U_tps]=tps_coeff(X,Y,U,rho)
+% [U_smooth,U_tps]=tps_coeff(ctrs,U,rho)
+%------------------------------------------------------------------------
+% OUPUT:
+% U_smooth: values of the quantity U at the N centres after smoothing
+% U_tps: tps weights of the centres
+
+%INPUT:
+% ctrs: Nxs matrix  representing the postions of the M centers, sources of the tps (s=space dimension)
+% U: Nx1 column vector representing the initial values of the considered scalar at the centres ctrs
+% rho: smoothing parameter: the result is smoother for larger rho.
+
+
+function [U_smooth,U_tps]=tps_coeff(ctrs,U,rho)
 %------------------------------------------------------------------------
 %rho smoothing parameter
-% ep = 1; 
-X=reshape(X,[],1);
-Y=reshape(Y,[],1);
-N=numel(X);
-rhs = reshape(U,[],1);
-rhs = [rhs; zeros(3,1)];
-ctrs = [X Y];% coordinates of measurement sites, radial base functions are located at the measurement sites
-%ctrs = dsites;%radial base functions are located at the measurement sites
+% X=reshape(X,[],1);
+% Y=reshape(Y,[],1);
+% N=numel(X);
+% rhs = reshape(U,[],1);
+U = [U; zeros(3,1)];
+% ctrs = [X Y];% coordinates of measurement sites, radial base functions are located at the measurement sites
 EM = tps_eval(ctrs,ctrs);
-% DM_data = DistanceMatrix(ctrs,ctrs);%2D matrix of distances between spline centres (=initial points) ctrs
-% IM_sites = tps(1,DM_data);%values of thin plate at site points
-% PM=[ones(N,1) ctrs];
-% EM = [IM_sites PM];
-%IM = IM_sites + rho*eye(size(IM_sites));%  rho=1/(2*omega) , omega given by fasshauer;
-
-%IM=[IM PM; [PM' zeros(3,3)]];
 RhoMat=rho*eye(N,N);%  rho=1/(2*omega) , omega given by fasshauer;
 RhoMat=[RhoMat zeros(N,3)];
 PM=[ones(N,1) ctrs];
-
 IM=[EM+RhoMat; [PM' zeros(3,3)]];
-%fprintf('Condition number estimate: %e\n',condest(IM))
-%DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
-%EM = tps(ep,DM_eval);%values of thin plate 
-%PM = [ones(size(epoints,1),1) epoints]; 
-%EM = [EM PM];
-U_tps=(IM\rhs);
-% PM = [ones(size(dsites,1),1) dsites]; 
-
+U_tps=(IM\U);
 U_smooth=EM *U_tps;