How to extract an arbitrary line of values from a numpy array?

Question

I have a numpy array that contains some image data. I would like to plot the 'profile' of a transect drawn across the image. The simplest case is a profile running parallel to the edge of the image, so if the image array is imdat, then the profile at a selected point (r,c) is simply imdat[r] (horizontal) or imdat[:,c] (vertical).

Now, I want to take as input two points (r1,c1) and (r2,c2), both lying inside imdat. I would like to plot the profile of the values along the line connecting these two points.

What is the best way to get values from a numpy array, along such a line? More generally, along a path/polygon?

I have used slicing and indexing before, but I can't seem to arrive at an elegant solution for such a where consecutive slice elements are not in the same row or column. Thanks for your help.

Solution 1

@Sven's answer is the easy way, but it's rather inefficient for large arrays. If you're dealing with a relatively small array, you won't notice the difference, if you're wanting a profile from a large (e.g. >50 MB) you may want to try a couple of other approaches. You'll need to work in "pixel" coordinates for these, though, so there's an extra layer of complexity.

There are two more memory-efficient ways. 1) use scipy.ndimage.map_coordinates if you need bilinear or cubic interpolation. 2) if you just want nearest neighbor sampling, then just use indexing directly.

As an example of the first:

import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt

#-- Generate some data...
x, y = np.mgrid[-5:5:0.1, -5:5:0.1]
z = np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

#-- Extract the line...
# Make a line with "num" points...
x0, y0 = 5, 4.5 # These are in _pixel_ coordinates!!
x1, y1 = 60, 75
num = 1000
x, y = np.linspace(x0, x1, num), np.linspace(y0, y1, num)

# Extract the values along the line, using cubic interpolation
zi = scipy.ndimage.map_coordinates(z, np.vstack((x,y)))

#-- Plot...
fig, axes = plt.subplots(nrows=2)
axes[0].imshow(z)
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')

axes[1].plot(zi)

plt.show()

The equivalent using nearest-neighbor interpolation would look something like this:

import numpy as np
import matplotlib.pyplot as plt

#-- Generate some data...
x, y = np.mgrid[-5:5:0.1, -5:5:0.1]
z = np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

#-- Extract the line...
# Make a line with "num" points...
x0, y0 = 5, 4.5 # These are in _pixel_ coordinates!!
x1, y1 = 60, 75
num = 1000
x, y = np.linspace(x0, x1, num), np.linspace(y0, y1, num)

# Extract the values along the line
zi = z[x.astype(np.int), y.astype(np.int)]

#-- Plot...
fig, axes = plt.subplots(nrows=2)
axes[0].imshow(z)
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')

axes[1].plot(zi)

plt.show()

However, if you're using nearest-neighbor, you probably would only want samples at each pixel, so you'd probably do something more like this, instead...

import numpy as np
import matplotlib.pyplot as plt

#-- Generate some data...
x, y = np.mgrid[-5:5:0.1, -5:5:0.1]
z = np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

#-- Extract the line...
# Make a line with "num" points...
x0, y0 = 5, 4.5 # These are in _pixel_ coordinates!!
x1, y1 = 60, 75
length = int(np.hypot(x1-x0, y1-y0))
x, y = np.linspace(x0, x1, length), np.linspace(y0, y1, length)

# Extract the values along the line
zi = z[x.astype(np.int), y.astype(np.int)]

#-- Plot...
fig, axes = plt.subplots(nrows=2)
axes[0].imshow(z)
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')

axes[1].plot(zi)

plt.show()

Solution 2

Probably the easiest way to do this is to use scipy.interpolate.interp2d():

# construct interpolation function
# (assuming your data is in the 2-d array "data")
x = numpy.arange(data.shape[1])
y = numpy.arange(data.shape[0])
f = scipy.interpolate.interp2d(x, y, data)

# extract values on line from r1, c1 to r2, c2
num_points = 100
xvalues = numpy.linspace(c1, c2, num_points)
yvalues = numpy.linspace(r1, r2, num_points)
zvalues = f(xvalues, yvalues)

Solution 3

I've been testing the above routines with galaxy images and think I found a small error. I think a transpose needs to be added to the otherwise great solution provided by Joe. Here is a slightly modified version of his code that reveals the error. If you run it without the transpose, you can see the profile doesn't match up; with the transpose it looks okay. This isn't apparent in Joe's solution since he uses a symmetric image.

import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt
import scipy.misc # ADDED THIS LINE

#-- Generate some data...
x, y = np.mgrid[-5:5:0.1, -5:5:0.1]
z = np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)
lena = scipy.misc.lena()  # ADDED THIS ASYMMETRIC IMAGE
z = lena[320:420,330:430] # ADDED THIS ASYMMETRIC IMAGE

#-- Extract the line...
# Make a line with "num" points...
x0, y0 = 5, 4.5 # These are in _pixel_ coordinates!!
x1, y1 = 60, 75
num = 500
x, y = np.linspace(x0, x1, num), np.linspace(y0, y1, num)

# Extract the values along the line, using cubic interpolation
zi = scipy.ndimage.map_coordinates(z, np.vstack((x,y))) # THIS DOESN'T WORK CORRECTLY
zi = scipy.ndimage.map_coordinates(np.transpose(z), np.vstack((x,y))) # THIS SEEMS TO WORK CORRECTLY

#-- Plot...
fig, axes = plt.subplots(nrows=2)
axes[0].imshow(z)
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')

axes[1].plot(zi)

plt.show()

Here's the version WITHOUT the transpose. Notice that only a small fraction on the left should be bright according to the image but the plot shows almost half of the plot as bright.

Here's the version WITH the transpose. In this image, the plot seems to match well with what you'd expect from the red line in the image.

Solution 4

For a canned solution look into scikit-image's measure.profile_line function.

It's built on top of scipy.ndimage.map_coordinates as in @Joe's answer and has some extra useful functionality baked in.

Solution 5

Combining this answer with the Event Handling example on MPL's documentation, here's the code to allow for GUI-based dragging to draw/update your slice, by dragging on the plot data (this is coded for pcolormesh plots):

import numpy as np 
import matplotlib.pyplot as plt  

# Handle mouse clicks on the plot:
class LineSlice:
    '''Allow user to drag a line on a pcolor/pcolormesh plot, and plot the Z values from that line on a separate axis.

    Example
    -------
    fig, (ax1, ax2) = plt.subplots( nrows=2 )    # one figure, two axes
    img = ax1.pcolormesh( x, y, Z )     # pcolormesh on the 1st axis
    lntr = LineSlice( img, ax2 )        # Connect the handler, plot LineSlice onto 2nd axis

    Arguments
    ---------
    img: the pcolormesh plot to extract data from and that the User's clicks will be recorded for.
    ax2: the axis on which to plot the data values from the dragged line.


    '''
    def __init__(self, img, ax):
        '''
        img: the pcolormesh instance to get data from/that user should click on
        ax: the axis to plot the line slice on
        '''
        self.img = img
        self.ax = ax
        self.data = img.get_array().reshape(img._meshWidth, img._meshHeight)

        # register the event handlers:
        self.cidclick = img.figure.canvas.mpl_connect('button_press_event', self)
        self.cidrelease = img.figure.canvas.mpl_connect('button_release_event', self)

        self.markers, self.arrow = None, None   # the lineslice indicators on the pcolormesh plot
        self.line = None    # the lineslice values plotted in a line
    #end __init__

    def __call__(self, event):
        '''Matplotlib will run this function whenever the user triggers an event on our figure'''
        if event.inaxes != self.img.axes: return     # exit if clicks weren't within the `img` axes
        if self.img.figure.canvas.manager.toolbar._active is not None: return   # exit if pyplot toolbar (zooming etc.) is active

        if event.name == 'button_press_event':
            self.p1 = (event.xdata, event.ydata)    # save 1st point
        elif event.name == 'button_release_event':
            self.p2 = (event.xdata, event.ydata)    # save 2nd point
            self.drawLineSlice()    # draw the Line Slice position & data
    #end __call__

    def drawLineSlice( self ):
        ''' Draw the region along which the Line Slice will be extracted, onto the original self.img pcolormesh plot.  Also update the self.axis plot to show the line slice data.'''
        '''Uses code from these hints:
        http://stackoverflow.com/questions/7878398/how-to-extract-an-arbitrary-line-of-values-from-a-numpy-array
        http://stackoverflow.com/questions/34840366/matplotlib-pcolor-get-array-returns-flattened-array-how-to-get-2d-data-ba
        '''

        x0,y0 = self.p1[0], self.p1[1]  # get user's selected coordinates
        x1,y1 = self.p2[0], self.p2[1]
        length = int( np.hypot(x1-x0, y1-y0) )
        x, y = np.linspace(x0, x1, length),   np.linspace(y0, y1, length)

        # Extract the values along the line with nearest-neighbor pixel value:
        # get temp. data from the pcolor plot
        zi = self.data[x.astype(np.int), y.astype(np.int)]
        # Extract the values along the line, using cubic interpolation:
        #import scipy.ndimage
        #zi = scipy.ndimage.map_coordinates(self.data, np.vstack((x,y)))

        # if plots exist, delete them:
        if self.markers != None:
            if isinstance(self.markers, list):
                self.markers[0].remove()
            else:
                self.markers.remove()
        if self.arrow != None:
            self.arrow.remove()

        # plot the endpoints
        self.markers = self.img.axes.plot([x0, x1], [y0, y1], 'wo')   
        # plot an arrow:
        self.arrow = self.img.axes.annotate("",
                    xy=(x0, y0),    # start point
                    xycoords='data',
                    xytext=(x1, y1),    # end point
                    textcoords='data',
                    arrowprops=dict(
                        arrowstyle="<-",
                        connectionstyle="arc3", 
                        color='white',
                        alpha=0.7,
                        linewidth=3
                        ),

                    )

        # plot the data along the line on provided `ax`:
        if self.line != None:
            self.line[0].remove()   # delete the plot
        self.line = self.ax.plot(zi)
    #end drawLineSlice()

#end class LineTrace


# load the data:
D = np.genfromtxt(DataFilePath, ...)
fig, ax1, ax2 = plt.subplots(nrows=2, ncols=1)

# plot the data
img = ax1.pcolormesh( np.arange( len(D[0,:]) ), np.arange(len(D[:,0])), D )

# register the event handler:
LnTr = LineSlice(img, ax2)    # args: the pcolor plot (img) & the axis to plot the values on (ax2)

This results in the following (after adding axis labels etc.), after dragging on the pcolor plot:

View more solutions

Comments

[Sven Marnach]

Nice answer. The only point I don't get is why the solution I proposed is slower (I didn't do timings, so I'm not even convinced it is).

[achennu]

Thanks for that fantastic answer, and +5 for the eyecandy. I have learnt several things (and new functions!) from this comprehensive answer. May the stack never overflow on thee. :)

[achennu]

@SvenMarnach Perhaps, it is not going to be particularly slower actually, given that both the methods are essentially running interpolation operations on the array. However, the nearest-neighbor approach comes closest to answering my question -- but now I see that interpolation is perhaps not a bad way to go. Thank you also for your response.

[Joe Kington]

@Sven - For what it's worth, I think your answer as-is doesn't do what you think... Your zvalues will be a 100x100 2D array, not a 100-element 1D array. That aside, though, I've usually found map_coordinates to be faster than FITPACK's BivariateSpline (which is what scipy.interpolate.interp2d uses) for the use case of interpolating at a few points from a large regular grid.

[Sven Marnach]

@Joe: Oh, yes, I certainly expected the interpolation function to follow the usual broadcasting rules. Thanks for the correction.

[Demis]

This only works for pcolormesh, due to lack of an API call that returns the original data array of the pcolormesh (in img). The use of img._meshWidth, img._meshHeight probably breaks use in other plots. See here: stackoverflow.com/questions/34840366/…

[gazzar]

I just ran across this too, and changed to zi = scipy.ndimage.map_coordinates(z, np.vstack((y,x)))

[Gilly]

Good answer with diagrams. One comment, you shouldn't have used a rotationally symmetrical example: I think (at least in the second NN example) the x and y need to be swapped.

[gradi3nt]

There is an error in the second and third examples. zi = z[x.astype(np.int), y.astype(np.int)] should be zi = z[y.astype(np.int), x.astype(np.int)] since x is the column index and y is the row index.

[K.-Michael Aye]

The error is also in the 1st example, as map_coordinates looking up values in the image needs y and x as coordinates, not x and y.

[Valeria]

@JoeKington would you know any efficient solution for sampling not only the central pixel along the line but also its, say, 5 neighbors from each side? Thanks.