3

I have a not-quite linear gradient at some angle to the horizontal as an image. Here's some toy data:

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)

diagonal gradient

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.

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.

enter image description here

Any idea how to achieve this? Thanks!

Joe Flip
  • 1,076
  • 4
  • 21
  • 37
  • possible duplicate of http://stackoverflow.com/questions/33085142/skewing-an-array-in-python ? – dnalow Sep 29 '16 at 18:01
  • @dnalow It's close but that solution doesn't interpolate. As I put in the description, it's not quite a linear gradient so there's more than just a skew that needs to be done. I guess each row would need to be interpolated to the values of the bottom row before they are translated in x. – Joe Flip Sep 29 '16 at 18:37

1 Answers1

1

You can intepolate to determine the skewness and interpolate again to correct it. 

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)

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:

demo

Full script:

import numpy as np
from scipy.ndimage.interpolation import map_coordinates

def fix(g):
    x = 1 if g[0,0] < g[0,-1] else -1
    y = 1 if g[0,0] < 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()
user6758673
  • 579
  • 2
  • 4