ach*_*nnu 53 python indexing numpy slice
我有一个包含一些图像数据的numpy数组.我想绘制在图像上绘制的横断面的"轮廓".最简单的情况是平行于图像边缘运行的轮廓,因此如果图像阵列是imdat,则所选点处的轮廓(r,c)简单地imdat[r](水平)或imdat[:,c](垂直).
现在,我要带作为输入两个点(r1,c1)和(r2,c2),卧都里面imdat.我想在连接这两点的线上绘制值的轮廓.
沿着这样一条线从numpy数组中获取值的最佳方法是什么?更一般地说,沿路径/多边形?
我之前使用过切片和索引,但是对于连续切片元素不在同一行或列中的情况,我似乎无法找到优雅的解决方案.谢谢你的帮助.
Joe*_*ton 81
@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()
Run Code Online (Sandbox Code Playgroud)

使用最近邻插值的等价物看起来像这样:
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()
Run Code Online (Sandbox Code Playgroud)

但是,如果你使用的是最近邻,你可能只想在每个像素处采样,所以你可能会做更像这样的事情,而不是......
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()
Run Code Online (Sandbox Code Playgroud)

acr*_*der 18
我一直用galaxy图像测试上面的例程,并认为我发现了一个小错误.我认为需要将转置添加到Joe提供的其他优秀解决方案中.这是他的代码的略微修改版本,揭示了错误.如果在没有转置的情况下运行它,您可以看到配置文件不匹配; 转置它看起来没问题.这在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()
Run Code Online (Sandbox Code Playgroud)
这是没有转置的版本.请注意,根据图像,左边只有一小部分应该是明亮的,但是图表显示几乎一半的图是明亮的.

这是带转置的版本.在此图像中,图表似乎与您对图像中红线的期望相匹配.

Sve*_*ach 17
可能最简单的方法是使用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)
Run Code Online (Sandbox Code Playgroud)
Dav*_*man 12
对于罐装的解决方案考虑scikit-image的measure.profile_line功能.
它建立scipy.ndimage.map_coordinates在@Joe的答案之上,并且有一些额外的有用功能.