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

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)

StackExchange Code Review Q#79381, answer score: 7