# Fast algorithm for interpolating data from polar coordinates to cartesian coordinates

I have a set of solution nodes generated over a polar grid. I would like to convert / interpolate these solution nodes onto a Cartesian grid:

That is, using the image above, for each node in the Cartesian grid I would interpolate a value from the closest existing nodes (red).

Currently, my approach is to generate a kd-tree for the original solution nodes, then use a nearest-neighbor search to obtain the three closest nodes. I then use barycentric interpolation to obtain a value from these three points. The problem, however, is that my polar grid is much finer along the radial direction than it is in the azimuthal direction, which means that my nearest-neighbor search almost always selects points from the same radial. This has the result of creating "striations" in my new solution, instead of smoothly interpolating along the azimuthal direction (i.e., the results look no different than if I had simply mapped the nearest point to the "interpolated" point).

Unfortunately, I don't know how to achieve a better sampling without sacrificing the kd-tree and losing a lot of the speed improvements. Am I being thick-headed and missing an obvious solution? Or does anyone know a better way to approach this problem?

• I can't tell what you're trying to do. To "interpolate" is to approximate the value of a function within some interval based on knowledge of what it does outside that interval. Your question doesn't seem to be about interpolation. Are you trying to find the nearest Cartesian grid point to each of your points given in polar co-ordinates? If so, what's wrong with just taking $x=r\cos\theta$, $y=r\sin\theta$ and rounding to the nearest grid point? Jun 19, 2015 at 7:34
• @David, I'm trying to map results from one mesh onto another. Each $(r,\theta)$ node has a solution $\mathbf{\phi}$ associated with it. I'm trying to re-represent these results onto a Cartesian mesh, such that I now have results $\mathbf{\phi'}$ for each pre-determined $(x, y)$ node, where $(x, y, \mathbf{\phi'})$ is interpolated from surrounding $(r,\theta,\mathbf{\phi})$ nodes. Jun 19, 2015 at 15:48
• So you have a function $f(r,\theta)$ whose values you know at the points of a polar grid and you want to interpolate to find the values of the function $g(x,y) = f(\sqrt{x^2+y^2}, \tan^{-1}(y/x))$ at the points of a Cartesian grid? (Modulo any mistakes I just made in converting from Cartesian to polar co-ordinates. ;-) ) Jun 19, 2015 at 16:03
• It certainly would be trivial to convert the Cartesian coordinates to polar and run them through the solver, but this is a computationally heavy process that we won't have the time for in the environment these models will be deployed in. Instead, one alternative is to precompute results for different sets of conditions, and then simply call those results and map them to the solution domain when the real-time conditions match the conditions of the precomputed solutions. Due to historical reasons, the precomputed solutions are generated over a polar grid, but the real-time grid is Cartesian. Jun 20, 2015 at 1:14

I don't know if this has been studied by others. Speculating off the top of my head, let me outline two possible approaches.

## Applying standard rectangular methods to polar coordinates

Here's one way to adapt more familiar interpolation methods (that work with Cartesian coordinates) to your setting. Let's start by looking at bilinear interpolation.

Let $P$ be the point on the Cartesian grid whose value you want to compute, using interpolation. Find the smallest "rectangle" of four nodes that surround the point $P$. A "rectangle" in polar coordinates is just like a rectangle in a Cartesian grid, except with up/down replaced with "inward"/"outward" and right/right replaced with "clockwise"/"counterclockwise". For instance, here is an example of a "rectangle" in a polar-coordinate grid:

This gives you 4 nodes that surround $P$. Finally, use these 4 nodes to interpolate and form the value for $P$. For instance, you could do bilinear interpolation, in the same way we do bilinear interpolation on a Cartesian grid, but now working with polar coordinates / barycentric coordinates.

This also makes clear how to generalize to other interpolation schemes, such as bicubic interpolation.

One possible approach you could try would be to pick the 3 closest nodes, then supplement them with a few other nodes from adjacent radials.

Let $P$ be your point. Start by picking the closest 3 nodes to $P$ (say). Then, supplement them. For instance, suppose all 3 nodes are on the radial $\theta$. Look at the radial just to the left of it, $\theta-\epsilon$, and out of all nodes on that radial, pick the 1 or 2 nodes that are closest to $P$. Similarly, look at the radial just to the right of it, $\theta+\epsilon$, and out of all nodes on that radial, pick the 1 or 2 nodes that are closest to $P$.

Now use those points to interpolate. You might use weighted least-squares regression to fit a linear (or bilinear/bicubic/etc.) model, where the nodes that are closer to $P$ receive a higher weight. If you experiment a bit with the number of points chosen, the type of model, and the weights, you might find a scheme that reduces the "striation" artifacts.

Maybe better schemes are known; I don't know. I'm not expert in this area.

I had the same problem when I was trying to simulate the Projection-Slice-Theorem.

A quick solution is to use generalized bivariate interpolations, e.g. from random sample points to rectangular grid (not specific to your problem , but will work, example in python here https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata.

For my problem which involved taking an FFT of a function sampled in cylindrical (rather polar) coordinates to transform to a Cartesian xy grid, there is a very efficient solution published by Averbuch, Donoho, Coifman and Elad in 2006: "Fast and Accurate Polar Fourier Transform". (https://elad.cs.technion.ac.il/software/ scroll down to Polar-FFT to get both paper and free code in Matlab).