# Remapping and interpolating 255x1024 array

I have a function that accepts a $255\times1024$ array and remaps it to another array in order to account for some hardware related distortion (lines that should be straight are curved by the lens). At the moment it does exactly what I want, but slowly (roughly 30 second runtime).

Specifically, I'm looking at the nested for loops that take 18 seconds to run, and the interpolation that takes 10s. Is there any way to optimize/speed up this process?

EDIT: Nested for loops have been optimized as per vps' answer. Am now only interested in optimizing the interpolation function (if that's even possible).

def smile(Z):
p2p = np.poly1d([ -3.08049538e-07,   3.61724996e-04,  -7.78775408e-02, 3.36876203e+00])
Y = np.flipud(np.rot90(np.tile(np.linspace(1,255,255),(1024,1))))
X = np.tile(np.linspace(1,1024,1024),(255,1))
for m in range(0,255):
for n in range(0,1024):
X[m,n] = X[m,n] - p2p(m+1)
x = X.flatten()
y = Y.flatten()
z = Z.flatten()
xy = np.vstack((x,y)).T
grid_x, grid_y = np.mgrid[1:1024:1024j, 1:255:255j]
newgrid = interpolate.griddata(xy, z,(grid_x,grid_y), method = 'linear',fill_value = 0).T
return newgrid

• Exemplary way for a first question to be asked by a new user. +1, I really hope you receive some quality reviews. Feb 2 '15 at 19:45

p2p(m + 1) does not depend on n. You may safely extract its calculation from the inner loop:

for m in range(0,255):
p2p_value = p2p(m + 1)
for n in range(0,1024):
X[m,n] = X[m,n] - p2p_value


If you call smile() multiple times, it is worthwhile to precompute p2p once.

Some further savings can be achieved by accounting for the sequentiality of p2p arguments. Notice that p2p(m+1) - p2p(m) is a polynomial of lesser degree; you may calculate it incrementally:

    p2p_value += p2p_delta(m)


Edit:

Some math:

Let $P(x) = ax^3 + bx^2 + cx + d$. You may see that $P(x+1) - P(x) = 3ax^2 + 3ax + a + 2bx + b + c = 3ax^2 + (3a + 2b)x +(a+b+c)$ is a second degree polynomial, a bit easier to calculate than the original third degree one. Which leads to the code (fill up the list of coefficients according to the above formula):

    p2p_delta = np.poly1d([...])
p2p_value = p2p_delta(0)
for m in range(0, 255)
for n in range(0, 1024)
X[m,n] -= p2p_value
p2p_value += p2p_delta(m)

• When you say precompute, do you mean make a list of all the p2p values and then do X[m,n] = X[m,n] - p2p_value[m+1]? Also I'm afraid I don't quite understand what's going on in the latter part of your answer, could you elaborate? Thanks. Feb 2 '15 at 21:01
• Re precompute, yes, I mean exactly that. Re deltas, will try to elaborate.
– vnp
Feb 2 '15 at 21:13
• Ah, gotcha. Thanks for the help. I managed to cut the runtime on the nested for loops to <0.4 seconds, which is absolutely fantastic. Although apparently doing the p2p(m+1) calculation is actually faster than indexing through p2p_values[m+1] (I'm not even including the time to initialize the p2p_values list)... I guess python doesn't like indexing? Lastly, I'm just going to assume that the interpolate.griddata(...) method can't be optimized, being a built-in function and all. Thanks again! Feb 2 '15 at 21:49