Skew a diagonal gradient to be vertical
<p>I have a not-quite linear gradient at some angle to the horizontal as an image. Here's some toy data:</p>

<pre><code>g = np.ones((5,20))
for x in range(g.shape[0]):
    for y in range(g.shape[1]):
        g[x,y] += (x+y)*0.1+(y*0.01)
</code></pre>

<p><a href="http://i.stack.imgur.com/ytcp6.png" rel="nofollow"><img src="http://i.stack.imgur.com/ytcp6.png" alt="diagonal gradient"></a></p>

<p>I want to essentially correct the skew in the gradient so that it is horizontal, i.e. the gradient increases to the right and all vertical slices are constant.</p>

<p>This will of course produce a parallelogram with a larger x-axis than the input image. Returning a masked Numpy array would be ideal. Here's a (terrible) cartoon to quickly illustrate.</p>

<p><a href="http://i.stack.imgur.com/QVtG7.png" rel="nofollow"><img src="http://i.stack.imgur.com/QVtG7.png" alt="enter image description here"></a></p>

<p>Any idea how to achieve this? Thanks!</p>

<p>You can intepolate to determine the skewness and interpolate again to correct it. </p>

<pre><code>import numpy as np
from scipy.ndimage.interpolation import map_coordinates

m, n = g.shape
j_shift = np.interp(g[:,0], g[0,:], np.arange(n))
pad = int(np.max(j_shift))
i, j = np.indices((m, n + pad))
z = map_coordinates(g, [i, j - j_shift[:,None]], cval=np.nan)
</code></pre>

<p>This works on the example image, but you have to do some additional checks to make it function on other gradients. It does not work on gradients that are nonlinear in the x-direction though. Demo:</p>

<p><a href="http://i.stack.imgur.com/hYyAb.png" rel="nofollow"><img src="http://i.stack.imgur.com/hYyAb.png" alt="demo"></a></p>

<p>Full script:</p>

<pre><code>import numpy as np
from scipy.ndimage.interpolation import map_coordinates

def fix(g):
    x = 1 if g[0,0] &lt; g[0,-1] else -1
    y = 1 if g[0,0] &lt; g[-1,0] else -1
    g = g[::y,::x]

    m, n = g.shape
    j_shift = np.interp(g[:,0], g[0,:], np.arange(n))
    pad = int(np.max(j_shift))
    i, j = np.indices((m, n + pad))
    z = map_coordinates(g, [i, j - j_shift[:,None]], cval=np.nan)

    return z[::y,::x]

import matplotlib.pyplot as plt

i, j = np.indices((50,100))
g = 0.01*i**2 + j

plt.figure(figsize=(6,5))
plt.subplot(211)
plt.imshow(g[::-1], interpolation='none')
plt.title('original')
plt.subplot(212)
plt.imshow(fix(g[::-1]), interpolation='none')
plt.title('fixed')
plt.tight_layout()
</code></pre>

python arrays numpy transform