Changeset 356 for trunk/src/tps_coeff.m


Ignore:
Timestamp:
Jan 3, 2012, 12:58:52 AM (9 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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/tps_coeff.m

    r246 r356  
    1 %'tps_uvmat': calculate thin plate shell coefficients
     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.
    210%------------------------------------------------------------------------
    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)
     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)
    724%------------------------------------------------------------------------
    825%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
     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
    1732EM = 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)]];
    2533RhoMat=rho*eye(N,N);%  rho=1/(2*omega) , omega given by fasshauer;
    2634RhoMat=[RhoMat zeros(N,3)];
    2735PM=[ones(N,1) ctrs];
    28 
    2936IM=[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 
     37U_tps=(IM\U);
    3838U_smooth=EM *U_tps;
Note: See TracChangeset for help on using the changeset viewer.