执行摘要 我们将探讨如何使用 dask-image 库创建一个基本的图像分割流程。
目录 只看代码
图像分割流程
自定义函数
扩展计算规模
附加内容：在 GPU 上使用数组
如何参与其中
这篇博文的内容最初是作为 2020 年的一次会议演讲 出现的。
只看代码 如果您想自己运行此示例，需要从 Broad Bioimage Benchmark Collection 下载示例数据：https://bbbc.broadinstitute.org/BBBC039
并安装这些依赖项
pip install dask-image>=0.4.0 tifffile
这是我们的完整流程
import numpy as np
from dask_image.imread import imread
from dask_image import ndfilters , ndmorph , ndmeasure
images = imread ( 'data/BBBC039/images/*.tif' )
smoothed = ndfilters . gaussian_filter ( images , sigma = [ 0 , 1 , 1 ])
thresh = ndfilters . threshold_local ( smoothed , blocksize = images . chunksize )
threshold_images = smoothed > thresh
structuring_element = np . array ([[[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ]], [[ 0 , 1 , 0 ], [ 1 , 1 , 1 ], [ 0 , 1 , 0 ]], [[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ]]])
binary_images = ndmorph . binary_closing ( threshold_image , structure = structuring_element )
label_images , num_features = ndmeasure . label ( binary_image )
index = np . arange ( num_features )
area = ndmeasure . area ( images , label_images , index )
mean_intensity = ndmeasure . mean ( images , label_images , index )
您可以继续阅读此图像分割流程的分步讲解，或直接跳至自定义函数 、扩展计算规模 或GPU 加速 部分。
图像分割流程 我们基本的图像分割流程包含五个步骤
读取数据 图像滤波 分割对象 形态学操作 测量对象 设置您的 Python 环境 在开始之前，我们需要设置 Python 虚拟环境。
至少您需要
pip install dask-image>=0.4.0 tifffile matplotlib
或者，您也可以安装 napari 图像查看器来可视化图像分割结果。
pip install "napari[all]"
要在 IPython 或 jupyter 中使用 napari，请在调用 napari 之前在单元格中运行 %gui qt magic 命令。更多详情请参阅 napari 文档 。
下载示例数据 我们将使用公开可用的图像数据集 BBBC039 (Caicedo 等人, 2018)，该数据集来自 Broad Bioimage Benchmark Collection (Ljosa 等人, Nature Methods, 2012 )。您可以在此处下载数据集：https://bbbc.broadinstitute.org/BBBC039
这些是荧光显微镜图像，我们可以看到单个细胞中的细胞核。
步骤 1：读取数据 我们图像分割流程的第一步是读取图像数据。我们可以使用 dask-image 的 imread 函数 来完成。
我们传入包含示例数据中所有 *.tif 图像的文件夹路径。
from dask_image.imread import imread
images = imread ( 'data/BBBC039/images/*.tif' )
默认情况下，磁盘上的每个 .tif 文件都成为了我们的 Dask 数组中的一个块。
步骤 2：图像滤波 通过少量模糊对图像进行去噪可以改善后续的分割效果。这是许多图像分割流程中常见的第一个步骤。我们可以使用 dask-image 的 gaussian_filter
from dask_image import ndfilters
smoothed = ndfilters . gaussian_filter ( images , sigma = [ 0 , 1 , 1 ])
步骤 3：分割对象 接下来，我们想将图像中的对象与背景分开。有很多不同的方法可以做到这一点。由于我们有荧光显微镜图像，我们将使用阈值分割方法。
绝对阈值 我们可以设置一个绝对阈值，低于此阈值强度的像素将被视为背景的一部分。
absolute_threshold = smoothed > 160
让我们用 napari 图像查看器来看一下这些图像。首先我们需要使用 %gui qt magic 命令
现在我们可以使用 napari 查看图像了
import napari
viewer = napari . Viewer ()
viewer . add_image ( absolute_threshold )
viewer . add_image ( images , contrast_limits = [ 0 , 2000 ])
但这有一个问题。
当我们查看不同图像帧的结果时，很明显无法使用一个“一刀切”的绝对阈值。数据集中的一些图像背景很亮，而另一些图像的荧光细胞核亮度较低。我们必须尝试一种不同的阈值分割方法。
局部阈值 我们可以使用局部阈值分割方法来改进分割效果。
如果我们为每个图像帧独立计算一个阈值，那么就可以避免帧间背景强度波动引起的问题。
thresh = ndfilters . threshold_local ( smoothed , images . chunksize )
threshold_images = smoothed > thresh
# Let's take a look at the images with napari
viewer . add_image ( threshold_images )
这里的结果看起来好多了，这是将细胞核与背景更清晰地分离，并且对所有图像帧都适用。
步骤 4：形态学操作 现在我们从阈值得到了二值掩膜，可以通过一些形态学操作进行清理。
形态学操作是我们对二值图像中的结构形状进行的更改。我们将在此简要介绍一些基本概念，但更详细的参考资料可以查阅 这篇优秀的 OpenCV 文档页面 。
腐蚀（Erosion） 是一种操作，其中二值图像中结构的边缘被侵蚀或“吃掉”。
图片来源：OpenCV 文档
膨胀（Dilation） 与腐蚀相反。膨胀会扩大二值图像中结构的边缘。
图片来源：OpenCV 文档
我们可以以不同的方式组合形态学操作以获得有用的效果。
形态学开运算（morphological opening） 是一种先腐蚀后膨胀的操作。
图片来源：OpenCV 文档
在上面的示例图像中，我们可以看到左侧背景有噪声和斑点。如果用于形态学操作的结构元素大于噪声斑点的大小，它们将在第一个腐蚀步骤中完全消失。然后进行第二个膨胀步骤时，背景中的噪声已经没有任何剩余可以膨胀了。这样我们就去除了背景噪声，同时我们感兴趣的主要结构（在这个例子中是 J 形）几乎完美地恢复了。
让我们利用这个形态学开运算技巧来清理我们的分割流程中的二值图像。
from dask_image import ndmorph
import numpy as np
structuring_element = np . array ([
[[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ]],
[[ 0 , 1 , 0 ], [ 1 , 1 , 1 ], [ 0 , 1 , 0 ]],
[[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ]]])
binary_images = ndmorph . binary_opening ( threshold_images , structure = structuring_element )
您会注意到这里我们需要对结构元素稍微小心一些。我们所有的图像帧都组合在一个 Dask 数组中，但我们只想独立地将形态学操作应用于每个帧。为此，我们将默认的二维结构元素夹在两层零之间。这意味着相邻的图像帧对结果没有影响。
# Default 2D structuring element
[[ 0 , 1 , 0 ],
[ 1 , 1 , 1 ],
[ 0 , 1 , 0 ]]
步骤 5：测量对象 任何图像处理流程的最后一步都是进行某种测量。我们将把二值掩膜转换为标签图像，然后测量这些对象的强度和大小。
为了让本教程的计算时间快速友好，我们将只测量数据的一个小 subset。让我们测量前三个图像帧中的所有对象（大约 300 个细胞核）。
from dask_image import ndmeasure
# Create labelled mask
label_images , num_features = ndmeasure . label ( binary_images [: 3 ], structuring_element )
index = np . arange ( num_features - 1 ) + 1 # [1, 2, 3, ...num_features]
这是根据我们的掩膜生成的标签图像的截图。
>>> print ( "Number of nuclei:" , num_features . compute ())
Number of nuclei : 271
测量图像中的对象 dask-image 的 ndmeasure 子包 包含了许多不同的测量函数。在这个例子中，我们将选择测量
每个对象的像素面积，以及 每个对象的平均强度。 area = ndmeasure . area ( images [: 3 ], label_images , index )
mean_intensity = ndmeasure . mean ( images [: 3 ], label_images , index )
运行计算并绘制结果 import matplotlib.pyplot as plt
plt . scatter ( area , mean_intensity , alpha = 0.5 )
plt . gca (). update ( dict ( title = "Area vs mean intensity" , xlabel = 'Area (pixels)' , ylabel = 'Mean intensity' ))
plt . show ()
自定义函数 如果您想做一些 dask-image API 中不包含的事情怎么办？我们可以使用几种选项来编写自定义函数。
Dask 的 map_overlap 和 map_blocks Dask 数组的 map_overlap map_blocks dask-image 中大多数函数的基础。您也可以自己使用它们。它们会将一个函数应用于 Dask 数组中的每个块。
import dask.array as da
def my_custom_function ( args ):
# ... does something really cool
result = da . map_overlap ( my_custom_function , my_dask_array , args )
您可以在此处阅读更多关于重叠计算的信息 。
Dask 的 delayed 如果您想对计算进行更灵活和细粒度的控制，那么可以使用 Dask delayed 。您可以在此处开始学习 Dask delayed 教程 。
scikit-image 的 apply_parallel 函数 如果您是经常在 Python 中进行图像处理的人，那么您可能已经在使用的工具之一是 scikit-image 。我想提请您注意 scikit-image 中提供的 apply_parallel map-overlap ，并且非常有用。
它不仅在大数据时有用，而且在数据可以放入内存但您想应用于数据的计算是内存密集型的情况下也很有用。这可能会导致您超出可用 RAM，而 apply_parallel
扩展计算规模 当您想从笔记本电脑扩展到超级计算集群时，可以使用 dask-distributed 来处理。
from dask.distributed import Client
# Setup a local cluster
# By default this sets up 1 worker per core
client = Client ()
client . cluster
请参阅此处的文档 ，以根据您的系统进行设置。
附加内容：在 GPU 上使用数组 我们最近一直在为 dask-image 添加 GPU 支持。
我们能够使用一个名为 CuPy 的库来添加 GPU 支持。CuPy 是一个具有类似 numpy API 的数组库，由 NVIDIA CUDA 加速。我们可以使用包含 cupy 块的 Dask 数组，而不是包含 numpy 块的 Dask 数组。这篇博文 解释了 GPU 加速的优势，并提供了 CPU、单 GPU 和多 GPU 计算的一些基准测试结果。
dask-image 中可用的 GPU 支持 从 dask-image 0.6.0 版本开始，以下六个子包中的四个支持 GPU 数组：
imread ndfilters ndinterp ndmorph 尚不支持 GPU 的 dask-image 子包有：
我们希望将来为这些子包添加 GPU 支持。
一个例子 这是使用 Dask 在 CPU 上进行图像卷积的示例
# CPU example
import numpy as np
import dask.array as da
from dask_image.ndfilters import convolve
s = ( 10 , 10 )
a = da . from_array ( np . arange ( int ( np . prod ( s ))). reshape ( s ), chunks = 5 )
w = np . ones ( a . ndim * ( 3 ,), dtype = np . float32 )
result = convolve ( a , w )
result . compute ()
这是使用 Dask 在 GPU 上进行图像卷积的相同示例。唯一需要更改的是输入数组的类型。
# Same example moved to the GPU
import cupy # <- import cupy instead of numpy (version >=7.7.0)
import dask.array as da
from dask_image.ndfilters import convolve
s = ( 10 , 10 )
a = da . from_array ( cupy . arange ( int ( cupy . prod ( cupy . array ( s )))). reshape ( s ), chunks = 5 ) # <- cupy dask array
w = cupy . ones ( a . ndim * ( 3 ,)) # <- cupy array
result = convolve ( a , w )
result . compute ()
您不能在同一计算中混合使用 CPU 数组和 GPU 数组。这就是为什么在上面第二个示例中，权重 w 必须是 cupy 数组的原因。
此外，您可以在 CPU 和 GPU 之间传输数据。因此，在 GPU 加速大于数据传输成本的情况下，这样做可能很有用。
将图像读取到 GPU 上 通常，我们希望通过从磁盘存储的图像中读取数据来开始图像处理。我们可以使用 imread arraytype=cupy 关键字参数来完成此操作。
from dask_image.imread import imread
images = imread ( 'data/BBBC039/images/*.tif' )
images_on_gpu = imread ( 'data/BBBC039/images/*.tif' , arraytype = "cupy" )
如何参与其中 使用 Dask 创建和分享您自己的分割或图像处理流程（加入当前关于分割的讨论 或在此提出新的博文主题 ）
贡献力量为 dask-image 添加 GPU 支持：https://github.com/dask/dask-image/issues/133
