fast Cartesian to Polar to Cartesian in Python

Question

I want to transform in Python 2d arrays/images to polar, process then, and subsequently transform them back to cartesian. The following is the result from ImajeJ Polar Transformer plugin (used on the concentric circles of the sample code):

The number and dims of the images is quite large so I was checking whether openCV has a fast and simple way to do this.

I read about cv. CartToPolar and PolarToCart but I failed to use it. I understand better the LogPolar where the input and output are arrays, and where you can set the center, interpolation,and inversion (i.e CV_WARP_INVERSE_MAP). Is there a way to use CartToPolar/PolarToCart in an similar fashion?

    import numpy as np
    import cv

    #sample 2D array that featues concentric circles
    circlesArr = np.ndarray((512,512),dtype=np.float32)
    for i in range(10,600,10): cv.Circle(circlesArr,(256,256),i-10,np.random.randint(60,500),thickness=4)

    #logpolar
    lp = np.ndarray((512,512),dtype=np.float32)
    cv.LogPolar(circlesArr,lp,(256,256),100,cv.CV_WARP_FILL_OUTLIERS)

    #logpolar Inverse
    lpinv = np.ndarray((512,512),dtype=np.float32)
    cv.LogPolar(lp,lpinv,(256,256),100, cv.CV_WARP_INVERSE_MAP + cv.CV_WARP_FILL_OUTLIERS)

    #display images
    from scipy.misc import toimage
    toimage(lp, mode="L").show()
    toimage(lpinv, mode="L").show()

This is for a tomography (CT) workflow where rings artifacts can be filtered out easier if they appear as lines.

Solution 1

the CV source code mentions a LinearPolar. it doesn't seem to be documented, but appears to be similar to LogPolar. have you tried that?

Solution 2

Latest versions of opencv supports a function cv2.linearPolar. This may be another solution that does not involve the use of opencv:

def polar2cart(r, theta, center):

    x = r  * np.cos(theta) + center[0]
    y = r  * np.sin(theta) + center[1]
    return x, y

def img2polar(img, center, final_radius, initial_radius = None, phase_width = 3000):

    if initial_radius is None:
        initial_radius = 0

    theta , R = np.meshgrid(np.linspace(0, 2*np.pi, phase_width), 
                            np.arange(initial_radius, final_radius))

    Xcart, Ycart = polar2cart(R, theta, center)

    Xcart = Xcart.astype(int)
    Ycart = Ycart.astype(int)

    if img.ndim ==3:
        polar_img = img[Ycart,Xcart,:]
        polar_img = np.reshape(polar_img,(final_radius-initial_radius,phase_width,3))
    else:
        polar_img = img[Ycart,Xcart]
        polar_img = np.reshape(polar_img,(final_radius-initial_radius,phase_width))

    return polar_img

Solution 3

Here's an example of the log-polar transform implemented using SciPy:

https://github.com/stefanv/supreme/blob/master/supreme/transform/transform.py#L51

Given that this is only a coordinate transformation, it should be easier to adapt to your problem than the OpenCV version.

Comments

[Papado]

Thank you very very much! Indeed LinearPolar does what it says. Unfortunately by using import cv it was not available, but I tried from opencv import cv and then cv.cvLinearPolar and works. Next days I'll try it's performance in large datasets. Thank you!

[andrew cooke]

cool. i wonder why it's not visible? i'll try filing a bug report.

[Papado]

Dear Stefan, thank you very much for you feedback. I'll check and benchmark your implementation the next days. Btw, I ended up browsing Supreme and seems very interesting. Have you published any article about it?

[cjm2671]

Hi Alessandro- I tried to solve this problem in this way too, and I wrote similar code, although I used a loop instead of meshgrid (I've never seen that before); do you know how well this performs? Mine was taking on the order of 1s for a VGA image - too long.

[cjm2671]

OK, I have tested your code and it is VERY fast compared with my iterative solution- plus I learned something new- thank you so much!

[Stefan van der Walt]

@Papado I never saw your comment, but yes--there's a paper on arXiv and a dissertation. By the way, the log polar transform can now be implemented on top of scikit-image in about 5 lines of code, using skimage.transform.warp.