Converting an image from polar to cartesian coordinates

Question

I have the following code which has too many loops inside it (almost 13 billion). This code is actually for image processing. I am first making a matrix of dimension 1024*360. This is the value of the projection data given to me. f is my image in the polar coordinate. The first loop gives me f in the polar-coordinate. f[r,phi] where r varies from 0 to 1024 and phi varies from 0 to 90. h is a filter that I have used in the first loop. In the second I am converting the image from polar coordinate to Cartesian coordinate. The image in the Cartesian coordinate is of dimension 725*725. I have set the range of x and y as (1,726).

I doubt that it is right too. But when I had it as (0,725) earlier, it showed me error "dividing with zero" as I use atany(y/x). Hence I changed it. And after that I am converting the image to grayscale. I am new to Python and I guess my algorithm seems to slow down the execution. Is there any other way to avoid these loops? I haven't been able to display the image until now.

from __future__ import division
import scipy.io 
import numpy as np
import pylab as pl
import matplotlib.pyplot as plt
import math
from math import sin
from math import pow
from math import degrees
from math import atan
import Image

q = np.zeros((1024,360))

x = scipy.io.loadmat('/home/prakriti/Project/A.mat')  
A = np.array((x.values()))
q[:,0:120] = A[0][240,:,:]
del x
del A

y = scipy.io.loadmat('/home/prakriti/Project/B.mat')  
B = np.array((y.values()))
q[:,120:240] = B[0][240,:,:]
del y
del B

z = scipy.io.loadmat('/home/prakriti/Project/C.mat')  
C = np.array((z.values()))
q[:,240:360] = C[0][240,:,:]
del z
del C

f = np.zeros((1024,90))  #In polar co-ordinate 

D = 0.28                #Source to isocenter distance
h = np.zeros(1024)
h[0:512] = np.linspace(0,1,512)
h[513:] = np.linspace(1,0,511) 

for r in range(0,1024) :
  for phi in range(0,90) :
    for b in range(0,360) :
      for s in range(0,1024) :
        U = (D + r*sin(b-phi))/D
        l = math.pow(U,-2)
        k = D/(math.pow((math.pow(D,2)+math.pow(s,2)),0.5))   
        f[r,phi] = 0.5*l*k*q[s,b]*h[s]

I =np.zeros((725,725))

for x in range(1,726) :
  for y in range(1,726) :
    r = math.trunc(math.pow((math.pow(x,2)+math.pow(y,2)),0.5))
    phi = math.trunc(math.degrees(math.atan(y/x)))
    I[x,y] = f[r,phi]

I8 = (((I - I.min()) / (I.max() - I.min())) * 255.9).astype(np.uint8)
img = Image.fromarray(I8)
img.save("Fanbeamreconstruction.png")
im = Image.open("Fanbeamreconstruction.png")
im.show()

Someone also pointed out that may be in looping over b and s I am overwriting f[r,phi] 360*1024 times and only the last one is used. What could this possible mean?

Answer 1

Well for starters, you can use xrange instead of range which will not create an entire list before you iterate, just create values as you need them. You can also move parts of the calculation to upper levels in your loops to decrease the work being done every cycle.

D_pow2 = math.pow(D,2)

for r in xrange(0,1024):
    for phi in xrange(0,90):
        for b in xrange(0,360):
            U = (D + r*sin(b-phi))/D
            l = math.pow(U,-2)

            for s in xrange(0,1024):
                k = D/(math.pow(D_pow2 + math.pow(s,2)),0.5))   
                f[r,phi] = 0.5*l*k*q[s,b]*h[s]