1

I would like to plot a 1D profile of a 2D image along an arbitrary line. The code below loads the image data hosted on github and plots it:

import urllib
import numpy as np
import matplotlib.pyplot as plt

url = "https://gist.github.com/andreiberceanu/7141843/raw/0b9d50d3d417b1cbe651560470c098700df5a1fc/image.dat"
f = urllib.urlopen(url)
data = np.loadtxt(f)

plt.imshow(data)

waves with line

The red line in the plot above was drawn by hand, as an example. I suppose one can parametrize it in the form a*x + b. I am also guessing some sort of interpolation is necessary, because the line passes though points which may not be part of the original 2D array of data.

Andrei Berceanu
  • 351
  • 2
  • 4
  • 10

1 Answers1

2

You want to use scipy.ndimage.map_coordinates. You need to build up a 2xn array that is the coordinates at which to sample and then do map_coordinates(im, samples).

I think this is it:

def sliceImage(I, a, b, *arg, **kws):
    from scipy import linspace, asarray
    from scipy.ndimage import map_coordinates
    from scipy.linalg import norm
    dst = norm(asarray(b) - a) + 1
    return map_coordinates(I, [linspace(strt, end, dst) 
                               for strt, end in zip(a, b)],
                           *arg, **kws)

Edit: On further consideration, I think this is more elegant:

def sliceImage(I, a, b, *arg, **kws):
    from scipy import linspace, asarray
    from scipy.ndimage import map_coordinates
    from scipy.linalg import norm
    a = asarray(a)
    b = asarray(b)
    dst = norm(b - a) + 1
    return map_coordinates(I, (a[:,newaxis] * linspace(1, 0, dst) +
                               b[:,newaxis] * linspace(0, 1, dst)),
                           *arg, **kws)

Edit: Thanks tcaswell: Added 1 to dst.

Ben
  • 9,184
  • 1
  • 43
  • 56
  • 1
    I feel there are clearer ways to generate the list of points to sample at. – tacaswell Oct 25 '13 at 00:44
  • I agree. There's this approach, which is more numpyish: a = asarray(a); b = asarray(b); dst = norm(b - a); return map_coordinates(I, (a[:,newaxis] * linspace(1, 0, dst) + b[:,newaxis] * linspace(0, 1, dst)), *arg, **kws) – Ben Oct 25 '13 at 02:02
  • 1
    I would use something like `linspace(0, 1, int(np.ceil(dst)) + 1)` to make sure you get enough points. (make sure `a, b =(0, 0), (0, 1)` gives back what you think it should) – tacaswell Oct 25 '13 at 02:42