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
  • 3
    \$\begingroup\$ Exemplary way for a first question to be asked by a new user. +1, I really hope you receive some quality reviews. \$\endgroup\$
    – syb0rg
    Commented Feb 2, 2015 at 19:45

1 Answer 1


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)


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)
  • \$\begingroup\$ 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. \$\endgroup\$
    – wes3449
    Commented Feb 2, 2015 at 21:01
  • \$\begingroup\$ Re precompute, yes, I mean exactly that. Re deltas, will try to elaborate. \$\endgroup\$
    – vnp
    Commented Feb 2, 2015 at 21:13
  • \$\begingroup\$ 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! \$\endgroup\$
    – wes3449
    Commented Feb 2, 2015 at 21:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.