A crash course on NumPy for images
A crash course on NumPy for images
scikit-image
的图像是NumPy数组。因此很大一部分图像操作与NumPy相关:
>>> from skimage import data
>>> camera = data.camera()
>>> type(camera)
<type 'numpy.ndarray'>
检索图像的几何图形和像素数量:
>>> camera.shape
(512, 512)
>>> camera.size
262144
检索有关灰度值的统计信息:
>>> camera.min(), camera.max()
(0, 255)
>>> camera.mean()
118.31400299072266
表示图像的NumPy数组可以是不同的浮点数字类型的整数。有关这些类型的更多信息以及scikit-image如何处理它们,请参阅图像数据类型及其含义。
NumPy索引
NumPy索引既可用于查看像素值,也可用于修改像素值:
>>> # Get the value of the pixel on the 10th row and 20th column
>>> camera[10, 20]
153
>>> # Set to black the pixel on the 3rd row and 10th column
>>> camera[3, 10] = 0
注意:在NumPy索引中,第一个维(camera.shape[0]
)对应于行,而第二个(camera.shape[1]
)对应于列,而camera[0, 0]
左上角的origin()。这匹配矩阵/线性代数符号,但与笛卡尔(x,y)坐标相反。有关更多详细信息,请参阅下面的坐标约定。
除了单个像素之外,还可以使用NumPy的不同索引可能性来访问/修改整组像素的值。
切片:
>>> # Set to black the ten first lines
>>> camera[:10] = 0
掩蔽(用布尔值掩码索引):
>>> mask = camera < 87
>>> # Set to "white" (255) pixels where mask is True
>>> camera[mask] = 255
花式索引(索引与索引)
>>> inds_r = np.arange(len(camera))
>>> inds_c = 4 * inds_r % len(camera)
>>> camera[inds_r, inds_c] = 0
尤其是,使用蒙版选择一组要执行进一步操作的像素非常有用。掩码可以是任何与图像形状相同的布尔阵列(或可以广播到图像形状的形状)。这对定义一个感兴趣的区域很有用,例如磁盘:
>>> nrows, ncols = camera.shape
>>> row, col = np.ogrid[:nrows, :ncols]
>>> cnt_row, cnt_col = nrows / 2, ncols / 2
>>> outer_disk_mask = ((row - cnt_row)**2 + (col - cnt_col)**2 >
... (nrows / 2)**2)
>>> camera[outer_disk_mask] = 0
布尔运算可用于定义更复杂的掩码:
>>> lower_half = row > cnt_row
>>> lower_half_disk = np.logical_and(lower_half, outer_disk_mask)
>>> camera = data.camera()
>>> camera[lower_half_disk] = 0
彩色图像
以上所有内容对于彩色图像也是如此:彩色图像是一个NumPy数组,并且为通道添加了一个附加的拖尾维度:
>>> cat = data.chelsea()
>>> type(cat)
<type 'numpy.ndarray'>
>>> cat.shape
(300, 451, 3)
这显示了cat
具有三个通道(红色,绿色和蓝色)的300x451像素图像。和以前一样,我们可以获取并设置像素值:
>>> cat[10, 20]
array([151, 129, 115], dtype=uint8)
>>> # set the pixel at row 50, column 60 to black
>>> cat[50, 60] = 0
>>> # set the pixel at row 50, column 61 to green
>>> cat[50, 61] = [0, 255, 0] # [red, green, blue]
我们也可以使用二维布尔模板来渲染二维彩色图像,就像我们上面的灰度图一样:
在2D彩色图像上使用2D蒙版
>>> from skimage import data
>>> cat = data.chelsea()
>>> reddish = cat[:, :, 0] > 160
>>> cat[reddish] = [0, 255, 0]
>>> plt.imshow(cat)
(源代码,png,pdf)
协调惯例
因为我们用numpy
数组表示图像,所以我们的坐标必须相应匹配。二维(2D)灰度图像(例如camera
上图)按行和列(缩写为(row, col)
或(r, c)
)进行索引,(0, 0)
最左边的元素位于左上角。在图书馆的各个部分,您还将看到rr
并cc
参考行列坐标列表。我们将它与(x, y)
通常表示笛卡尔坐标的标准区分开来,其中x
水平坐标,y
垂直坐标和原点在左下角。(例如,Matplotlib使用这个约定。)
在彩色(或多通道)图像的情况下,最后一个维度包含颜色信息并用channel
或表示ch
。
最后,对于诸如视频,磁共振成像(MRI)扫描或共聚焦显微镜之类的3D图像,我们称主要维度为plane
,缩写为pln
或p
。
这些约定总结如下:
图像类型 | 坐标 |
---|---|
2D灰度 | (行,列) |
2D多通道(例如RGB) | (row,col,ch) |
3D灰度 | (pln,row,col) |
3D多声道 | (pln,row,col,ch) |
scikit-image中的许多功能都直接在3D图像上运行:
>>> im3d = np.random.rand(100, 1000, 1000)
>>> from skimage import morphology
>>> from scipy import ndimage as ndi
>>> seeds = ndi.label(im3d < 0.1)[0]
>>> ws = morphology.watershed(im3d, seeds)
在许多情况下,第三个成像维度的分辨率低于其他两个。一些scikit-image函数提供了一个spacing
关键字参数来处理这些图像:
>>> from skimage import segmentation
>>> slics = segmentation.slic(im3d, spacing=[5, 1, 1], multichannel=False)
其他时候,处理必须按平面进行。当飞机是主要维度时,我们可以使用以下语法:
>>> from skimage import filters
>>> edges = np.zeros_like(im3d)
>>> for pln, image in enumerate(im3d):
... # iterate over the leading dimension (planes)
... edges[pln] = filters.sobel(image)
有关阵列顺序的说明
尽管轴的标记看起来是任意的,但它对操作速度可能有显着影响。这是因为现代处理器从不从内存中检索一个项目,而是从整个相邻项目中取出一个项目。(这称为预取)。因此,即使操作次数相同,在内存中彼此相邻的处理元素也会比以不同顺序处理元素更快。
>>> def in_order_multiply(arr, scalar):
... for plane in list(range(arr.shape[0])):
... arr[plane, :, :] *= scalar
...
>>> def out_of_order_multiply(arr, scalar):
... for plane in list(range(arr.shape[2])):
... arr[:, :, plane] *= scalar
...
>>> import time
>>> im3d = np.random.rand(100, 1024, 1024)
>>> t0 = time.time( x = in_order_multiply(im3d, 5 t1 = time.time()
>>> print("%.2f seconds" % (t1 - t0))
0.14 seconds
>>> im3d_t = np.transpose(im3d).copy() # place "planes" dimension at end
>>> im3d_t.shape
(1024, 1024, 100)
>>> s0 = time.time( x = out_of_order_multiply(im3d, 5 s1 = time.time()
>>> print("%.2f seconds" % (s1 - s0))
1.18 seconds
>>> print("Speedup: %.1fx" % ((s1 - s0) / (t1 - t0)))
Speedup: 8.6x
当您迭代的维度更大时,加速更加戏剧性。编写算法时值得考虑这个数据局部性
。特别是,除非另有说明,否则scikit-image会使用C连续数组,因此应该沿着计算最内层循环中的最后/最右边的维度进行迭代。
关于时间的说明
虽然scikit-image目前不提供函数来专门处理随时间变化的3D数据,但我们与numpy数组的兼容性使我们能够非常自然地使用形状的5D阵列(t,pln,row,col,ch ):
>>> for timepoint in image5d:
... # each timepoint is a 3D multichannel image
... do_something_with(timepoint)
然后我们可以补充上述表格如下:
图像类型 | 坐标 |
---|---|
2D彩色视频 | (t,row,col,ch) |
3D多声道视频 | (t,pln,row,col,ch) |