【问题标题】:How to extract an arbitrary line of values from a numpy array?如何从 numpy 数组中提取任意行值?
【发布时间】:2011-12-14 06:30:38
【问题描述】:

我有一个包含一些图像数据的 numpy 数组。我想绘制在图像上绘制的横断面的“轮廓”。最简单的情况是与图像边缘平行运行的配置文件,因此如果图像数组是imdat,那么在选定点(r,c) 的配置文件就是imdat[r](水平)或imdat[:,c](垂直) .

现在,我想将两个点 (r1,c1)(r2,c2) 作为输入,它们都位于 imdat 内。我想沿着连接这两点的线绘制值的轮廓。

沿着这条线从 numpy 数组中获取值的最佳方法是什么?更一般地,沿着路径/多边形?

我以前使用过切片和索引,但对于连续切片元素不在同一行或同一列中的情况,我似乎无法找到一个优雅的解决方案。感谢您的帮助。

【问题讨论】:

  • 哪一行?不能保证数组中的两个任意条目之间存在唯一的“线”。只有两个结尾条目位于同一行、同一列、同一对角线或反对角线时才会存在这样一条唯一的行。
  • 确实如此,因为“线”必须以不均匀的方式穿过像素,并且在不同的计算中可能会生成不同的线。但是,我主要对从起点 (r1,c1) 到 (r2,c2) 沿这个给定“方向”的图像值的趋势感兴趣。选择线路的特殊性对我的需求并不重要。

标签: python numpy indexing slice


【解决方案1】:

这是一种不使用 scipy 包的方法。它应该运行得更快并且易于理解。基本上,点 1 (pt1) 和点 2 (pt2) 之间的任何坐标对都可以转换为 x 和 y 像素整数,因此我们不需要任何插值。

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

def euclideanDistance(coord1,coord2):
    return np.sqrt((coord1[0]-coord2[0])**2+(coord1[1]-coord2[1])**2)

def getLinecut(image,X,Y,pt1,pt2):
    row_col_1, row_col_2 = getRowCol(pt1,X,Y), getRowCol(pt2,X,Y)
    row1,col1 = np.asarray(row_col_1).astype(float)
    row2,col2 = np.asarray(row_col_2).astype(float)
    dist = np.sqrt((pt1[0]-pt2[0])**2+(pt1[1]-pt2[1])**2)
    N = int(euclideanDistance(row_col_1,row_col_2))#int(np.sqrt((row1-row2)**2+(col1-col2)**2))
    rowList = [int(row1 + (row2-row1)/N*ind) for ind in range(N)]
    colList = [int(col1 + (col2-col1)/N*ind) for ind in range(N)]
    distList = [dist/N*ind for ind in range(N)]
    return distList,image[rowList,colList]#rowList,colList

def getRowCol(pt,X,Y):
    if X.min()<=pt[0]<=X.max() and Y.min()<=pt[1]<=Y.max():
        pass
    else:
        raise ValueError('The input center is not within the given scope.')
    center_coord_rowCol = (np.argmin(abs(Y-pt[1])),np.argmin(abs(X-pt[0])))
    return center_coord_rowCol

image = np.asarray(Image.open('./Picture1.png'))[:,:,1]
image_copy = image.copy().astype(float)

X = np.linspace(-27,27,np.shape(image)[1])#[::-1]
Y = np.linspace(-15,15,np.shape(image)[0])[::-1]

pt1, pt2 = (-12,-14), (20,13)
distList, linecut = getLinecut(image_copy,X,Y,pt1,pt2)
plt.plot(distList, linecut)

plt.figure()
plt.pcolormesh(X,Y,image_copy)
plt.plot([pt1[0],pt2[0]],[pt1[1],pt2[1]],color='red')
plt.gca().set_aspect(1)

使用的Picture1.png图: 有关更多详细信息,请参见此处: https://github.com/xuejianma/fastLinecut_radialLinecut

代码还有另一个功能:取几条角度均匀分布的线的平均值。

【讨论】:

    【解决方案2】:

    对于固定解决方案,请查看 scikit-imagemeasure.profile_line 函数。

    它建立在 scipy.ndimage.map_coordinates 之上,就像在 @Joeanswer 中一样,并且包含一些额外有用的功能。

    【讨论】:

      【解决方案3】:

      将此答案与Event Handling example on MPL's documentation 相结合,这是允许基于 GUI 的拖动以通过拖动绘图数据来绘制/更新切片的代码(这是为 pcolormesh 绘图编码的):

      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)
      

      这会导致以下结果(在添加轴标签等之后),在 pcolor 图上拖动后:

      【讨论】:

      • 这仅适用于pcolormesh,因为缺少返回pcolormesh(在img 中)的原始数据数组的API 调用。 img._meshWidth, img._meshHeight 的使用可能会破坏其他情节的使用。见这里:stackoverflow.com/questions/34840366/…
      【解决方案4】:

      我一直在用星​​系图像测试上述例程,并认为我发现了一个小错误。我认为需要将转置添加到乔提供的其他出色解决方案中。这是他的代码的一个稍微修改的版本,它揭示了错误。如果你在没有转置的情况下运行它,你会看到配置文件不匹配;使用转置它看起来还可以。这在 Joe 的解决方案中并不明显,因为他使用的是对称图像。

      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()
      

      这是没有转置的版本。请注意,根据图像,左侧只有一小部分应该是亮的,但该图显示几乎一半的图是亮的。

      这是带有转置的版本。在这张图片中,情节似乎与您对图片中红线的预期非常匹配。

      【讨论】:

      • 我也刚好遇到这个,改成zi = scipy.ndimage.map_coordinates(z, np.vstack((y,x)))
      【解决方案5】:

      @Sven 的答案是简单的方法,但对于大型数组来说效率相当低。如果您正在处理一个相对较小的数组,您不会注意到差异,如果您想要一个较大的配置文件(例如 > 50 MB),您可能需要尝试其他几种方法。但是,您需要在“像素”坐标中处理这些问题,因此会有额外的复杂性。

      还有两种更节省内存的方法。 1) 如果您需要双线性或三次插值,请使用scipy.ndimage.map_coordinates。 2)如果你只是想要最近邻采样,那么直接使用索引。

      以第一个为例:

      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()
      

      使用最近邻插值的等效项如下所示:

      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()
      

      但是,如果您使用最近邻,您可能只需要每个像素的样本,所以您可能会做更多类似的事情,而不是...

      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()
      

      【讨论】:

      • 不错的答案。我唯一不明白的是为什么我提出的解决方案较慢(我没有计时,所以我什至不相信它是)。
      • 感谢您的精彩回答,并感谢您的赏心悦目 +5。我从这个全面的答案中学到了几件事(和新功能!)。愿堆栈永远不会在你身上溢出。 :)
      • @SvenMarnach 也许,它实际上不会特别慢,因为这两种方法本质上都是在数组上运行插值操作。然而,最近邻方法最接近回答我的问题——但现在我看到插值可能不是一个坏方法。也感谢您的回复。
      • 图表的好答案。一条评论,您不应该使用旋转对称示例:我认为(至少在第二个 NN 示例中)x 和 y 需要交换。
      • 第二个和第三个例子有错误。 zi = z[x.astype(np.int), y.astype(np.int)] 应该是 zi = z[y.astype(np.int), x.astype(np.int)],因为 x 是列索引,y 是行索引。
      【解决方案6】:

      可能最简单的方法是使用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)
      

      【讨论】:

        猜你喜欢
        • 2017-01-19
        • 1970-01-01
        • 2023-02-02
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        相关资源
        最近更新 更多