[julia-users] Re: Thin Plate Splines

[Jan Kybic]

>
> I have a set of irregularly gridded data (x,y,z) and I am trying to create 
> an interpolating surface using Thin Plate Splines. I couldn't find any 
> existing Julia routines so I thought I'd just do it my self. Here is my 
> implementation 
> <http://nbviewer.ipython.org/urls/gist.githubusercontent.com/lstagner/89964e63557ddb831e72/raw/49049103742450d80ee28d8876fa1a8b8037f970/tps.ipynb>.
>  
> As you can see its wrong. I been staring at it for a while now and I am 
> beginning to think I must be hitting some sort of bug or quirk of the 
> language. It either that or I did something wrong. If I get this to work I 
> was thinking about incorporating it into one of the existing interpolation 
> packages.
>
> Can anyone figure out why this is not working?
>
>
 Hello. I have never seen the formulas you are using to find 'a' and 'b' 
and I do not have the time to check. Provided they are correct, I would 
suspect numerical problems - it is often not a good idea to calculate a 
matrix inverse explicitely and the TPS matrix is ill-conditioned for even a 
moderate number of points. I am attaching my Matlab implementation - you 
see that I am solving the linear system by SVD, which is more stable. It 
should be trivial to translate to Julia. Hope this helps.

Jan

function [c,cnd]=lndthscoef(x,z) ;
% Calculate coefficients for thin-plate spline interpolation
% Calculates coefficients (usable e.g. by lndthseval) to interpolate
% function having values z at points x
%
% Usage: c=lndthscoef(x,z) 
%
% $Id: lndthscoef.m,v 1.3 2000/01/27 18:11:53 jkybic Exp $
% Jan Kybic, 1999

n=length(z) ;
A=zeros(n+3) ;

for i=1:n,
  for j=1:n,
    r2=(x(i,1)-x(j,1))^2+(x(i,2)-x(j,2))^2 ;
    if r2==0,
      A(i,j)=0 ;
    else
      A(i,j)=r2*log(r2) ; % actually twice as much as r^2 ln(r) 
    end ;
  end ;
  A(i,n+1:n+3)=[ x(i,1:2) 1 ] ; 
end ;

A(n+1,1:n)=x(:,1)' ;
A(n+2,1:n)=x(:,2)' ;
A(n+3,1:n)=ones(1,n) ;

% Solve using SVD techniques
[u,s,v]=svd(A) ;
s=diag(s) ; 
cnd=s(1)/s(n+1) ;
q=(s(1)*1e-20+1e-50) ; % the q-value is sensitive !!!
for i=1:n+3,
  if s(i)<q, s(i)=0 ; 
  else s(i)=1/s(i) ; end ;
end ;

c=v*diag(s)*u'*[z(:) ; 0 ; 0 ; 0 ] ;