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[240,:,:] del x del A y = scipy.io.loadmat('/home/prakriti/Project/B.mat') B = np.array((y.values())) q[:,120:240] = B[240,:,:] del y del B z = scipy.io.loadmat('/home/prakriti/Project/C.mat') C = np.array((z.values())) q[:,240:360] = C[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?