polar_transform.py

# http://stackoverflow.com/questions/9924135/fast-cartesian-to-polar-to-cartesian-in-python

import cv2
import numpy as np
import scipy.ndimage as nd

import irtk

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, radius, theta_step=300,zoom=2):
    theta, R = np.meshgrid( np.linspace(0, 2*np.pi, theta_step),  # x
                            np.linspace(0, radius,zoom*radius))   # y

    print "theta-R shape", theta.shape, R.shape

    Xcart, Ycart = polar2cart(R, theta, center)
    polar_img = nd.map_coordinates(img, [Ycart,Xcart], order=3, mode='nearest')
    polar_img = np.reshape(polar_img, (zoom*radius,theta_step))

    print "polar_img shape", polar_img.shape
    
    return polar_img

def spherical2cart(r, theta, phi, center):
    """ theta in [0,2pi] and phi in [0,pi]"""
    x = r * np.sin(phi) * np.cos(theta) + center[0]
    y = r * np.sin(phi) * np.sin(theta) + center[1]
    z = r * np.cos(phi) + center[2]
    
    return x, y, z

def img2spherical(img, center, radius, theta_step=200, phi_step=100, zoom=2, order=3):
    # see
    # http://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html
    # for indexing='ij'
    R, theta, phi = np.meshgrid( np.linspace(0, radius,zoom*radius),
                                 np.linspace(0, 2*np.pi, theta_step),
                                 np.linspace(0, np.pi, phi_step),      
                                 indexing='ij' )                 

    print "theta-R shape", theta.shape, R.shape
    
    Xcart, Ycart, Zcart = spherical2cart(R, theta, phi, center)
    polar_img = nd.map_coordinates(img, [Zcart,Ycart,Xcart], order=order, mode='constant')
    polar_img = np.reshape(polar_img, (zoom*radius,theta_step,phi_step))

    print "polar_img shape", polar_img.shape
    
    return polar_img

f = "5066_4.nii"
brain_center = [-10.1,75.3,-52.8]

img = irtk.imread(f, dtype='float32').saturate().rescale()

img1 = img[42]
brain_center1 = img1.WorldToImage(brain_center)

center = brain_center1[:2]
radius = img1.shape[0]/2
img_polar = img2polar(img1, center, radius, theta_step=300)

cv2.imwrite("test.png", img_polar)

img = img.resample(1.0,interpolation='bspline')

print img.shape

center = img.WorldToImage(brain_center)
radius = img.shape[2]/2
img_polar = img2spherical(img, center, radius, theta_step=400, phi_step=200,order=1)

irtk.imwrite("res.nii.gz", img_polar)

print img_polar.shape