4. 关于图像的NumPy速成课程#

scikit-image 中的图像由 NumPy ndarrays 表示。因此，许多常见的操作可以通过用于操作数组的标准 NumPy 方法来实现:

>>> import skimage as ski
>>> camera = ski.data.camera()
>>> type(camera)
<type 'numpy.ndarray'>

备注

标记的类似数组的数据类型，如 pandas.DataFrame 或 xarray.DataArray，在 scikit-image 中并不原生支持。然而，存储在这些类型中的数据可以通过某些假设转换为 numpy.ndarray``（参见 ``pandas.DataFrame.to_numpy() 和 xarray.DataArray.data）。特别是，这些转换会忽略采样坐标（DataFrame.index、DataFrame.columns 或 DataArray.coords），这可能导致数据表示错误，例如，当原始数据点间距不规则时。

获取图像的几何形状和像素数量:

>>> camera.shape
(512, 512)
>>> camera.size
262144

获取图像强度值的统计信息:

>>> camera.min(), camera.max()
(0, 255)
>>> camera.mean()
118.31400299072266

表示图像的 NumPy 数组可以是不同的整数或浮点数类型。有关这些类型的更多信息以及 scikit-image 如何处理它们，请参阅图像数据类型及其含义。

4.1. NumPy 索引#

NumPy 索引既可以用于查看像素值，也可以用于修改它们:

>>> # Get the value of the pixel at the 10th row and 20th column
>>> camera[10, 20]
153
>>> # Set to black the pixel at the 3rd row and 10th column
>>> camera[3, 10] = 0

小心！在 NumPy 索引中，第一个维度 (camera.shape[0]) 对应于行，而第二个维度 (camera.shape[1]) 对应于列，原点 (camera[0, 0]) 位于左上角。这与矩阵/线性代数表示法相匹配，但与笛卡尔坐标系 (x, y) 相反。有关更多详细信息，请参见下面的坐标约定。

除了单个像素之外，还可以使用NumPy的不同索引功能来访问/修改整组像素的值。

切片:

>>> # Set the first ten lines to "black" (0)
>>> camera[:10] = 0

掩码（使用布尔掩码进行索引）:

>>> mask = camera < 87
>>> # Set to "white" (255) the pixels where mask is True
>>> camera[mask] = 255

花式索引（使用索引集进行索引）:

>>> import numpy as np
>>> 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

../_images/sphx_glr_plot_camera_numpy_001.png

NumPy 的布尔运算可以用来定义更复杂的掩码:

>>> lower_half = row > cnt_row
>>> lower_half_disk = np.logical_and(lower_half, outer_disk_mask)
>>> camera = data.camera()
>>> camera[lower_half_disk] = 0

4.2. 彩色图像#

对于彩色图像，上述所有内容仍然适用。彩色图像是一个具有额外尾随维度的 NumPy 数组，用于通道:

>>> cat = ski.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 (50th row, 60th column) to "black"
>>> cat[50, 60] = 0
>>> # set the pixel at (50th row, 61st column) to "green"
>>> cat[50, 61] = [0, 255, 0]  # [red, green, blue]

我们也可以对2D多通道图像使用2D布尔掩码，就像我们上面处理灰度图像那样：

(Source code, png, hires.png, pdf)

../_images/numpy_images-1.png — 在2D彩色图像上使用2D掩码#

包含在 skimage.data 中的示例彩色图像的通道沿着最后一个轴存储，尽管其他软件可能遵循不同的约定。scikit-image 库支持彩色图像的函数有一个 channel_axis 参数，可以用来指定数组的哪个轴对应于通道。

4.3. 坐标约定#

因为 scikit-image 使用 NumPy 数组表示图像，坐标约定必须匹配。二维（2D）灰度图像（如上面的 camera）通过行和列进行索引（简写为 (row, col) 或 (r, c)），最低元素 (0, 0) 位于左上角。在库的各个部分，你还会看到 rr 和 cc 表示行和列坐标的列表。我们将这种约定与 (x, y) 区分开来，后者通常表示标准笛卡尔坐标，其中 x 是水平坐标，y 是垂直坐标，原点位于左下角（例如，Matplotlib 轴使用这种约定）。

在多通道图像的情况下，任何维度（数组轴）都可以用于颜色通道，并由 channel 或 ch 表示。在 scikit-image 0.19 之前，此通道维度始终是最后一个，但在当前版本中，通道维度可以通过 channel_axis 参数指定。需要多通道数据的函数默认使用 channel_axis=-1。否则，函数默认使用 channel_axis=None，表示没有轴被假定为对应于通道。

最后，对于体积（3D）图像，如视频、磁共振成像（MRI）扫描、共聚焦显微镜等，我们将主要维度称为``plane``，缩写为``pln``或``p``。

这些约定总结如下：

*scikit-image 中的维度名称和顺序约定*#
图像类型	坐标
2D 灰度	（行，列）
2D 多通道（例如 RGB）	（行，列，字符）
3D 灰度	(行, 列, 列)
3D 多通道	(pln, 行, 列, 字符)

注意 ch 的位置由 channel_axis 参数控制。

scikit-image 中的许多函数可以直接对 3D 图像进行操作:

>>> import numpy as np
>>> import scipy as sp
>>> import skimage as ski
>>> rng = np.random.default_rng()
>>> im3d = rng.random((100, 1000, 1000))
>>> seeds = sp.ndimage.label(im3d < 0.1)[0]
>>> ws = ski.segmentation.watershed(im3d, seeds)

然而，在许多情况下，第三空间维度的分辨率低于其他两个维度。一些 scikit-image 函数提供了一个 spacing 关键字参数来帮助处理这种类型的数据:

>>> slics = ski.segmentation.slic(im3d, spacing=[5, 1, 1], channel_axis=None)

其他时候，处理必须逐层进行。当层沿着前导维度堆叠时（符合我们的约定），可以使用以下语法：：

>>> edges = np.empty_like(im3d)
>>> for pln, image in enumerate(im3d):
...     # Iterate over the leading dimension
...     edges[pln] = ski.filters.sobel(image)

4.4. 数组维度顺序的注意事项#

尽管轴的标签可能看起来是任意的，但它对操作速度有显著影响。这是因为现代处理器从不只从内存中检索一个项目，而是检索一整块相邻的项目（这种操作称为预取）。因此，处理内存中相邻的元素比处理分散的元素更快，即使操作的数量相同:

>>> 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
>>> rng = np.random.default_rng()
>>> im3d = rng.random((100, 1024, 1024))
>>> t0 = time.time(); x = in_order_multiply(im3d, 5); t1 = time.time()
>>> print("%.2f seconds" % (t1 - t0))  
0.14 seconds
>>> 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连续数组。在使用嵌套循环时，数组的最后一个/最右边的维度应该在计算的最内层循环中。在上面的例子中，*= numpy 运算符遍历了所有剩余的维度。

4.5. 关于时间维度的说明#

尽管 scikit-image 目前没有提供专门处理随时间变化的3D数据的函数，但它与NumPy数组的兼容性使我们能够非常自然地处理形状为 (t, pln, row, col, ch) 的5D数组:

>>> for timepoint in image5d:  
...     # Each timepoint is a 3D multichannel image
...     do_something_with(timepoint)

然后我们可以如下补充上述表格：

*scikit-image 中维度名称和顺序的附录*#
图像类型	坐标
2D 彩色视频	(t, 行, 列, 字符)
3D 彩色视频	(t, pln, 行, 列, ch)