# Methods for interpolating from points that are not on a regular grid?

I'm working on a project where I'll need to be able to interpolate scalar potential values at arbitrary points in a 3-D box from a large (potentially millions to billions of points) collection of values that were calculated at points of convenience (roughly randomly and uniformly distributed) that are not on a regular grid. I'll be computing a trajectory and interpolating values of this potential at points along along the trajectory, so my queries will consistently be at points near the most recently used point.

Without much background in this area, my initial idea would be to put the known points into a k-d tree and access the tree each time I need to interpolate a value. An alternative would be to sort and bin the points and then use them to interpolate values onto a rectangular grid, and then interpolate from the regular grid as needed.

Are there other more specialized data structures and interpolation methods that might be useful for this task?

• As method RBF can be used for interpolating scattered data. Maybe is a suitable option for your problem. – Mauro Vanzetto Dec 2 '16 at 8:19
• RBF will be fine for the interpolation once I've identified a set of nearby points- the problem is finding the nearby points. – Brian Borchers Dec 2 '16 at 15:23
• As you suggested, kd-tree can be both fast and efficient to find the nearest points. You can implement it with a very small memory cost if you reorder the points (then you no longer need to keep the reference). In addition if you do a balanced kd-tree, you further save memory. There is an implementation in my geogram lib: alice.loria.fr/software/geogram/doc/html/index.html – BrunoLevy Dec 2 '16 at 18:22
• I also recommend a balanced KD-Tree approach to identify neighbors to some location you want to find a value at. I would recommend then using something like local regression and building a local linear model with your nearest neighbors dataset. This tends to provide fairly accurate results, relative to some typical least square fit, handles unstructured data easily, and in 3D would be quite efficient. – spektr Dec 2 '16 at 19:32
• Since you have a roughly uniform distribution of points, a $kd$ tree seems like overkill. Just sort the points into the cells of a uniform mesh of size, say, $(N^{1/3}/3)^3$ if you have $N$ points. Then every one of these cells has $3^3=27$ points, which you can easily search through for nearest neighbors. Furthermore, it is an $O(1)$ effort to identify the cell your trajectory is current in, along with the neigbhors of the current evaluation point if you need more interpolation points. – Wolfgang Bangerth Dec 3 '16 at 1:45

I have some experience working in particle trajectory simulations using interpolated flow field information. In our case, we needed to simulate the motion of numerous particle using a single time dependent flow. Consequently, we had to carry out a complex interpolation (position and time) during the trajectories.

Our conclusion ended up being that it was much more cost effective to map the flow information on a Cartesian grid. Then, using integer arithmetic and the particle position, we could quickly find the position of the particle within the grid. Afterward, we would use trilinear interpolation in order to obtain the flow properties at the position of the particle. This was by far the most robust and efficient approach. Furthermore (and this is important I find) it could be parallelized in distributed memory environment (MPI) quite easily since you only have to decompose the Cartesian grid (and ensure the correct transfer of information for the particles). Since we had to load a considerable amount of data, distributed memory parallelism was essential.

Summary of why we went to a Cartesian mapping:

• Mapping is expensive, but you only need to do it once per data
• Allows for easy distributed memory parallelism
• Allowed us to use trivial trilinear interpolation on the grid

Some issues you might encounter:

• Your data structure could become very sparse if you have a complex 3D geometry.