Thanks but I just figured out the problem. It was a comma...all day for a comma.
The problem was with this line tmp = tmp + a[j]*G([xout[i],yout[i]],X[j,:]) It should of been this tmp = tmp + a[j]*G([xout[i] yout[i]],X[j,:]) so when the G function did norm(x1 .- x2) it was silently broadcasting to a 2x2 matrix which is wrong. I have updated the notebook <http://nbviewer.ipython.org/urls/gist.githubusercontent.com/lstagner/89964e63557ddb831e72/raw/05f9ef40c09854fd83a2267bf47c15f286d2c649/tps.ipynb> to reflect the fix and make it a bit more efficient. Thin plate splines are just a special case of Polyharmonic Splines. Would there be interest of expanding this script into a package? -Luke On Wednesday, May 13, 2015 at 12:26:12 AM UTC-7, Jan Kybic wrote: > > 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 >
