brokenseas brokenseas - 28 days ago 19
Python Question

Calculating a 3D gradient with unevenly spaced points

I currently have a volume spanned by a few million every unevenly spaced particles and each particle has an attribute (potential, for those who are curious) that I want to calculate the local force (acceleration) for.

np.gradient only works with evenly spaced data and I looked here: Second order gradient in numpy where interpolation is necessary but I could not find a 3D spline implementation in Numpy.

Some code that will produce representative data:

import numpy as np
from scipy.spatial import cKDTree

x = np.random.uniform(-10, 10, 10000)
y = np.random.uniform(-10, 10, 10000)
z = np.random.uniform(-10, 10, 10000)
phi = np.random.uniform(-10**9, 0, 10000)

kdtree = cKDTree(np.c_[x,y,z])
_, index = kdtree.query([0,0,0], 32) #find 32 nearest particles to the origin
#find the gradient at (0,0,0) by considering the 32 nearest particles?


(My problem is very similar to Compute 3D gradient with unevenly spaced sample locations but there seemed to have been no solution there, so I thought I'd ask again.)

Any help would be appreciated.

lhk lhk
Answer

Intuitively, for the derivate wrt one datapoint, I would do the following

  • Take a slice of the surrounding data: data=phi[x_id-1:x_id+1, y_id-1:y_id+1, z_id-1:z_id+1]. The approach with the kdTre looks very nice, of course you can use that for a subset of the data, too.
  • Fit a 3D polynomial, you might want to look at polyvander3D. Define the point in the middle of the slice as the center. Calculate the offsets to the other points. Pass them as coordinates to the polyfit.
  • Derive the polynomial at your position.

This would be a simple solution to your problem. However it would probably be very slow.

EDIT:

In fact this seems to be the usual method: http://scicomp.stackexchange.com/questions/480/how-can-i-numerically-differentiate-an-unevenly-sampled-function

The accepted answer talks about deriving an interpolating polynomial. Although apparently that polynomial should cover all the data (Vandermonde matrix). For you that is impossible, too much data. Taking a local subset seems very reasonable.