Fast algorithm for polar -> cartesian conversion

Question

I have an image on a polar grid. This image should be transformed into a cartesian grid, but the only algorithm I know of is really slow for this. Now I use the cartesian grid, for each point I find the r and theta values, and then I look in two vectors to find the smallest error defined by:

min{(th_vec - theta)^2 + (range - r)^2}

This gives a nested for-loop inside of the outer nested for-loop, so I have a complexity of O(N^4). A 512x512 image uses a whole minute to complete. Of course, a complexity like that can not be used, so I'm wondering if anyone know of any faster algorithms to do this?

I have the image, and the two vectors. The X-axis of the image is the angle, while the Y-axis of the image is the length from the center. The angle is always from 0-2pi, and the range goes from 0 to r_max.

Thank you in advance.

EDIT: The range goes from 0 to r_max, not -r_max to r_max as it stood before. I see that there have been some missunderstandings. I have used the normal, inverse, conversion with;

r=sqrt(x^2 + y^2);
theta=atan2(y,x);

The problem is that I have to first convert the x and y values to x' and y' values, since the grid is from -r_max to r_max in the resulting image, but in pixels in the data. So I have a 512x512 image, but r_max can be something like 3.512. So I have to convert each pixel value into the grid value, then find the r and theta values. When I have found the r and theta values I have to run trough two vectors, range and th_vec, to find the pixel in the original image that matches:

min{(range - r)^2 + (th_vec - theta)^2}

This gives me a complexity of O(n^4), since the th_vec and range vectors are the same size as the image. So if I have a square matrix of 512x512 elements, I have to run trough 68 719 476 736 elements, which is way slow. So I'm wondering if there is a faster algorithm? I can't change the input data, so as far as I know, this is the only way to do it if you don't start with triangulation and stuff, but this is to expensive in times of memory.

Spillage

What's this for? Also, why don't you have either angle from 0 to pi or range from 0 to r_max? 2*pi gives a complete circle, so why would you need negative distance? – Solemnize

Is your polar grid uniformly partitioned with respect to the polar coordinates? – Zemstvo

If you find r_0 and th_0 as some floating point value from your x,y then you only have to look at four pairs (r,th) in your polar image, i.e. the four nearest neighbours of (r_0,th_0), so the four combinations of floor(r_0),ceil(r_0) and floor(th_0),ceil(th_0) where floor() and ceil() produce something that is rounded to your polar grid. – Beneficiary

Answer 1

How about

x=r*cos(angle)
y=r*sin(angle)

This is the standard way of converting from polar to Cartesian, and unless you're going to use some kind of table lookup, there's not really a faster option.

Edit: wrang wrang has a good point. If you're trying to transform an image in polar coordinates I(angle, r) to an image in Cartesian coordinates I_new(x, y), you are definitely better off using the inverse transformation, as follows:

for x=1,...,width
    for y=1,...,height
        angle=atan2(y, x)
        r=sqrt(x^2+y^2)
        I_new(x, y)=I(angle, r)
    end
end

As a rule, angle and r will not be integer, so you have to do some kind of interpolation in the image I. The easiest way to do this is simply to round angle and r; this will give you nearest-neighbour interpolation. If you need better quality, try more sophisticated types of interpolation such as bilinear or bicubic interpolation.

Uranic

You have to describe how to fill in the entire image, so the inverse transform is more likely useful unless you're using some clever interpolation method. – Solemnize

Answer 2

You could loop over each pixel in the polar image map and then render the resulting arc section in the cartesian image plane:

polar to cartesian conversion http://img24.imageshack.us/img24/4635/polartocartesian.png

const float dR = 2*r_max / polar_image_height;
const float dA = 2*pi / polar_image_width;

float angle;
float radius;
for (int polar_x = 0; polar_x < polar_image_width; polar_x++)
{
    for (int polar_y = 0; polar_y < polar_image_height; polar_y++)
    {
        angle = polar_x * dA;
        radius = polar_y * dR - r_max;
        DrawArcSection(radius, radius+dR, angle, angle+dA);
    }
}

Many drawing libraries have built-in functions for drawing that arc section, but you could always just approximate it with a simple polygon:

void DrawArcSection(float minRadius, float maxRadius,
                    float minAngle, float maxAngle)
{
    point P1 = MakePoint(minRadius * cos(minAngle) + image_width/2,
                         minRadius * sin(minAngle) + image_height/2);
    point P2 = MakePoint(minRadius * cos(maxAngle) + image_width/2,
                         minRadius * sin(maxAngle) + image_height/2);
    point P3 = MakePoint(maxRadius * cos(minAngle) + image_width/2,
                         maxRadius * sin(minAngle) + image_height/2);
    point P3 = MakePoint(maxRadius * cos(maxAngle) + image_width/2,
                         maxRadius * sin(maxAngle) + image_height/2);

    DrawPolygon(P1, P2, P3, P4);
}

Cauliflower

Answer 3

If you don't care about smoothing, why don't you just compute the polar coordinate for each destination Cartesian pixel coordinate and read the color value? See http://en.wikipedia.org/wiki/Polar_coordinate_system#Converting_between_polar_and_Cartesian_coordinates if you need help doing that.

Solemnize

yes - If you want the final image to be filled it's better to do the reverse transform and interpolate in the source image. – Cult

Answer 4

If your grid is uniformly partitioned with respect to the polar coordinates then your algorithm can be reduced to O(N^2) if you take advantage of the fact that the closest point to (r, theta) will be one of the four corners of the grid element in which it is contained.

In the more general case where the grid is the product of arbitrary partitions of the r and theta dimensions, that might grow to O( (N log N)^2) if you have to search for the location of the point in each partition. However, if the partitions were systematically constructed, you should be able to get back down to O(N^2).

Zemstvo

Answer 5

If all your images are 512x512 then I'd use a lookup table that maps a weighted set of pixels in your polar image to the cartesian image. This is a lot of work upfront but makes your final calculation O(n^2). If a LUT is not an option then I'd use:

x=r*cos(angle)
y=r*sin(angle)

On each pixel in the polar image to map it to "a" pixel in the cartesian image, where the output pixel is the average of all input pixels that fall on it. Then apply repeated dilations until there are no uninitialised pixels left. For the dilation you use a 3x3 structuring element and only replace the value of the output pixel with the value of the centre pixel if it previously had no value. Then, as a final measure, apply a gaussian filter to the entire image to smooth out the hard edges. This is the fastest method I can think of that will produce an image that is pleasant to look at in a reasonable amount of time.

Beneficiary

Answer 6

Memory fails, but there might be a fast version of this algorithm that involves the FFT. Once upon a time I took a class on medical imaging and it seems like this kind of stuff popped up when untransforming/rasterizing CT scans. Some keywords to search for would be the radon transform, the filtered backprojection algorithm and CT scans. I looked these up briefly on wikipedia and nothing jumped out, but maybe a more thorough review would turn up some gold.

Pice

Answer 7

O(N²log(N)) algorithm:

• Array S will be used for closest source (polar) coord per cartesian coord.
• S starts filled with a "not-initialized yet" value. (Python: None, Haskell: Nothing, etc.)
• O(N²) - Iterate all your polar coordinates.
  • Translate to cartesian coord
  • Find closest cartesian coord in your destination image. (by rounding and applying borders)
  • Fill in the relevant cell in S with this coordinate
• O(N²log(N)) - Perform a modified Dijkstra algorithm as described below:
  • The "Graph" for our search algorithm is as follows:
    • All cells of S are nodes
    • A cell's neighbors are those the King in chess can move to from it
  • The "Score" of a cell is infinite if it isn't initialized, and distance from the untouched cartesian coord of the polar coord its pointing to
  • When updating a neighbor of cell N we put the value from cell N in it (but like in Dijkstra, only if it makes its score better than its current score)
  • The starting point is the array S initialized as described above

Supercharger