title: Results from my thin plate spline interpolation implementation are dependant of the independent variables
body: <p>I implemented the <a href="https://en.wikipedia.org/wiki/Thin_plate_spline" rel="nofollow noreferrer">thin plate spline</a> algorithm (see also <a href="http://profs.etsmtl.ca/hlombaert/thinplates/" rel="nofollow noreferrer">this description</a>) in order to interpolate scattered data using Python.</p>
  
  <p>My algorithm seems to work correctly when the bounding box of the initial scattered data has an aspect ratio close to 1. However, scaling one of the data points coordinates changes the interpolation result. I created a minimal working example that is representative of what I am trying to accomplish. Below are two plots showing the results of the interpolation of 50 random points.</p>
  
  <p>First, the interpolation of <code>z = x^2</code> on the domain <code>x = [0, 3], y = [0, 120]</code>:</p>
  
  <p><a href="https://i.stack.imgur.com/ePxf9.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ePxf9.png" alt="failed interpolation"></a></p>
  
  <p>As you can see, the interpolation fails. Now, executing the same process but after scaling the <code>x</code> values by a factor of 40, I get:</p>
  
  <p><a href="https://i.stack.imgur.com/7X3ZV.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/7X3ZV.png" alt="successful interpolation"></a></p>
  
  <p>This time, the result looks better. Choosing a slightly different scaling factor would have resulted in a slightly different interpolation. This shows that something is wrong in my algorithm but I can't find what exactly. Here is the algorithm:</p>
  
  <pre><code>import numpy as np
  import numba as nb
  
  # pts1 = Mx2 matrix (original coordinates)
  # z1   = Mx1 column vector (original values)
  # pts2 = Nx2 matrix (interpolation coordinates)
  
  def gen_K(n, pts1):
      K = np.zeros((n,n))
  
      for i in range(0,n):
          for j in range(0,n):
              if i != j:
                  r = ( (pts1[i,0] - pts1[j,0])**2.0 + (pts1[i,1] - pts1[j,1])**2.0 )**0.5
                  K[i,j] = r**2.0*np.log(r)
  
      return K
  
  def compute_z2(m, n, pts1, pts2, coeffs):   
      z2 = np.zeros((m,1))
  
      x_min = np.min(pts1[:,0])
      x_max = np.max(pts1[:,0])
      y_min = np.min(pts1[:,1])
      y_max = np.max(pts1[:,1])
  
      for k in range(0,m):
          pt = pts2[k,:]
  
          # If point is located inside bounding box of pts1
          if (pt[0] &gt;= x_min and pt[0] &lt;= x_max and pt[1] &gt;= y_min and pt[1] &lt;= y_max):
              z2[k,0] = coeffs[-3,0] + coeffs[-2,0]*pts2[k,0] + coeffs[-1,0]*pts2[k,1]
  
              for i in range(0,n):
                  r2 = ( (pts1[i,0] - pts2[k,0])**2.0 + (pts1[i,1] - pts2[k,1])**2.0 )**0.5
  
                  if r2 != 0:
                      z2[k,0] += coeffs[i,0]*( r2**2.0*np.log(r2) )
  
          else:
              z2[k,0] = np.nan
  
      return z2
  
  gen_K_nb = nb.jit(nb.float64[:,:](nb.int64, nb.float64[:,:]), nopython = True)(gen_K)
  compute_z2_nb = nb.jit(nb.float64[:,:](nb.int64, nb.int64, nb.float64[:,:], nb.float64[:,:], nb.float64[:,:]), nopython = True)(compute_z2)
  
  def TPS(pts1, z1, pts2, factor):
      n, m = pts1.shape[0], pts2.shape[0]
  
      P = np.hstack((np.ones((n,1)),pts1))
      Y = np.vstack((z1, np.zeros((3,1))))
  
      K = gen_K_nb(n, pts1)
      K += factor*np.identity(n)
  
      L = np.zeros((n+3,n+3))
      L[0:n, 0:n] = K
      L[0:n, n:n+3] = P
      L[n:n+3, 0:n] = P.T
  
      L_inv = np.linalg.inv(L)
      coeffs = L_inv.dot(Y)
  
      return compute_z2_nb(m, n, pts1, pts2, coeffs)
  </code></pre>
  
  <p>Finally, here is the code snippet I used to create the two plots:</p>
  
  <pre><code>import matplotlib.pyplot as plt
  import numpy as np
  
  N = 50 # Number of random points
  
  pts = np.random.rand(N,2)
  pts[:,0] *= 3.0 # initial x values
  pts[:,1] *= 120.0 # initial y values
  
  z1 = (pts[:,0])**2.0
  
  for scale in [1.0, 40.0]:
  
      pts1 = pts.copy()
      pts1[:,0] *= scale
  
      x2 = np.linspace(np.min(pts1[:,0]), np.max(pts1[:,0]), 40)
      y2 = np.linspace(np.min(pts1[:,1]), np.max(pts1[:,1]), 40)
      x2, y2 = np.meshgrid(x2, y2)
  
      pts2 = np.vstack((x2.flatten(), y2.flatten())).T
      z2 = TPS(pts1, z1.reshape(z1.shape[0], 1), pts2, 0.0)
  
      # Display
      fig = plt.figure(figsize=(4,3))
      ax = fig.add_subplot(111)
      C = ax.contourf(x2, y2, z2.reshape(x2.shape), np.linspace(0,9,10), extend='both')
      ax.plot(pts1[:,0], pts1[:,1], 'ok')
      ax.set_xlabel('x')
      ax.set_ylabel('y')
      plt.colorbar(C, extendfrac=0)
      plt.tight_layout()
  
  plt.show()
  </code></pre>