Avizo用户使用手册-6
[TOC]
Chapter 6
6 教程:高级图像处理、分割与分析
以下教程需要Avizo许可证(Avizo Lite Edition许可证是不够的)。
Getting Started with Advanced Image Processing and Quantitative Analysis高级图像处理与定量分析入门——使用Avizo进行图像处理和分析的基础- 示例:
Measuring a Catalyst测量催化剂——掩膜与距离图Separating, Measuring and Reconstructing分离、测量与重建——分水岭分离Distribution of Pore Diameters in Foam泡沫中孔径的分布——自定义测量Average Thickness of Material in Foam泡沫中材料的平均厚度——分离厚度
Watershed Segmentation分水岭分割——高级分割More about Image Filtering关于图像滤波的更多内容——降低图像噪声或伪影,或增强感兴趣的特征More about label measures关于标签测量的更多内容——管理多组测量Cavity Analysis Tutorial空腔分析教程——使用环境光遮蔽进行空腔分析
6.1 高级图像处理与定量分析入门
在这个逐步教程中,你将学习使用Avizo进行图像处理和分析的基础知识。下面所用的示例可以很容易地扩展到其他应用中,并且遵循图像分析的典型工作流:
- 图像增强
- 特征提取
- 数据测量与分析
本节包含以下部分:
Processing grayscale images处理灰度图像3D versus 2D stack interpretation3D与2D堆栈解释Binarization of grayscale images灰度图像的二值化Separation分离Analysis and measures分析与测量Interactive selection交互式选择Measure filters测量滤波器Sieves筛分Label images标签图像Processing images in memory (in-core) and on disk (out-of-core)在内存中(核内)和在磁盘上(核外)处理图像Scripting脚本编写
要学习本教程,你应该熟悉Avizo的基本概念。特别是,你应该能够加载文件、与3D查看器交互,以及将模块连接到数据模块。所有这些问题都在Avizo第2章——入门中讨论过。为了日常使用Avizo,熟悉Avizo的图像滤波和分割也会对你有所帮助(参见第3章——图像与体积的可视化和处理)。
6.1.1 处理图像
首先,你需要:
- 在Avizo中从
AVIZO_ROOT目录加载3D显微断层扫描体积data/images/foam/foam.am。该数据对象会出现在Project View中。 - 如果
Auto-Display被禁用,请附加一个Ortho Slice模块以可视化数据。为此,右键单击绿色数据图标。在对象弹出菜单中,选择名为Display的类别条目,然后双击Ortho Slice条目(或单击Create)。你也可以使用搜索字段并键入Ortho Slice的首字母,然后在列表中选择该模块,它就会自动在Project View中被创建。 - 再附加两个
Ortho Slice模块。在第二个Ortho Slice的Properties面板中选择xz方向,在第三个Ortho Slice的Properties面板中选择yz方向。所得到的显示见图6.1。
图 6.1:初始数据。
加载项目data/tutorials/image-processing-advanced/GettingStartedBasics-1-LoadData.hx即可完成上述步骤。
为了通过图像分析获得最佳结果,通常需要改善图像质量。下一步将说明如何在Avizo中使用图像滤波器来处理图像——图像滤波器通常用于平滑或降噪。你可以在教程第6.7节——关于图像滤波的更多内容中了解更多关于图像滤波器的信息。
- 将一个
Median Filter模块附加到数据。(使用搜索字段并键入Median Filter的首字母,然后在列表中选择该模块。) - 在
Median Filter的Interpretation端口中选择3D。 - 按下
Apply按钮。
计算完成后,结果被存储在一个新的图像对象foam.filtered中(见图6.2)。
- 将这三个
Ortho Slice附加到所得到的图像。为此,在Properties面板中更改每个切片的Data。所得到的显示见图6.3。
加载项目data/tutorials/image-processing-advanced/GettingStartedBasics-2-ImageProcessing.hx即可完成上述步骤。
图 6.2:Median Filter网络。
图 6.3:使用中值滤波器去除噪声之后。
6.1.2 解释为3D图像或2D图像堆栈
有时,把图像处理算法的输入数据解释为一个3D体积,或者解释为一系列2D平面,会很有用。例如,许多图像滤波器和图像处理算法既可以使用2D内核在每个XY切片上执行,也可以使用3D内核在整个体积上执行。在某些情况下,可能更倾向于使用2D算法——出于性能考虑,或者为了根据数据和期望的结果获得更合适的效果。
在许多Avizo模块中,interpretation解释端口显示了当前模块的状态(即XY planes或3D)。如果该端口的状态是XY planes,意味着该模块将在每个XY切片上执行。如果该端口的状态是3D,意味着该模块将一次性在整个三维图像上执行。
在某些情况下,解释端口无法被更改(它是灰色的),例如,当处理只能在XY平面上应用时。
6.1.3 获取更多帮助
按下模块Properties面板中的问号按钮可显示该模块的帮助。该帮助可能是上下文相关的,取决于解释模式或模块中所选的处理类型。
6.1.4 二值化
二值化意味着把灰度图像转换为二值图像(即只有内部和外部两种材料的标签图像)。当灰度图像中的相关信息对应于某个特定的灰度值区间时,就会使用阈值二值化。阈值处理是一种简单的segmentation分割方法——Avizo中也提供了更复杂的自动、半自动或手动分割工具。阈值二值化可以用Interactive Thresholding模块完成,该模块会提示你在可视化反馈下设置阈值级别。
- 将一个
Interactive Thresholding模块附加到已滤波的数据。
你可以交互式地修改阈值,并立即获得2D或3D的视觉反馈。被选中的像素在显示的图像中呈蓝色。
- 在该模块的
Properties面板中,你可以把Intensity range端口设置为例如0-30的范围。 - 勾选和取消勾选
Preview Type端口中的3D选项,以获得仅2D或3D的预览(见图6.4)。 - 按下
Apply按钮以启动该模块。
名为foam.thresholded的输出二值图像会在Project View中生成(见图6.5)。
图 6.4:Interactive Thresholding的2D预览。
在输出的二值图像中,所有初始灰度值位于两个边界之间的像素都被设置为1,其他所有像素都被设置为0。
Interactive Thresholding模块创建了一个二值图像。对于二值图像,Avizo以蓝色显示强度为1的体素。如果你把一个Ortho Slice附加到所得到的图像,默认会选择一个合适的颜色图。
图 6.5:Interactive Thresholding网络。
- 通过切换
Project View中模块图标里的橙色可见性按钮,或Properties面板中模块名称旁边的按钮,隐藏Interactive Thresholding的预览。 - 把第一个
Ortho Slice连接到阈值处理后的数据。所得到的显示见图6.6。
图 6.6:用Ortho Slice显示的二值图像。
加载项目data/tutorials/image-processing-advanced/GettingStartedBasics-3-Binarization.hx即可完成上述步骤。
6.1.5 关于二值图像的更多内容
在二值图像中,所有满足某组条件的像素(这里的条件是像素强度位于Interactive Thresholding所设的两个边界之内)都被设置为值1(感兴趣的像素),所有其他像素都被设置为0(背景)。在Avizo中,二值图像没有特定的类型——Avizo的二值图像就是只有一个标签的标签图像(每体素8/16/32位,值为1;外部值为0)。不过,Avizo的某些模块可能明确要求二值图像作为输入数据。
6.1.6 关于二值化的更多提示
(本节是可选的,完成后续教程并不需要阅读本节。)
Avizo中提供了大量工具,用于有效地进行数据二值化和图像分割。在许多情况下,该过程可以自动化,可能需要组合多个步骤,有时需要用户输入。在某些情况下,进行半自动或手动分割可能是必要的或更容易的:特别是,Segmentation Editor就是为此目的而设计的(注意Segmentation Editor仅支持8位标签图像)。另外请记住,改善图像采集可能比分析质量差的图像要容易得多。
以下是一些可以促进自动二值化的Avizo模块:
Auto Thresholding自动计算一个阈值。你可以选择最适合你数据的准则,通常是factorisation因子化。Interactive Top-Hat是一个强大的工具,用于分割背景不均匀的区域——当简单的阈值处理无法在不引入不需要的噪声的情况下捕获所需特征时。顶帽变换可以看作是一种”局部阈值处理”。顶帽的结果常常使用逻辑运算OR Image与阈值处理的结果相结合。Hysteresis Thresholding用于在低阈值和高阈值之间实现中间的二值化,这两个阈值分别定义了一个安全保留区域和一个被拒绝的区域。例如,你可以使用Interactive Thresholding界面交互式地选择阈值,然后再使用Hysteresis Thresholding。另请参阅Canny Edge Detector模块。- 许多图像滤波器(如梯度或拉普拉斯算子)可用于辅助二值化,例如用于
Edge Detection边缘检测。 - 基于测量来过滤区域(如本教程后面所示)也可以成为一种强大的分割技术。
二值化之后,可能需要分离某些对象,如下一节所示。
6.1.7 分离
在示例数据集中,泡沫中的一些孔隙看起来是相互接触的,但理想情况下,为了进行正确的分析,应当把它们分离开。当采集数据过于粗糙或噪声过大时,阈值处理无法避免这类输出——因为所考虑对象的灰度值在整个体积中不够均匀,或者因为分辨率太低而无法区分某些对象的边界。你可以使用Separate Objects模块来分离相连的颗粒。
- 将一个
Separate Objects模块附加到阈值处理后的数据。 - 把
Marker Extent端口的值设置为1,而不是默认值4。这是一个对比因子,控制着标记待分离对象的种子的大小。增大该值可能会合并某些标记,从而减少被分离出的对象数量。 - 然后按下
Apply按钮。Project View中会生成一个foam.separate数据。 - 把第一个
Ortho Slice附加到新数据。所得到的显示见图6.7。
Separate Objects模块的原理是在距离图上计算分水岭线。Separate Objects模块是分水岭、距离图和H-Maxima的高层次组合。它可以作为一个简单直接的分离工具使用,在许多情况下都能令人满意。不过,你可能会注意到某些分离可能缺失或不理想,特别是对于非凸形状(也要考虑3D的情况)。有关更多细节和高级分离,请参阅Example 2: Separating, Measuring, and Reconstructing Individual Objects - Pores in Foam。
图 6.7:颗粒现已被分离。
加载项目data/tutorials/image-processing-advanced/GettingStartedBasics-4-Separation.hx即可完成上述步骤。
6.1.8 分析
然后,你可以使用分析模块来分别获取每个被分离颗粒的体积、表面积、平均值、体素数量等。对图像堆栈的这种分析是通过使用Label Analysis模块来完成的,该模块用于提取统计和数值信息,包括对对象的测量。
- 将一个
Label Analysis模块附加到已分离的数据。 - 在该模块的专用端口中,把
foam.am设置为Intensity Image。 - 按下
Apply按钮。
Project View中会创建一个新的label图像数据对象foam.label,并显示Tables面板,其中以电子表格样式的表格展示结果:分析结果foam.Label-Analysis也会在Project View中被创建(见图6.8和图6.9)。
工具栏提供了以下功能:
- 复制表格的部分内容
- 以多种格式导出电子表格
图 6.8:分析网络。
图 6.9:分析电子表格。
- 按升序或降序对列进行排序
- 绘制与某个测量相对应的直方图(见下文)
- 执行标签查找
label seek(见下面第6.1.9节——交互式选择)
basic基本测量(在该模块的Measures端口中选择)会显示在Tables面板中,如下所示(Volume3d、Area3D等)。
- 在下方的电子表格中选择
Volume3d列。 - 按下工具栏中的直方图按钮(见图6.10)。
一个窗口会打开,显示如图6.11所示的Volume3d直方图。
在Measures端口中,basic是一组包含最常用测量的测量组:
图 6.10:分析面板工具栏:直方图按钮。
图 6.11:Volume3d直方图。
Volume3d、Area3d、BaryCenterX、BaryCenterY、BaryCenterZ和Mean。也可以定义新的测量组,这些测量组既由预定义的测量组成,也可以包含用户自定义的测量。要了解更多相关内容,请参阅第6.4节——进一步的图像分析中的专门教程。
注意:默认情况下,如果Avizo的单位管理未启用,结果将以指定体素大小时所用的相同单位给出。你可以启用Avizo的单位管理,以带单位的方式显示数据和测量结果。关于如何使用Avizo中的单位管理的所有细节,请参阅第10.2.9节——Avizo中的单位。
加载项目data/tutorials/image-processing-advanced/GettingStartedBasics-5-Analysis.hx即可完成上述步骤。
6.1.9 交互式选择
Avizo允许你把3D查看器中的图像与分析面板中相应的行关联起来,以便定位具有相应测量值的单个对象。
图 6.12:分析面板工具栏:标签查找按钮。
- 单击分析面板工具栏中的标签查找按钮(见图6.12),或按下键盘上的
L键。一个新的Ortho Slice会自动附加到已分离的图像并显示在3D查看器中,同时还会显示一个点拖动器(point dragger)(见图6.13)。 - 选择分析面板下方表格中的一个单元格。点拖动器将移动到3D视图中相应的对象位置。如果显示了某个分析测量的直方图,则会出现一条垂直线,显示所选行的位置和数值,如图6.14所示。
- 在查看器窗口中,你可以使用矩形手柄移动拖动器,然后在释放按钮时,分析表格会高亮显示相应的对象行。为了移动拖动器,你必须把查看器设置为交互模式(按
ESC键)。然后把鼠标移到拖动器的某个十字准线上并按下鼠标左键。被拾取的十字准线颜色会发生变化。拖动器的移动被限制在相应的平面内。 - 你也可以用鼠标中键单击3D查看器中所显示场景里的某个可拾取对象,例如显示切片上的某个特定孔隙:拖动器将移动到被拾取的点,相应的电子表格行也会被高亮显示。
- 单击标签查找按钮或按下
L键可退出标签查找模式。
图 6.13:3D查看器中对应于第92行的点拖动器(棕色交叉线)。
6.1.10 基于测量的过滤
你可以对3D查看器中显示的颗粒进行过滤。例如,你可以决定只可视化Volume3d属于某个指定范围的颗粒。
- 将一个
Analysis Filter附加到foam.Label-Analysis。 - 把
Image端口连接到foam.label。 - 在
Filter端口中输入Volume3d >= 30000来创建一个新的过滤器。要在公式中插入Volume3d,你可以直接键入它,或者在公式字段下方显示的列表中双击它。 - 按下
Apply按钮。
图 6.14:直方图中对应于第92行的线。
这会创建一个包含较少对象的新分析。
- 你可以通过把
Ortho Slice连接到新的标签图像foam.label-filtering来验证创建的对象变少了,如图6.16所示。
图 6.15:Analysis Filter工作流。
提示:由测量驱动的过滤可以成为数据分割的强大工具。它允许你基于例如尺寸、形状因子、方向或多个准则的组合来选择或剔除区域。
加载项目data/tutorials/image-processing-advanced/GettingStartedBasics-6-FilterAnalysis.hx即可完成上述步骤。
图 6.16:过滤结果。
6.1.11 使用筛分对测量进行分类
你可以定义一组数值范围,然后用它来显示使用该分布的直方图,或者创建一个显示该分类的新标签图像。
- 将一个
Sieve Analysis模块附加到foam.Label-Analysis。 - 把
Data端口连接到foam.label。 - 通过更改
Number of values来添加一个数值,设为4。 - 通过编辑建议值或移动直方图中相应的标记来修改这些值。你也可以按下
Detect按钮以获得规则的间隔。 - 按下
Apply按钮。一个新的标签图像foam.Sieved已被创建。 - 使用
Volume Rendering模块显示该标签图像,如图6.18所示。
加载项目data/tutorials/image-processing-advanced/GettingStartedBasics-7-SieveAnalysis.hx即可完成上述步骤。
6.1.12 标签图像
- 隐藏
Volume Rendering - 把第一个
Ortho Slice附加到标签图像foam.label。
在与结果电子表格一起创建的label图像中,每个颗粒都已被识别并分配了一个唯一的索引。在本例中,该标签数据被存储为16位标签图像。此类图像默认使用循环颜色图显示,以便相邻的颗粒更有可能以不同的颜色呈现,如图6.19所示。
注意:来自Segmentation Editor和Multi-Thresholding模块的Avizo标签图像是每体素8位的标签图像。
图 6.17:Sieve Analysis工作流。
图 6.18:用Volume Rendering显示的筛分分析结果。
6.1.13 在磁盘上处理数据
在某些情况下,你可能会决定使用Visilog(.im6)数据格式,以便在打开数据文件时能够”把数据保留在磁盘上”,而不是把数据完全加载到内存中。这样,某些Avizo模块就可以例如逐切片地加载和处理图像数据(核外out-of-core),从而避免把完整数据加载到内存中(核内in-core)。这使得处理远大于系统可用内存的数据成为可能,代价是处理时间。
默认情况下,模块输出将以与模块输入相同的方式创建。因此,为了完全在磁盘上处理数据(输入和输出),你需要在打开输入数据时选择Stay on disk。
除了Visilog格式的文件之外,你还可以对其他数据格式(例如.lda)使用Stay on disk。你也可以把未压缩的Avizo文件、Raw文件作为Large Disk Data加载。
图 6.19:每个颗粒都已被识别并分配了一个唯一索引。
6.1.14 脚本编写
可以把一个完整的处理序列放入脚本中,以便为例行任务自动化分析过程。关于使用Avizo进行脚本编写的详细信息,请参阅第6.3节——示例3:分离、测量与重建。
6.1.15 结论
本教程向你介绍了:
- Avizo界面模块和帮助菜单,
- 灰度图像处理模块,
- 二值图像处理模块,
- 用于存储索引(已分割)图像的标签图像,
- 3D与2D堆栈处理模块,
- 核内与核外处理,
- 如何计算测量分析,
- 使用过滤器精简已分割的数据,
- 使用筛分来解释你的测量结果,
- 脚本编写。
这些概念可以通过无数种方式扩展,以应对新的挑战。试着把你的图像处理问题与这个简单的工作流联系起来:
- 图像处理,使数据更易于二值化,
- 通过阈值处理或顶帽等工具进行二值化(有时结合两种技术),
- 打标签,为所有不相连的对象建立索引,
- 为所有已索引的对象测量关键属性,
- 通过视觉检查以及过滤或筛分来分析所测量的数据。
本介绍重点说明了如何使用Avizo对3D数据执行复杂的分割和分析,但除了这里介绍的内容之外,还有更多的处理操作和测量方法。
Avizo为你提供了这套广泛的工具集,以便你能够针对面临的任何处理挑战选用合适的工具。
6.2 示例1:测量催化剂
本教程说明了使用Avizo的更多技术:
- 使用掩膜隔离感兴趣的对象。
- 使用距离图。
- 使用图像算术和分布直方图。
要学习本教程,你应该已经阅读了第一个教程第6.1节——高级图像处理与定量分析入门,并熟悉Avizo的基本操作。显示模块的可见性可以用Avizo的常规方式管理:单击Project View中图标里的橙色方形按钮,或Properties面板中模块标题旁边的按钮。关于Project Views,请参阅第10.1.9节。
本示例中使用的3D图像是通过显微断层扫描采集的:一个近似球形的载体包含催化剂和孔隙。催化剂在图像中呈暗色(低强度体素)。孔隙和背景呈亮色(高强度体素)。中间的灰度值对应于载体。
本示例的目标是获得催化剂体素与背景(外部)之间的距离分布。这里的一个难点在于,外部强度与孔隙内部的强度接近或相同,这就无法使用简单的阈值处理来隔离外部。此外,某些孔隙与外部相连,这就无法使用像Segmentation Editor的magic wand魔术棒或Reconstruction from Markers模块那样的”漫水填充”(flood fill)方法。
注意:在本教程中,你将会看到关于如何管理任意感兴趣区域的一些提示。另一个类似问题的常见示例是:把岩芯样品的孔隙空间与岩芯外部隔离开,以便例如计算岩石的孔隙率。
图 6.20:催化剂与孔隙的显微断层扫描图像。
整个过程被拆分为若干步骤/小节,描述了一个逐步的测量工作流:
Object Detection and Masks对象检测与掩膜More about Region of Interest and Masks关于感兴趣区域与掩膜的更多内容Using Distance Map使用距离图More about Distance Maps关于距离图的更多内容Measurement Distribution测量分布
6.2.1 对象检测与掩膜
- 首先从
AVIZO_ROOT目录加载data/tutorials/image-processing-advanced/Catalyst.am。 - 如果
Auto-Display被禁用,请把一个Ortho Slice模块附加到Project View中的Catalyst.am图像图标以显示该图像。
加载项目data/tutorials/image-processing-advanced/CatalystDistribution-1-LoadData.hx即可完成本教程步骤(见图6.21)。
现在,你可以开始本示例中用于检测对象的第一个步骤:thresholding阈值处理,然后是closing闭运算,以便”填充”对象并准备掩膜。下一节将给出更多关于为你的数据创建掩膜和任意感兴趣区域的可能方法的提示。
阈值处理Thresholding
- 把一个
Interactive Thresholding模块附加到Catalyst.am模块。
为了寻找合适的阈值,你可以直接更改Intensity Range端口。根据Preview Type端口的设置,2D或3D预览会被交互式地渲染。记住要通过更改Preview Slice Number和Preview Orientation端口来检查整个体积。其他Avizo模块对这项任务也会有帮助(参见第3.2节 可视化3D图像)。
图 6.21:初始图像。
这里,把图像在0到225之间进行阈值处理会得到一个二值图像,其中:
- 强度级别 = 1 -> 载体或催化剂(材料),
- 强度级别 = 0 -> 孔隙或外部背景。
图 6.22:项目视图。
应用Interactive Thresholding模块将创建一个二值图像(只有内部和外部材料的图像标签),如下所述:
- 在
Project view中选择Interactive Thresholding模块,然后使用该端口的滑块手柄或文本区域,把Intensity Range值更改为0-225的范围。 - 通过拖动该端口的滑块,把
Preview Slice Number端口更改为45。 - 按下
Apply按钮开始处理。 - 勾选
Preview Type端口中的3D以获得3D预览,取消勾选3D则撤销该预览。 - 把当前链接到
Catalyst.am的Ortho Slice附加到所得到的图像上——方法是单击并把Ortho Slice的连接线拖动到Catalyst.thresholded;默认会选择一个合适的颜色图。 - 通过单击
Project view中该模块的橙色方块,隐藏Interactive Thresholding的预览。
加载项目CatalystDistribution-2-InteractiveThresholding.hx即可完成本教程步骤(见图6.23)。
图 6.23:Thresholding之后的材料二值图像。
形态学:闭运算对象Morphological: Closing object
为了检测对象的形状,你现在可以应用形态学模块。数学形态学模块是基于形状和尺寸准则的变换。
应用于二值图像的形态学Closing模块会给出另一个二值图像,其中:
- 对象内部的小孔洞被填充,
- 对象边界被平滑,
- 靠近的对象被连接起来。
Closing模块实际上先对二值化区域进行膨胀,然后进行腐蚀:直观地说,膨胀填充孔洞并重新连接分离的区域,然后腐蚀恢复原始的外部形状。
图 6.24:Closing模块对二值图像的作用。
以下是填充对象孔隙的步骤:
- 把一个
Closing模块附加到Catalyst.thresholded。 - 为了填充对象内部的所有孔洞,你必须把
Size端口设置为6(结构元素的大小)。此类特定数值可以通过几次尝试找到,例如通过滑动Ortho Slice来检查整个体积。 - 然后按下
Apply按钮以创建所得到的二值图像。 - 把已有的
Ortho Slice附加到Catalyst.closing。
加载项目CatalystDistribution-3-Closing.hx即可完成本教程步骤(见图6.25)。
图 6.25:Closing之后的对象二值图像。
注意:你可能在上图右侧注意到的伪影,是由于膨胀太靠近图像边界所致。为防止这种情况,对象周围的背景边界应大于闭运算的尺寸。使用图像Crop Editor可以很容易地解决这个问题:在本示例中,你可以把Adjust设置为10,然后取消勾选Add mode的Replicate,并把Pixel value设置为背景强度(即0);然后按下Enlarge按钮就会添加一个10体素的边界。不过,在本教程中你可以忽略这个扩大步骤。
6.2.2 关于感兴趣区域与掩膜的更多内容
把测量或处理限制在数据的一个子集上,往往是必要或有用的。
如果该子集是一个轴对齐的长方体,你可以使用以下工具:
- 许多模块支持
Region Of Interest(ROI)输入。你可以把一个ROI Box模块附加到你的数据,然后把显示或计算模块的ROI输入连接到该ROI Box模块。 - 图像
Crop Editor可以裁切或扩展你的数据。 Extract Subvolume模块在内存中复制你数据的一部分,可以进行下采样。
如果你需要任意形状的掩膜或感兴趣区域——例如圆柱形ROI——可以使用以下工具:
- 每个二值图像的弹出菜单中归入
Image Processing/Image Morphology的许多模块,都可用于创建或组合上面所示的那种掩膜。其他示例包括Convex Hull(逐切片应用)、Fill Holes、Reconstruction from Markers。 Volume Edit模块用于使用像圆柱这样的交互式工具修改体积。它也可以通过脚本使用。Segmentation Editor有许多有用的工具可用于快速创建掩膜,例如画笔、成形套索、Selection/Interpolate。
6.2.3 使用距离图
第二步是计算催化剂的距离图。下一节将给出关于距离图的更多提示。
在二值图像上应用距离图算法会得到一个灰度图像,其中每个体素的强度代表以体素为单位到对象边界的最小距离。对于对象距离图的给定体素强度:
- 强度级别 = 0 -> 背景
- 强度级别 = 1 -> 对象包络
- 其他低强度级别 -> 对象中靠近对象包络的部分
- 其他高强度级别 -> 对象中远离对象包络的部分
现在,你可以按如下方式创建这个距离图:
- 把一个新的
Chamfer Distance Map模块附加到Catalyst.closing。 - 把
Interpretation设置为3D并单击Apply。 - 把一个新的
Ortho Slice附加到结果(Catalyst.distmap)以查看对象的距离图。
加载项目CatalystDistribution-4-DistanceMap.hx即可完成本教程步骤(见图6.26)。
图 6.26:对象距离图。
掩膜处理Masking
- 对初始图像
Catalyst.am应用一个Interactive Thresholding模块。 - 把
threshold设置在0到100之间,并单击Apply。
这会得到一个二值图像(Catalyst2.thresholded),其中:
- 强度级别 = 1 -> 催化剂,
- 强度级别 = 0 -> 载体、孔隙或背景。
加载项目CatalystDistribution-5-InteractiveThresholding.hx即可完成本教程步骤(见图6.27)。
掩膜处理将用于计算催化剂的受限距离图。掩膜操作以一个灰度图像作为第一个输入,一个二值图像作为第二个输入(掩膜图像),并提供一个灰度图像作为输出,其中:
图 6.27:催化剂二值图像。
- 掩膜图像中每个黑色体素在输出图像中被设置为0,
- 掩膜图像中每个蓝色体素被设置为来自灰度图像的初始级别。
用催化剂图像对距离图图像进行掩膜处理,会得到一个灰度图像,其中:
- 每个非零强度代表催化剂的一个体素,
- 该强度值等于以体素为单位到对象包络的距离。
你可以按如下方式创建这样的掩膜:
- 应用一个
Mask模块:第一个输入是Catalyst.distmap(距离图图像),第二个输入是先前得到的催化剂二值图像Catalyst2.thresholded。 - 把一个新的
Ortho Slice附加到结果以查看对象的距离图。为获得更好的渲染效果,你可以把Ortho Slice的颜色图范围映射到Catalyst.masked的完整范围(min-max):在Colormap端口中把最小-最大值设置为0…93。
加载项目CatalystDistribution-6-Masking.hx即可完成本教程步骤(见图6.28)。
图 6.28:催化剂距离图。
6.2.4 关于距离图的更多内容
距离图(也称为距离变换)是许多图像处理技术的强大工具。为了更快的计算速度,所计算的距离可能是一个离散近似(倒角图chamfer map)。Avizo提供了多个版本的距离图,你可以针对特定用途进行检查。大多数可用的距离模块位于模块弹出菜单的Image Processing子部分:
Chessboard Distance Map(2D/3D棋盘距离),2D/3D Chamfer Distance Map(棋盘和对角线),Geodesic Distance Map(基于掩膜来隐藏特定部分),2D/3D Closest Boundary Points,- 作为一个特殊情况,
Image Processing/Image Morphology组中的Propagation Distance及相关模块, Distance Map,使用倒角距离或欧氏距离,Image Processing/Distance Maps/Distance Map for Skeleton。有关距离图和3D骨架化的更多细节,请参阅第7.5章 骨架化用户指南。Image Processing/Distance Maps/Distance Map on Disk Data,只能操作旧版LDD磁盘数据。Propagation Distance(已浏览体素距离)被归入模块弹出菜单的Image Processing/Image Morphology子部分。
6.2.5 测量分布
你计算了一个基于体素单位测量距离的倒角图。为了获得一致的结果,你必须考虑体素校准:把距离图像乘以体素大小,可将图像强度转换为公制系统。
模块组Image Processing/Arithmetics Operations可用于两个图像之间或一个图像与一个常数之间的运算。
提示:Arithmetic模块也提供了一种灵活的方式来对图像执行计算。
现在你可以获得距离的分布:
- 在
Catalyst.masked上使用Multiply by Value作为第一个输入(Input Image 1端口)。把Value端口设置为5(假设体素大小为5微米)。 - 单击Apply以得到
Catalyst.mult作为结果。 - 在
Project View中选择该图像图标时,可通过Properties面板中显示的Info端口获取最大值。 - 把一个
Histogram模块附加到Catalyst.mult,以针对每个灰度级i计算并绘制强度为i的体素数量。每个级别的点数将以直方图的形式绘制出来。把Range端口设置为{3,400},把Max Num Bins端口设置为80。取消勾选Options端口中的logarithmic,然后单击Apply以显示直方图窗口。使用File菜单,你可以对直方图拍摄快照,或把直方图数据保存到csv文件。 - 在催化剂距离图上应用
Histogram模块,会生成一个图形,显示位于距对象包络给定距离处的催化剂体素数量。
加载项目CatalystDistribution-7-Histogram.hx即可完成本教程步骤(见图6.29)。
图 6.29:距离的分布。
6.3 示例2:分离、测量与重建
本教程给出了关于对象分离和几何信息提取的更多细节:
Principle of the Watershed Algorithm, an essential tool for image processing分水岭算法的原理——图像处理的一个基本工具Prior Segmentation前置分割Object Separation using Watershed step-by-step使用分水岭逐步进行对象分离Separation Troubleshooting分离故障排查Filtering Individual Objects过滤单个对象Geometry Reconstruction几何重建
要学习本教程,你应该已经阅读了第一个教程第6.1节——高级图像处理与定量分析入门,并熟悉Avizo的基本操作。
本示例中使用的3D图像是用泡沫的若干切片数据生成的。
本示例的目的是以比Getting Started示例更详细的方式来隔离和量化这些气泡,以便更好地控制结果。在许多情况下,分离对于补偿过低的图像分辨率、噪声或整个图像中的强度变化是必不可少的。
6.3.1 分水岭算法的原理
分水岭算法是一种强大的方法,在图像处理中有许多应用,例如用于自动的对象分割或分离。
该算法模拟了在2D或3D图像中从一组已标记区域开始的漫灌过程。它根据一个优先级图来扩展这些区域,直到这些区域到达分水岭线。该过程可以看作是在一片地貌中逐步浸没。
图 6.30:泡沫图像。
图 6.31:分水岭。
该算法依赖于两个输入:
- 一个标签图像,包含被标记的标记区域(marker regions),用作漫灌的种子区域。在过程结束时,被分离出的对象数量将与被不同标记的标记数量一样多。这些标记就像是我们想要获取其汇水区域的河流。
- 一个灰度图像,起到地貌高度场或海拔图的作用,它控制着漫灌的推进过程,并最终决定分水岭分隔线的位置。这些分隔线位于我们地貌中山谷之间的脊线上。
通过仔细选择这两个输入(标记和优先级图),可以实现不同的应用。
一个示例是Getting Started教程中用于分离泡沫孔隙的Separate Objects模块。它首先计算一个距离图(有关距离图的更多内容,请参阅第6.2节——测量催化剂教程)。该距离图为分水岭过程提供了优先级图输入。距离图的极大值区域——孔隙最内部的区域——提供了用于分水岭的标记输入。该过程将在下一节中详细描述。
6.3.2 前置分割
首先,你需要对数据进行分割,以获得孔隙的二值图像。
- 首先从
AVIZO_ROOT目录在Avizo中加载图像堆栈data/images/foam/foam.am。 你可能希望对数据执行某种降噪处理。你通常可以应用像带
3D Interpretation的Median Filter这样的滤波器;中值滤波器是一种非线性数字滤波技术,常用于在去除噪声的同时保留边缘。这种降噪是一个典型的预处理步骤,用于改善后续处理(如分割或边缘检测)的结果。不过,在本示例中你可以跳过这个阶段,直接开始通过阈值处理创建二值图像。
- 把一个
Interactive Thresholding模块附加到项目视图中的foam.am对象。 - 使用
Intensity Range端口滑块的游标,把低阈值和高阈值级别设置为0和38,然后单击Apply。在预览中,你可以看到图像中的分辨率和强度分布不允许你直接分割出分离的孔隙:无论选择什么阈值,某些孔隙仍然是相连的,除非在孔隙中留下过多噪声。 - 把一个
Ortho Slice附加到foam.thresholded以可视化结果。
加载项目WatershedSeparation-1-Thresholding.hx即可完成本教程步骤(见图6.32)。
图 6.32:二值图像。
通过拖动Slice Number端口的游标,你可能会注意到一些伪影孔洞(例如在切片4、5、6中),它们大多与灰度图像上出现的穿过孔隙的较高强度环有关。把一个Ortho Slice附加到灰度图像,并把Mapping type设置为histogram以突出显示这些伪影。隐藏或移除之前的Ortho Slice。
与其对灰度图像甚至采集过程进行上游校正,在某些情况下,校正二值图像(例如填充孔洞)可能更为有效。
可选地,建议填充待分离对象内部的孔洞,因为基于距离图的分离方法可能对这些伪影孔洞敏感。
- 你可以把一个
Fill Holes模块附加到foam.thresholded。 - 把
Interpretation设置为3D并单击Apply以得到foam.filled作为结果。
6.3.3 使用分水岭逐步进行分离
从二值图像开始,你现在可以继续进行孔隙分离。你将计算一个距离图,从距离图的极大值区域创建标记,然后应用快速分水岭算法。
- 把一个
Chamfer Distance Map模块附加到二值图像对象(foam.thresholded)。 - 把
Interpretation设置为3D并单击Apply以得到foam.distmap作为结果。 - 附加一个
Ortho Slice来检查结果。孔隙内的每个体素都获得了一个对应于它到黑色背景(泡沫)距离的值。
加载项目WatershedSeparation-2-DistanceMap.hx即可完成本教程步骤(见图6.33)。
- 附加一个
H-Maxima模块。把Contrast端口保留为默认值(4)并单击Apply。 - 把一个
Ortho Slice附加到foam.hMaxima,并把Transparency类型设置为Alpha。
加载项目WatershedSeparation-3-MergedMaxima.hx即可完成本教程步骤(见图6.34)。
该模块创建一个二值图像,包含输入距离图图像的区域极大值,这些极大值在作为参数给出的对比度变化范围内被”合并”。由于使用距离图作为输入,因此结果是对象内部最内层区域的集合。
分水岭算法要求为最终被分离出的每个区域提供一个唯一的标签。输入图像中具有相同值的两个区域会被合并。
- 把一个
Labeling模块附加到foam.hMaxima并单击Apply。 - 把一个
Ortho Slice附加到foam.labels,并把Transparency类型设置为Alpha。
加载项目WatershedSeparation-4-Markers.hx即可完成本教程步骤(见图6.36)。
图 6.33:倒角距离图。
图 6.34:距离图的合并极大值。
距离图还需要被反转,因为分水岭算法会朝着输入优先级图(即地貌海拔)数值增加的方向扩展标记。
- 对
foam.distmap应用一个NOT模块。 - 把一个
Ortho Slice附加到foam.not。
图 6.35:在一个图像剖面样本上展示的H-Maxima(合并极大值)原理。
图 6.36:由Labeling模块创建的标记。
加载项目WatershedSeparation-5-ReversedDistanceMap.hx即可完成本教程步骤(见图6.37)。
现在你可以计算分水岭分隔线的图像了。
- 把一个
Marker-Based Watershed模块附加到反转距离数据对象(foam.not)。把foam.labels设置为Input Label Image,把Type设置为Watershed。按下Apply。 - 把一个
Ortho Slice附加到结果以查看分水岭线。把Alpha或Binary设置为Transparency类型,以便把线叠加到灰度图像上。 - 把一个
Ortho Slice附加到foam.am。把Colormap最小值设置为0。
加载项目WatershedSeparation-6-SeparationLines.hx即可完成本教程步骤(见图6.38)。
图 6.37:使用Line Probe模块显示的反转距离图。
图 6.38:分水岭分隔线。
为完成分离,你可以从孔隙的二值图像中减去分隔线:
- 对
foam.thresholded应用一个AND NOT Image模块作为第一个输入,foam.watershed作为第二个输入。 - 附加一个
Ortho Slice以查看结果foam.sub。
加载项目WatershedSeparation-7-SeparatedPores.hx即可完成本教程步骤(见图6.39)。
图 6.39:已分离的孔隙。
图 6.40:分水岭分离工作流。
6.3.4 分离故障排查
你可能会在结果中看到未被分离的孔隙,或者相反,看到不需要的孔隙分离。由于该方法基于几何距离准则,因此分离对于凸的、近球形的对象效果最好。以下是一些改进分离的指导原则:
- 如果你查看某个特定切片上的标记图像,某些孔隙中的标记可能看起来缺失了,这只是因为它们位于3D中的其他位置。
- 不过在某些情况下,标记可能真的缺失了——因为从距离图区域极大值的角度看,两个对象被认为是合并的,那么分离也会缺失。你可以尝试降低
H-Maxima(或Separate Objects)的对比因子,以使标记更小、更加分离。 - 如果某个分离缺失了,意味着某个标记缺失了。
- 分离”线”可能看起来太粗或不正确,这只是因为你所查看的切片与某个分离面有些相切。
- 如果对象形状是非凸的,可能会出现不需要的分离:这会导致多个局部极大值,因而产生被分开的对象。对象中的小凹陷可能导致跨越它的分离。一种解决方案是增大
H-Maxima的对比因子,以使标记更大并被合并。 - 基于距离图的分离可能会遵循最短路径,即直线而不是期望的形状。这是因为分水岭是由距离驱动的:分离由几何准则驱动。
Chamfer Distance Map是一个离散的倒角图。你可以改用更精确的欧氏距离图(参见对象弹出菜单中Image Processing/Distance Maps下引用的Distance Map或其他距离图类型);不过在大多数情况下这影响甚微。
作为主要的指导原则,被分离出的对象数量将与标记数量一样多。
如果结果不令人满意,可能有两种情况:
- 对象数量不令人满意:在这种情况下,你必须在标记上下功夫。
- 分离线不令人满意:你必须在优先级图(即”地貌高度场”)上下功夫。
对于标记,Separate Objects的解决方案是在距离图上做工作。在这种情况下,你考虑的是几何形状(大多是最内部的区域),但标记也可以通过其他方式获得。有时,如果颗粒的中心更暗或更亮,使用灰度图像(强度)可能会很有意义。
如果对象相当均匀,你必须坚持使用几何信息。一种方法是调整H-Maxima的对比因子设置。基本上,H-Maxima参数对应于两个极大值之间的最小深度。用地理来类比:即鞍部与峰顶之间的高度差,这样你就只保留两个不同的峰。
在一个显示了不应被分开的对象的裁剪数据集上进行试验,并调整设置以确保每个对象内部只有一个标记,这样可能更容易。一种可能的改进方法是先做一次H-Maxima,然后用距离图的阈值图像对结果进行掩膜处理。这可以避免为小的连通部分保留标记,因为你只保留了对应于距离图极大值且距离值最小的标记。
一旦对标记(对应于对象数量)感到满意,你还可以改进分离线。同样,上面描述的过程(Separate Objects)依赖于几何形状,你也可以把强度或强度梯度作为深度的函数来使用。把几何信息和强度信息结合起来可能有点困难。在某些情况下,你可以通过距离函数与强度函数的组合来实现,例如使用Blend with Image、Blend with Value或Arithmetic模块,以某种方式模拟人眼组合这两种信息的复杂方式。请注意,分水岭寻找的是峰值分隔线,因此是局部极大值——就像距离图那样,你可能需要使用函数的负值来适应这种情况。
另一种基于某些准则来改进分离的强大技术是过滤由分水岭方法添加的分隔:
- 隔离出被添加的分隔:例如使用
AND NOT Image模块获得原始图像与分离后图像之间的差异, - 用
Labeling模块为这些分隔打标签, - 使用
Label Analysis模块对每个分隔进行一些测量,例如,某个分隔区域内距离图值的最大值可以指示该分隔是否深入到了对象内部过深的位置。
你将在接下来的教程中找到更多分水岭的应用示例,它们使用不同的方法来创建标记和优先级图。
6.3.5 过滤单个对象
你可以对已被分离的单个对象进行测量。
首先,你需要把已分离孔隙的二值图像转换为一个标签图像,其中每个孔隙都被唯一标识:
- 对已分离孔隙图像(
foam.sub)应用一个Labeling模块。 - 你可以附加一个
Voxelized Rendering模块来查看已打标签的图像。
加载项目WatershedSeparation-8-SegmentedPores.hx即可完成本教程步骤(见图6.41)。
其次,你想要过滤掉不需要的小对象,最后测量单个孔隙的体积。
为了移除小对象,你可以使用Getting Started教程中描述的测量和测量过滤器。这涉及两个步骤:
图 6.41:已分割的孔隙。
- 把一个
Filter by Measure模块附加到标签图像。 - 把
foam.am设置为Input Intensity Image。 - 在
Measure端口的大量可用测量列表中选择Volume3d。 - 把
Number of Objects端口设置为15,以只保留15个最大的体积。按下Apply。 - 可以使用
Voxelized Rendering或Boundary Rendering模块来查看最大的体积。 - 你也可以使用
Label Analysis模块在Tables面板中显示和管理特定的测量。
在接下来的某个教程(第6.4节——进一步的图像分析)中,你将找到更多测量和高级分析的应用示例。
加载项目WatershedSeparation-9-FilteredPores.hx即可完成本教程步骤(见图6.42)。
图 6.42:已分割并过滤的孔隙。
6.3.6 几何重建
最后,你可以重建气泡的几何形状:
- 把一个
Generate Surface模块附加到结果标签图像(foam3.labels),取消勾选Border端口中的Adjust Coords,然后按下Apply以生成foam3.surf。 - 通过附加一个
Surface View模块来显示所得到的曲面。
加载项目WatershedSeparation-10-Reconstruction.hx即可完成本教程步骤(见图6.43)。
由于结果图像已经包含用于材料标签的整数值,因此它可以直接用于曲面重建。其他图像类型可能需要转换为Avizo标签数据(例如使用Convert Image Type模块)。请注意,Generate Surface模块可以处理超过256个标签。
由于曲面多边形数量庞大,用Surface View模块显示所得到的曲面可能会非常慢。在这种情况下,建议先对曲面进行简化,以便在你的硬件上更快地显示。
图 6.43:已过滤孔隙的曲面重建。
Avizo还允许你把曲面导出为各种文件格式,或者生成并导出一个四面体模型,例如适用于使用某些外部求解器进行有限元模拟。
对应于整个本教程的一个演示脚本可以在以下位置找到:data/tutorials/image-processing-advanced/PorositySurfaceReconstruction.hx
它使用位于data/tutorials/image-processing-advanced中的一个脚本对象来自动化该处理过程。它与图6.40中描述的工作流相匹配。
脚本加载后,你可以选择性地更改若干端口值,然后单击Action端口的Apply按钮以开始处理。
6.4 示例3:进一步的图像分析——泡沫中孔径的分布
本示例展示了如何计算泡沫样品中孔径的分布,以及如何定义一个自定义测量来计算孔隙的球形度。
要学习本教程,你应该已经阅读了第6.1节——高级图像处理与定量分析入门,并熟悉Avizo的基本操作。
本节分为以下步骤:
pore detection孔隙检测,pore post-processing孔隙后处理,custom measure group definition to determine the distribution of pore diameters自定义测量组定义,用于确定孔径的分布,custom measure definition to compute the sphericity of pore自定义测量定义,用于计算孔隙的球形度。
本示例中使用的图像是通过显微断层扫描采集的。它代表由材料和孔隙组成的泡沫。孔隙在图像中呈暗色(低强度体素)。材料呈明亮色(高强度体素)。
- 首先从
AVIZO_ROOT目录加载data/tutorials/image-processing-advanced/FoamPoro.am。 - 如果
Auto-Display被禁用,请连接一个Ortho Slice以在3D查看器中可视化数据,如图6.44所示。
6.4.1 第一步:孔隙检测
- 把一个
Interactive Thresholding连接到数据。 - 使用
Threshold端口把图像在0到50之间进行阈值处理。 - 按下Apply。
- 隐藏
Interactive Thresholding模块,并把Ortho Slice连接到输出FoamPoro.thresholded。
如图6.45所示,把图像在0到50之间进行阈值处理会得到一个二值图像,其中:强度级别1 = 孔隙率,强度级别0 = 载体(材料)。
6.4.2 第二步:孔隙后处理
应用于二值图像的形态学Opening算子会给出另一个二值图像,其中:小对象被移除,对象边界被平滑,某些对象可能被断开连接。
图 6.44:泡沫孔隙的初始显微断层扫描图像。
- 把一个
Opening模块连接到二值图像。 - 把
Size [px]端口设置为1。 - 按下Apply。
- 把
Ortho Slice连接到输出FoamPoro.opening。
在先前计算出的孔隙二值图像上应用形态学开运算,会得到一个噪声和伪影被减少的滤波图像,如图6.46所示。
Separate Objects模块检测出分隔团聚颗粒的曲面。这些曲面从初始图像中被减去,见图6.47。
- 把一个
Separate Objects模块连接到已滤波的二值图像。 - 把
Marker Extent设置为1。 - 按下Apply。
- 把
Ortho Slice连接到输出FoamPoro.separate(见图6.48)。
图 6.45:孔隙。
6.4.3 第三步:自定义测量组定义,用于确定孔径的分布
Avizo的Label Analysis模块允许为3D图像的每个颗粒计算一组测量。执行完单个分析后,可以绘制给定测量的直方图,以便产生该测量分布的表示。
- 把一个
Label Analysis模块连接到已分离的二值图像。 - 把
FoamPoro.am设置为Intensity Image。
在该模块的Measures端口中,basic是一组预先选定的原生测量。可能会出现这种情况:你并不需要basic测量组中的所有测量,或者你想在分析表中捆绑一组不同的测量。对于这些情况,你可以创建自己的测量组。
对于给定的颗粒,等效直径测量计算的是相同体积的球形颗粒的直径。因此等效直径由以下公式给出:
- 按下
Measures端口的配置按钮(带3个点的按钮)。会打开一个面板用于选择测量组(见图6.49)。
图 6.46:已滤波的孔隙。
图 6.47:孔隙后处理。
图 6.48:已分离的孔隙。
- 通过按下测量组选择器旁边的专用按钮(1)来创建一个新的测量组。
- 在弹出窗口中,把新组命名为
diameter并按下OK。 - 在原生测量列表(2)中选择
EqDiameter,并使用箭头按钮把它添加到该组(3)。 - 按下OK。
新的diameter组只包含等效直径测量,现在已在Label Analysis模块的Measures端口中被选中。
- 按下
Apply按钮。
Project View中会创建一个新的label图像数据对象FoamPoro.label,并显示Tables面板,其中以电子表格样式的表格展示结果:分析结果FoamPoro.Label-Analysis也会在Project View中被创建(见图6.50和图6.51)。
图 6.49:测量组的选择。
图 6.50:分析项目。
- 在下方的电子表格中选择
EqDiameter列。 - 按下工具栏中的直方图按钮。
一个窗口会打开,显示EqDiameter直方图,如图6.52所示。
加载项目data/tutorials/image-processing-advanced/CustomerMeasurements-1-EquivalentDiameter.hx即可完成上述步骤。
图 6.51:分析电子表格。
图 6.52:EqDiameter直方图。
6.4.4 第四步:自定义测量定义,用于计算孔隙的球形度
Avizo提供了一组预定义的原生测量,但也可以保存用户自定义的测量。
图 6.53:自定义测量定义按钮。
- 选择
Label Analysis模块。 - 在
Properties面板中,按下Measures端口的配置按钮(带3个点的按钮)。 - 如果尚未选中,请选择
diameter组。 - 通过按下专用按钮创建一个自定义测量(见图6.53)。
- 把它命名为
Sphericity并按下OK。
Measure Edition面板会打开。
球形度是衡量一个对象有多接近球形的测量,其表达式为:
其中V是颗粒的体积,A是它的表面积。
它是球体表面积(该球体与给定颗粒具有相同体积)与该颗粒表面积之比。球体的球形度为1,任何非球体的球形度都小于1。几种典型几何体的球形度参考值:正二十面体≈0.939、正十二面体≈0.910、正八面体≈0.846、立方体≈0.806、圆柱(h=2r)≈0.874、圆锥(h=2√2r)≈0.794、圆环(R=r)≈0.894、正四面体≈0.671。
用Avizo的测量变量表达该公式为:
1 | (pi**(1/3)*(6*Volume3d)**(2/3))/Area3d |
- 在面板的专用字段中输入球形度公式(见图6.54)。
- 按下
Close。 - 在自定义测量列表中选择
Sphericity,用箭头按钮把它添加到diameter组。 - 按下OK。
- 按下
Label Analysis模块的Apply按钮。
分析面板会更新,显示EqDiameter和Sphericity两项测量。
图 6.54:测量编辑。
提示:使用自定义组只选择需要的测量,有助于加快分析过程。
注意1:自定义测量组和自定义测量这类用户定义数据,会作为本地设置在工作会话结束时持久保存,因此重启Avizo后仍可取用。这些自定义数据也会保存在项目脚本中。
注意2:对于小孔隙(即由很少体素组成的孔隙),计算出的球形度可能大于1。这是因为Area3d测量使用弦近似计算(这通常能更好地近似面积),而Volume3d测量则不使用任何近似。
加载项目data/tutorials/image-processing-advanced/CustomerMeasurements-2-Sphericity.hx即可完成上述步骤。
6.5 示例4:进一步的图像分析——泡沫中材料的平均厚度
目标:计算泡沫样品中材料的平均厚度。数据仍用data/tutorials/image-processing-advanced/FoamPoro.am(孔隙暗、材料亮)。
算法分4步:孔隙检测 → 检测分隔面 → 材料的距离图 → 计算材料平均厚度。
6.5.1 孔隙检测
复用6.4的前几步:阈值处理 → 形态学开运算 → Separate Objects,得到已滤波并分离的孔隙二值图像FoamPoro.separate。
对应项目:PorosityThickness-1-SeparationObjects.hx(见图6.55、6.56)
图 6.55:泡沫孔隙的显微断层扫描图像。
图 6.56:孔隙率的二值图像。
6.5.2 检测分隔面
目标:找出穿过材料、且与两个孔隙等距的那些面。
Influence Zones模块(影响区骨架,SKIZ):输入二值图像,输出二值图像,其中蓝色体素更靠近区域中心的那个对象,黑色体素则与至少两个最近对象等距。
- 对
FoamPoro.separate应用Influence Zones→ 得到FoamPoro.zones - 再对
FoamPoro.zones应用NOT(反转)→ 得到FoamPoro.not
Influence Zones + NOT 的组合,就给出了穿过材料分隔各孔隙的那些面的二值图像。
对应项目:PorosityThickness-2-InfluenceZones.hx、PorosityThickness-3-SeparationSurfaces.hx
图 6.57:由Influence Zones得到的图像/骨架。
图 6.58:分隔面图像。
6.5.3 材料的距离图
- 对孔隙二值图像
FoamPoro.separate应用NOT→ 得到材料二值图像FoamPoro2.not(黑=孔隙/背景,蓝=材料) 对材料二值图像应用
Chamfer Distance Map(Interpretation设为3D)→FoamPoro2.distmap距离图含义:强度0=孔隙或背景,强度1=材料包络,低强度=靠近包络的材料,高强度=远离包络的材料。
对距离图应用
Mask,把Input Binary Image设为分隔面图像FoamPoro.not掩膜后得到的灰度图像中,每个非零强度代表一个分隔面体素,其强度值等于两个最近对象之间距离的一半。
对应项目:PorosityThickness-4-DistanceMap.hx、PorosityThickness-5-Mask.hx
图 6.59:材料二值图像。
图 6.60:材料的距离图。
图 6.61:分隔面的距离图。
6.5.4 计算材料平均厚度
其中:
- $V_{size}$ = 以微米为单位的体素尺寸
- $\Sigma_i$ = 分隔面距离图图像中体素强度之和
- $NbVoxSep$ = 二值分隔面图像中的体素数量
取值方法:
Volume Fraction模块(3D Interpretation)的Label Voxel Count列给出3D图像中每个标签的已标记体素数。二值图像只有一个标签,所以该列返回$NbVoxSep$。Intensity Integral模块(3D Interpretation)的Volume列给出3D图像中体素强度之和,即$\Sigma_i$。
对应项目:PorosityThickness-6-Thickness.hx
6.6 分水岭分割
为什么需要分水岭:由于图像采集和重建带来的伪影,简单阈值分割常常不准甚至错误,尤其在相变过渡区域:
- 边缘被噪声和部分体积效应模糊时,正确或唯一的阈值难以准确确定。部分体积效应源于采集分辨率限制,它模糊了相与特征之间的过渡(即单个体素并不只对应一种均匀物理体积)。
- 分割超过两相时,高低强度相之间的过渡可能引入不需要的中间”镀层”(coating)伪影相。
- 图像上光照或强度的变化会导致不同区域需要不同阈值。
本教程涵盖两种Avizo工具:Segmentation Editor中的交互式分水岭工具;以及用于多材料/多相的Watershed Segmentation向导。
6.6.1 在Segmentation Editor中用分水岭工具分割砂样
数据:压实二氧化硅砂样(部分水饱和),X射线显微断层扫描,体素约11.2 µm。数据集中可清楚区分空气、水、硅酸盐颗粒三相;本教程只关心孔隙空间 vs 颗粒。
图 6.62:部分水饱和的压实二氧化硅砂样。
准备种子标记
- 打开
data/sandpack/sandpack128-filtered.am(该图像已对原始图像应用Non-Local Means Filter降噪) - 用
Edit New Label Field进入Segmentation Editor - 在
Materials列表中把Inside重命名为Pore_Space;再Add一个新材料并命名为Grains - 选
Threshold Tool,Masking边界设为0-6000,勾选All slices,单击Select Masked Voxels,选中体素显示为红色 - 选中
Pore_Space,在Selection组按大”+”(或按A键)指派 - 把
Threshold Tool边界改为10000-51000,All slices保持勾选,Select Masked Voxels - 选中
Grains,同样按”+”指派
关键思路:这里的阈值看起来偏保守、会留下未选区域,但更高的阈值可能捕获颗粒内部噪声或跨越真实相界。这里要的是一个“安全”选择,后续靠分水岭工具把它填充并扩展到真实相界。
图 6.63:在Segmentation Editor中粗略指派Pore_Space相。
图 6.64:在Segmentation Editor中粗略指派Grains相。
执行分水岭
- 选
Watershed Tool - 在
Landscape image选项中按Create a new gradient image—— 用快速Canny方法计算梯度幅值图像,作为控制标记扩展的地貌图像 Marker列可选择哪些材料用作标记,此处保持全部勾选Output catchment basins保留默认side-by-side,以获得连续相邻的标签(而非被外部体素分隔的标签)- 按
Apply and create new label field
新标签图像中,标记被扩展到地貌图像的边缘(即梯度幅值的极大值处)。边界会基于种子标记和梯度图像,落在相强度过渡之间的最优位置。
图 6.65:应用Segmentation Editor分水岭。
实用提示
- 结果不满意时,删除当前标签即可取消分水岭步骤,回到原始标记标签,调整或补全后再次应用(见图6.66)。
- 用
Label field选择器菜单可在标记与分水岭结果、或不同分割结果之间来回切换。 - 用
Image选择器菜单可临时更改作为背景的分割图像,再切回实际被分割的图像。这也便于把Create a new gradient image创建的地貌图像显示为标签的背景(需要调整数据以获得更好的可视化)。 - 颗粒可能未按预期分离,原因可能是:某些情况下确实存在带中间强度相的固结颗粒;或图像分辨率与部分体积效应所致。后者可能需要改进种子标记(例如用
TopHat工具)和/或梯度(见Image Gradient模块),也可以用Separate Objects来分离颗粒解决。 - 分水岭之后可以继续修饰分割结果,例如用
Segmentation菜单中的Smooth labels...或Remove islands...。
图 6.66:删除当前标签或切换到另一个标签。
6.6.2 用Watershed Segmentation向导分割多相
数据:data/tutorials/chocolate-bar.am,含慕斯、焦糖、巧克力和空气等区域。
先看直接阈值为何不够
- 附加
Auto Thresholding到chocolate-bar.am,Type设为Auto Segment 3 Phases,Criterion设为factorisation(即Otsu方法),Apply,再附加Ortho Slice
结果中一层慕斯(浅蓝)似乎包裹着巧克力顶层(深蓝)—— 但外部与巧克力之间的这层慕斯实际上是部分体积效应造成的伪影(见图6.67、6.68)。调整阈值无法解决该问题。如图6.69所示,图像梯度能给出真实边界的有用线索。因此下一节采用更稳健的、利用图像梯度的方法。
图 6.67:用Auto Thresholding自动分割图像。
图 6.68:简单阈值处理会选中从高强度到低强度过渡区中因部分体积效应产生的不需要的”镀层”体素,这一点由使用Line Probe模块得到的强度剖面显示出来。
图 6.69:Image Gradient的峰值有助于更好地定位相界。
Watershed Segmentation向导
该向导是一个脚本模块,逐步引导完成以下分割流程:
- 定义”保守的”阈值,为空气、慕斯和巧克力确立初始种子标记
- 计算映射材料边缘”锐度”的图像梯度幅值
- 标记空气与钢件之间的锐利边缘,通过掩膜种子标记来防止”镀层效应”
- 完成种子标记朝对象边缘的分水岭扩展
操作步骤:
- 移除
Auto Thresholding模块及其结果数据,Project View中只留chocolate-bar.am - 把
Watershed Segmentation向导模块附加到chocolate-bar.am
选中向导模块后,Properties面板中的端口会显示当前步骤及可选项与参数(见图6.70)。任何步骤都可以回退和修正。
提示:向导会保留中间数据以便回退步骤。对于可能超出可用内存的大数据,可在进入下一步时显示并移除中间数据(用对象菜单的
show),但此后就无法回退步骤了。
- 任何时候都可以移动切片或改变其方向。当前视图是XY切片第147号。
- 第1步:输入要分割的独立材料/相数量 —— 此处设为3(空气(外部)、焦糖、巧克力)。按
Action端口的Apply。 - 第2步:可选地附加一个预先计算好的梯度图像(
Gradient端口)。直接按Apply即可触发图像梯度幅值的计算。 - 第3步(Threshold Gradient Magnitude):用滑块调整阈值,基于上一步得到的梯度来标记锐利边缘。
- 取消勾选
auto复选框以显示Gradient Threshold端口 Gradient Threshold可设为707,用于掩膜空气-巧克力过渡区域- 此时可通过改变切片号或方向来验证标记是否恰当
- 可选中
MainOrthoSlice模块,调整边缘区域的颜色图范围或透明度以便更好地观察边缘(见图6.71) - 然后重新选中
Watershed Segmentation模块(见图6.72),按Action端口的Apply
- 取消勾选
图 6.70:Watershed Segmentation向导第1步。
图 6.71:MainOrthoSlice模块属性。
图 6.72:Watershed Segmentation向导第3步。
接下来三个”Threshold Phase…”步骤用于为每种材料定义标记标签。
- 第4步(相0:空气):此处不追求精确贴合对象边缘,而是要安全地避开可能发生实际材料过渡的区域。慕斯中的孔洞应被标记进相0,否则它们会被归入慕斯区域。范围可用0-435(见图6.73),按
Apply。 - 第5步(相1:慕斯):范围可用629-674。同样,目标是标记给定材料的内部区域(见图6.74)。如果选择更低的下限,空气与巧克力之间的区域会被选中,那就会产生与图6.67相同的伪影。按
Apply。 - 第6步(相2:巧克力):最后,在下一步中你可以设置相2(巧克力)的范围。这里范围可用1041-1910(见图6.75)。应尽可能避开慕斯区域。按
Apply完成该步骤。
图 6.73:Watershed Segmentation向导第4步。
图 6.74:Watershed Segmentation向导第5步。
图 6.75:Watershed Segmentation向导第6步。
在下一步中,你可以触发分水岭计算并得到最终的标签图像结果。
- 直接按
Apply即可触发计算(见图6.76)。 - 你可以通过改变切片编号或方向来验证结果,必要时也可以回到之前的步骤去修改某些参数。
- 然后你就可以移除
Watershed Segmentation模块,这同时也会清理掉那些辅助的显示模块。
图 6.76:Watershed Segmentation向导执行完成。
你可能想去掉被向导分割为标签1的那个空气相。为此,你可以直接使用Subtract Value模块(值设为1),或者使用更通用的Arithmetic模块(表达式为”A-1”)(见图6.77)。
图 6.77:使用Arithmetic模块移除空气相。
更多分割提示More segmentation hints
上面展示的分水岭技术可以用在许多不同的工作流中,以便根据所期望的结果实现要求更高或更具针对性的分割。
例如,在chocolate-bar数据集中,可以略微区分出由焦糖构成的第4个相。图像噪声使得无法通过阈值处理直接标记出正确的区域。不过,先应用一个平滑滤波器(例如Non-Local Means Filter)有助于进一步的分割。使用Watershed Segmentation向导做一次简单的尝试(4个相,梯度阈值100,各相标记阈值为0-435、630-750、800-1210、1250-1459),会得到图6.78所示的结果。
图 6.78:在Non-Local Means滤波之后使用Watershed Segmentation向导的示例。
如果结果不令人满意,可以用分割编辑器Segmentation Editor通过若干步骤来校正或细化该分割。你可以使用分割编辑器交互式地调整标记区域,然后像前面教程中那样应用分水岭工具。
再举一个例子:小孔隙可以被单独分割出来(例如使用Interactive Top-Hat模块),然后用Segmentation Editor或Arithmetic模块与先前的分割结果合并。
一个示例结果是data/tutorials/chocolate-bar-labels-4-phases.am(见图6.79)。
图 6.79:4个相的结果示例。
6.7 关于图像滤波的更多内容
在分割之前,往往有必要降低图像噪声或伪影,并增强感兴趣的特征。数字图像滤波器就是用来增强图像或突出图像特征的工具。这些工具可能基于复杂的算法,其效果和性能会随输入图像和控制参数而变化,因此可能需要通过试验才能达到期望的结果。
在本教程中,你将学习如何:
- 区分Avizo中可用的各类图像滤波器,以及通常常用的主要滤波器。
- 有效地使用图像滤波器来调整参数并比较结果。
- 组合应用多个图像滤波器。
要学习本教程,你应该已经阅读了第一个教程第6.1节——图像处理与分析入门,并熟悉Avizo的基本操作。特别是在入门教程中,你可以看到如何应用一个图像滤波器,以及之后完整的定量分析工作流。
6.7.1 选择图像滤波器
为方便起见,图像滤波器按其主要功能分类整理。你可以在对象弹出菜单的浏览面板中、Image Processing这个顶层类别下浏览各个滤波器类别。你也可以使用对象弹出菜单的搜索工具,按名称快速找到某个滤波器。以下各节简要描述了不同的类别,并推荐了一部分图像滤波器,同时给出一些提示。
6.7.1.1 平滑Smoothing
这是为图像分割准备数据时最重要的一个类别。这些滤波器有助于平滑有噪声的图像。有些用户把它们称为”降噪”滤波器。它们能有效降低噪声,但可能需要小心,以免改变图像中所含的信息,尤其是用于定量分析的场合。Avizo中最常用的平滑滤波器有:
Median Filter中值滤波——一种保留边缘的基本滤波器。对椒盐噪声(散点)非常有效。Bilateral Filter双边滤波——在平滑与保留边缘(特别是锐角)之间取得平衡的滤波器。需要Avizo许可证。Non-Local Means非局部均值(GPU加速)。该滤波器对有噪声的数据极为有效,同时能保留边缘(对白噪声效果最佳)。它通常是有噪声图像的首选。不过它可能非常耗时。缩小搜索窗口会减少计算时间,但也可能降低平滑效率(取决于噪声分布)。最好不要在该滤波器之前应用边缘增强。需要Avizo许可证。Edge-Preserving Smoothing保边平滑——一种保留边缘的扩散滤波器。Anisotropic Diffusion各向异性扩散——一种GPU加速的扩散滤波器。需要Avizo许可证。Curvature-Driven Diffusion曲率驱动扩散——一种可能更好地保留细小结构的扩散滤波器。需要Avizo许可证。
提示:请确保诸如保边扩散的对比度阈值这类滤波器参数是按照你的数据范围来设置的。默认范围通常是针对8位数据的,对于动态范围更高的16位图像,需要大幅提高这些取值。
6.7.1.2 锐化Sharpening
这些滤波器有助于加强边缘处的对比度,让细节看起来更清晰。常用的锐化滤波器有:
Unsharp masking非锐化掩蔽Delineate——同时也起平滑滤波器的作用。需要Avizo许可证。
6.7.1.3 边缘检测Edge detection
这些滤波器会突出不同材料或相之间的边界。它们可以用来直接提取特征轮廓和边缘,或用于基于分水岭的分割(参见教程第6.6.1节——高级分割)。
以下是该类别中最常用的模块:
Sobel Filter——快速的基本边缘检测,可用作梯度幅值的近似。Image Gradient——功能全面的模块,支持快速近似(Canny)或降噪型梯度(Canny Deriche、Gaussian、Sobel、Prewitt)。
6.7.1.4 频域Frequency domain
这些滤波器通过在频域中进行变换来工作。
FFT——傅里叶变换,是许多图像滤波技术的基础模块。Deconvolution——用于对3D光学显微图像进行反卷积的专用模块。
6.7.1.5 灰度变换Grayscale transforms
严格来说,该类别中的模块并不是图像滤波器。上面提到的滤波器通常是通过检查每个像素周围的一个强度值邻域来对像素进行处理的。而灰度变换则独立地作用于各个像素,即不考虑相邻的像素值。以下模块可用于对图像灰度或明暗不均进行全局校正:
Shading Correction、Shading Correction Wizard、Correct Z Drop、Background Detection Correction——这些模块有助于补偿图像中不均匀的背景。Match Contrast——根据一个参考来调整图像的动态范围。需要Avizo许可证。
6.7.2 调整图像滤波器
先在数据的一个有限子集上调整处理工作流,始终是一个好习惯。图像滤波器也不例外,带视觉反馈地交互式调整参数非常有用。以下是一些有助于更快调整图像滤波器的方法:
Crop Editor——移除数据中无用的部分。Extract Subvolume——提取一个子集用于试验。Slice——把某些图像滤波器应用到所显示的切片上。- 图像滤波器的2D/3D解释端口——2D处理(逐切片)通常明显更快,而且可能给出相似的结果。
- 图像滤波器的内核尺寸、迭代次数、窗口大小——滤波器性能往往取决于此类参数。
Filter Sandbox——在某个图像区域上预览滤波器效果。需要Avizo许可证。
以下是Filter Sandbox模块的使用方法。Filter Sandbox对于选择滤波器并调整其参数非常有帮助。请注意,这个便捷的脚本模块只从Avizo中可用的滤波器里提供了最常用的一部分。
- 新建一个项目(Ctrl+N)并打开
data/tutorials/chocolate-bar.am。 - 附加一个
Filter Sandbox模块。你可以在右键对象弹出菜单的Image Processing类别中找到它,或者在搜索字段中键入Filter Sandbox的若干字母。
3D查看器中会显示一个Ortho Slice,上面叠加着一个预览框,预览框被一个带蓝色标签的拖动器所环绕。
图 6.80:启动Filter Sandbox。
- 在属性面板中,把
Filter端口设置为Median。你可以看到该滤波器已被应用于预览区域。 - 你可以拾取并拖动预览框,或调整其大小(先按
ESC键把3D查看器设为交互模式)。一开始保持较小的尺寸可以在尝试滤波器时节省时间。 - 你可以更改诸如迭代次数之类的滤波器参数。请注意,当把解释设置为3D时,该滤波器会被应用到一个3D板块(slab)上,其深度取决于参数。即使在有限的预览区域上,3D滤波也可能明显更慢。
- 你可以临时隐藏预览(
Preview端口),以便对比原图与滤波后的效果。 - 通过
Histograms端口,你可以显示在原始预览区域和滤波后预览区域上计算的直方图,从而查看该滤波器对更好地分离峰值有何帮助。你可以在直方图图形上右键单击,以在线性和对数刻度之间切换。属性面板中会显示一份统计摘要。 - 把
Preview type改为Interactive Thresholding:然后你就可以调整一个阈值,并验证该滤波器对这种阈值分割的影响。 - 按下
Apply按钮会把该滤波器应用到整个输入图像,并创建或更新一个结果图像。
图 6.81:在Interactive Thresholding模式下使用Filter Sandbox。
关于如何比较图像滤波器结果的信息,另请参阅教程第8.2节 数据融合、比较与合并数据,特别是其中演示如何同步视图与显示模块的示例。
6.7.3 组合图像滤波器
有时,依次应用两个或更多滤波器可以得到最好的结果。Slice模块允许用户选择性地把一系列滤波器(2D)应用到所显示的切片上:这提供了一种简便的方式来尝试单个滤波器或滤波器组合。大多数图像滤波器是不可交换的,因此你应当注意操作的顺序。
另一种增强图像的方式是把若干滤波器组合起来。你可以简单地把图像滤波器模块链式地附加到中间结果上。提示:勾选auto-refresh开关,结果会在输入或任何模块端口发生变化时随之更新。
在某些情况下,更复杂的组合可能会有用。有些滤波器以检测或保留边缘而闻名,另一些则用于平滑或降噪。这些滤波器可能相互矛盾,因而并不总是可以简单地依次应用。
下面你将看到如何逐步构建一个滤波器组合,使你能够:
- 保留边缘(有限的边缘平滑),
- 平滑那些本身均匀但有噪声的区域。
以下示例说明了滤波器组合的可能技术,但它并不是要展示某个特定实际场景下的首选方案,也不是一个通用解法。一般来说,像Bilateral或Non-Local Means这样的滤波器在平滑的同时保留边缘方面就已经能做得不错,仍应优先尝试。
整个过程分为几个步骤:
Strong Smoothing using Median Filter使用Median Filter进行强平滑Slight filter on the edges using Bilateral filter使用Bilateral filter对边缘做轻微滤波(需要Avizo许可证)Edge detection and masking using Sobel filter使用Sobel filter进行边缘检测与掩膜Compositing filters with mask用掩膜组合各滤波器
图 6.82:初始图像。
6.7.3.1 使用Median Filter进行强平滑
- 加载
AVIZO_ROOT/data/tutorials/image-processing-advanced/Catalyst.am。 - 把一个
Median Filter模块附加到Catalyst.am。 - 把
Interpretation设为3D并按Apply。 - 把一个
Ortho Slice模块附加到Catalyst.filtered。
加载项目data/tutorials/image-processing-advanced/FilterCompositing-1-MedianFilter.hx即可完成本教程步骤(见图6.83)。
图像均匀区域中的噪声已被降低,但边缘随之变得非常模糊。
图 6.83:中值滤波Median filter。
6.7.3.2 使用Bilateral Filter对边缘做轻微滤波
- 把一个
Bilateral Filter模块附加到Catalyst.am。 - 把
Interpretation设为3D,Kernel尺寸设为3,然后按Apply。 - 把一个
Ortho Slice模块附加到Catalyst2.filtered。
加载项目data/tutorials/image-processing-advanced/FilterCompositing-2-BilateralFilter.hx即可完成本教程步骤(见图6.84)。
边缘上的一些噪声已被去除,边缘也被清晰地保留了下来。不过,残余的噪声主要在图像的均匀区域上仍然可见。
图 6.84:双边滤波Bilateral filter。
6.7.3.3 使用Sobel Filter进行边缘检测与掩膜
- 把一个
Sobel Filter模块附加到Catalyst.am。 - 把
Filter设为3D并按Apply。 - 把一个
Ortho Slice模块附加到Catalyst-filtered。
加载项目data/tutorials/image-processing-advanced/FilterCompositing-3-SobelFilter.hx即可完成本教程步骤(见图6.85)。
该滤波器突出了边缘:边缘区域为白色,均匀区域为黑色,有噪声的区域为灰色。
图 6.85:Sobel滤波Sobel filter。
6.7.3.4 用掩膜组合各滤波器
现在,你将使用一个Arithmetic模块来合成滤波图像,其基本思路是:边缘区域使用双边滤波后的图像(B),均匀区域使用中值滤波后的图像(C)。这两个图像通过归一化的Sobel滤波图像(A)按如下方式线性混合:
其中(A)应当是一个取值范围为0到1的归一化浮点图像。
- 把一个
Convert Image Type模块附加到Catalyst-filtered.am。 - 由于
Catalyst-filtered.am的类型是范围为0-255的8位无符号字节图像,请把output type设为32-bit float,并把Scaling: scale设置为使输出范围变为0…1(使用0.00392157,几乎等于1/256)。你也可以改为在下面的表达式中把A写成A/256。 - 按
Apply以创建Catalyst-filtered.to-float。 - 把一个
Arithmetic模块附加到Catalyst-filtered.to-float。 - 把
Input B设为Catalyst.filtered(Median Filter的结果),把Input C设为Catalyst2.filtered(Bilateral Filter的结果)。 - 把
Expr设为B*A + C*(1-A)并按Apply。 - 把一个
Ortho Slice模块附加到Result。 - 把
Colormap范围设为0-255。
加载项目data/tutorials/image-processing-advanced/FilterCompositing-4-Arithmetic.hx即可完成本教程步骤(见图6.86)。
这种组合同时利用了中值滤波和双边滤波两者的优点:边缘被保留,均匀部分被平滑。你可以通过以下方式改进这类组合:
- 对初始图像进行滤波,以去除圆形的采集伪影,
- 应用明暗校正和/或直方图均衡化,
- 使用其他一些滤波器,
- 对掩膜进行滤波,
- 调整组合公式,
- 基于新的掩膜添加一系列算术运算。
图 6.86:带掩膜的滤波器组合Filters with mask。
6.8 关于标签测量的更多内容
6.8.1 测量组选择对话框
该对话框让你能够管理多个测量组。这些测量组可以被独立修改。它们会自动存储在用户设置中,因此修改在重启应用程序后依然保留。你可以从应用程序提供的测量列表中选择任意测量。关于每个测量的更多细节,可以在可用测量列表(list of available measures)中找到。也可以基于现有测量来创建自定义测量。
管理测量组Managing the measure groups
该对话框的上半部分显示所有测量组的列表以及一组工具。该列表可用于浏览测量组并选择某个组。
当从诸如Label Analysis这样的模块打开该对话框时,测量选择端口中当前所选的组会直接显示在对话框中。反过来,当该对话框被确认时,测量选择端口会被更新为对话框中所选的测量组。
利用测量组工具,你可以:
- 创建一个新的空组,
- 保存从项目文件或脚本加载的组(从GUI创建的组会自动保存),
- 把一个现有组复制为一个新组,
- 重命名一个现有组,
- 移除一个现有组。
以下默认组是不可编辑的。不过,每个默认组都可以用一个新名称复制出来,然后再编辑这个副本。
basic组,用于3D图像,包含测量Volume3d、Area3d、BaryCenterX/Y/Z和Mean。basic2D组,用于2D图像,包含测量Area、BaryCenterX/Y和Mean。Standard Shape Analysis组,包含标准形状参数测量。Weighted Shape Analysis组,包含加权形状参数测量。
图 6.87:测量组选择对话框The Measures Group Selection Dialog Box。
选择测量Selecting measures
左侧面板列出用户测量和原生测量。右侧面板列出当前组中已选中的测量。要把某个测量添加到选择组或将其移除,只需在列表中双击它即可。要一次添加或移除多个测量,请选中它们并使用中间的按钮 [>] 和 [<]。也可以使用 <Delete> 键快捷方式来移除所选的测量。
6.8.2 自定义测量
自定义测量可以通过把现有测量组合到一个数学公式中来创建。要创建一个新测量,请单击自定义测量列表上方的图标。在提示输入新测量名称之后,会打开测量编辑器Measure Editor。
创建完成后,新的自定义测量会列在自定义测量列表中。每个自定义测量旁边都有一些工具按钮,可用于:
- 重新编辑该测量的公式
- 把它从列表中移除
- 保存它——如果该自定义测量是从项目文件或脚本加载而来的(从GUI创建的自定义测量会自动保存)
测量编辑器The Measure editor
所有自定义测量都列在顶部的下拉菜单中。公式区域会根据所选测量而更新。所选测量可以用自定义测量列表旁边的工具按钮进行复制、重命名或删除。请注意,自定义测量在被删除之前,应先从所有测量组中显式地移除。
图 6.88:自定义测量编辑对话框The Custom measure Edition Dialog Box。
该对话框的主体部分用于编辑所选测量。第一个参数是该测量所产生结果值的单位量纲(unit dimension)。这个单位量纲将与被分析标签图像的单位相结合,以确定输出结果值的确切单位。例如,如果测量量纲是Area(面积),而输入图像的坐标单位是$cm$,那么输出值的单位就是$cm^2$。
第二个参数是该测量的数学公式。公式语法支持一组基本运算符和数学函数,可与现有测量配合使用。受支持的运算符列在公式区域的左侧。可用的函数和测量列在公式区域的下方。作为快捷方式,你可以双击底部列表中的关键字,或把它们拖放到公式区域中。
注意:公式中应使用工作单位(working unit)。关于工作单位如何设置的更多信息,请参阅自动确定或手动设置工作坐标单位。
公式定义会被实时检查,并根据其有效性着色:公式正确时为绿色,否则为红色。
6.8.3 可配置的原生测量
某些测量可以用额外的参数进行配置。这些测量在原生测量列表中以 […] 工具按钮标识。单击该工具按钮可打开标签测量属性编辑器Label Measures Attributes Editor。通过该编辑器可以编辑五类属性:
图 6.89:编辑Histogram与Feret测量的属性。
Feret angles,用于Feret 2D测量——这类测量在XY平面上沿每个单元的给定数量的直径方向进行测量。角度在[0,180]范围内均匀采样。默认情况下,有10个角度,即每18度一个。Feret 3D angles,用于Feret 3D测量——这类测量在每个单元周围的3D空间中进行测量。默认使用31个3D采样。Co-occurrence directions,用于基于共生矩阵计算的测量,按灰度级对像素对中给定方向(dx,dy)进行分类。共生矩阵的各分量由下式给出:
其中$I(x,y)$是坐标$(x,y)$处的图像灰度级。该等式意味着:对于给定的一对$(i,j)$,$M(i,j)$包含满足$I(x,y) = i$且$I(x+dx,y+dy) = j$的像素数量。该矩阵是对称的,并按如下方式归一化:
其中N是图像的灰度级数量。这些操作使其能够与图像尺寸无关,并保持关于某个方向及其对称性的性质。
Histogram parameters,用于基于每个标签的灰度直方图计算的测量。可配置的属性包括:所要考虑的灰度值范围,以及所生成的箱(bin)的大小。默认情况下,这些设置会被自动计算。Quantile values,用于可配置的HistoQuantile测量。另请参阅直方图的分位数(Quantile of an histogram)。Breadth 3D sampling,用于Breadth 3D测量——这类测量会在与主轴正交的方向上搜索最大的正交Feret直径,并与Feret 3D配合使用。Breadth 3D的值定义了在使用Feret 3D之后,在每个正交平面上所采用的采样。在使用该测量之前,应先把Feret 3D的采样设置为所需的值。
重要:属性值在各个测量之间是共用的。这意味着,例如当你为测量FeretShape编辑Feret角度的数量时,所有其他基于Feret角度的测量都会使用这个新值。
请注意,对话框中只有被所编辑测量支持的那些属性才会被启用。
6.8.4 关于测量组备份
当你保存项目时,重做该分析所需的全部信息都会被存储在项目文件中。这包括为该分析所选的测量组,以及自定义测量的公式。之后,当该项目在另一台计算机上被加载时,这些新的测量组定义和自定义测量也会被加载到新环境中。把这些新定义在新机器上进行本地存储是可选的,此时这些未保存的条目旁边会出现一个新的Save按钮。你不必保存这些新条目;它们是完全可用的,并且可以一直使用直到应用程序关闭。请注意,如果某个自定义测量或组在该机器上已经定义过但取值不同,那么来自项目的定义会被重命名,以避免任何冲突。
图 6.90:从项目加载的未保存测量Unsaved measures loaded from a Project。
6.8.5 脚本编写提示
要在脚本内部定义自定义测量,你可以使用labelMeasure create Tcl命令。该命令是项目文件中所用语法的一个简单替代方案。项目文件确实使用了一套更复杂的机制,以避免从某个项目加载多个自定义测量时出现冲突和重复,但那套机制并不适合用于脚本编写。
labelMeasure create Tcl命令以一种更简单的方式处理冲突:如果该测量名称已被使用,就会生成一个新名称,并把这个名称作为命令的结果返回。返回的名称应当在依赖它的测量公式中使用。使用硬编码的名称可能会把依赖关系替换成一个错误的对象。示例:
1 | labelMeasure create mymeasure length "Volume/Area" |
这在控制台中是可行的,但如果用在脚本中、而脚本环境里已经存在一个名为”mymeasure”的测量时,就可能引发隐藏的依赖问题。在脚本中,请优先使用:
1 | set msr [ labelMeasure create mymeasure length "Volume/Area" ] |
6.9 使用环境光遮蔽进行空腔分析
以下教程需要Avizo许可证(Avizo Lite Edition许可证是不够的)。
本教程说明如何使用Compute Ambient Occlusion模块,按照文献[1, 2]中所描述的方法在二值图像中提取空腔。这里,空腔被假定为物体内部的中空空间。通常,空腔属于背景,而物体是前景。尽管在本教程中我们关注的是空腔的分割,但该概念同样可以应用于孔隙结构。
如果空腔与物体外部的背景相连通,Compute Ambient Occlusion模块会特别有用。该模块利用了环境光遮蔽(ambient occlusion)的思想:对于埋在物体内部的区域,环境光遮蔽值很高;对于物体外部的区域,环境光遮蔽值很低。环境光遮蔽的计算方式是:从每个背景体素向各个方向投射射线,射线的数量是一个用户自定义的值。如果某条射线击中了前景,就意味着投射该射线的那个体素在这个方向上被遮蔽了。击中前景的射线数量与所投射射线总数之比,就给出了环境光遮蔽的值。
在本教程中,我们将为一个石头数据集(数据由英国的Richard Watson提供)计算空腔,该数据集带有需要被分割出来的生物侵蚀痕迹。展示本教程结果的一个示例可以在此处找到。
6.9.1 前提条件
6.9.1.1 输入图像的二值化
该算法要求以一个二值图像作为输入。因此,第一步我们需要对图像进行二值分割,以区分前景和背景。灰度图像可以通过阈值处理进行二值化。另请参阅高级图像处理与定量分析入门中的”灰度图像的二值化”一节。我们通常使用以下方法之一:
Auto Thresholding- 如
Watershed Segmentation向导中所描述的分水岭分割 Segmentation Workroom分割工作室中的阈值工具
要学习接下来的步骤,请从data/tutorials/cavityanalysis/文件夹加载数据集Stone_binary-segmentation.am。
6.9.1.2 硬件
Compute Ambient Occlusion模块需要一块支持CUDA计算能力2.0或更高版本的NVIDIA显卡。为防止出现CUDA错误和超时,强烈建议使用高端显卡。
6.9.2 工作流
从环境光遮蔽场计算出阈值分割结果,所建议的工作流分为以下几个步骤:
Computation of the ambient occlusion field环境光遮蔽场的计算Cavity segmentation空腔分割Postprocessing后处理
6.9.2.1 环境光遮蔽场的计算
为了计算环境光遮蔽场,请执行以下步骤。
- 用鼠标右键单击对象池中的
Stone_binary-segmentation.am对象,并在弹出对话框的搜索行中键入单词Ambient。选择Compute Ambient Occlusion条目。结果,对象池中会出现一个新的红色对象(见图6.91)。 Max Distance端口的默认值被设置为其最大值。这确保了所投射的射线在到达某个前景体素或边界框之前不会终止。为了加快计算速度,应当减小该值。- 把
Number of Rays端口设置为100以提高采样,不过通常50就已足够(见图6.92)。为了计算该二值图像背景的环境光遮蔽场,请按下Apply按钮并等待。根据你的显卡情况,计算可能需要几秒到一分钟。
图 6.91:池中环境光遮蔽场的计算工作流。
图 6.92:环境光遮蔽属性。
结果,对象池中会出现一个名为Stone_binary-segmentation.ambientOcclusion的新数据对象(见图6.93)。这是一个浮点标量场,包含-0.1到1.0之间的值。-0.1这个值是保留给前景体素的。对于背景体素,取值在0.0到1.0之间变化,其中1.0表示完全被遮蔽,0.0表示完全没有被遮蔽。
图 6.93:环境光遮蔽结果。
6.9.2.2 空腔分割
为了以环境光遮蔽场作为输入来计算阈值分割,请执行以下步骤。
- 在对象池中右键单击该环境光遮蔽场,并为它附加一个
Interactive Thresholding模块。 - 转到
Intensity Range端口,把下限阈值设置为0.7。我们的经验表明,对于许多数据集而言,这都是一个适合从环境光遮蔽场计算阈值分割的范围。 - 按下
Apply按钮。
这样就生成了空腔结构的二值分割结果,例如可以使用Voxelized Rendering模块将其可视化(见图6.94)。
图 6.94:空腔分割结果Cavity Segmentation Result。
6.9.2.3 后处理
在大多数情况下,该分割的结果仍会包含一些噪声,也就是一些互不相连的小标签。你可以在Segmentation Workroom分割工作室中进行后处理。一种做法是使用魔术棒Magic Wand工具(见图6.95)。更多细节请参阅分割工具(Segmentation Tools)一节。
图 6.95:分割编辑器Segmentation Editor。
参考文献References
- D. Baum, J. Titschack, Cavity and Pore Segmentation in 3D Images with Ambient Occlusion, EuroVis 2016 - Short Papers, 2016, DOI: dx.doi.org/10.2312/eurovisshort.20161171
- J. Titschack, D. Baum, K. Matsuyama, K. Boos, C. Färber, W.-A. Kahl, K. Ehrig, D. Meinel, C. Soriano, S. R. Stock, Ambient occlusion - a powerful algorithm to segment shell and skeletal intrapores in computed tomography data, Computers and Geosciences, Vol.115, pp. 75-87, 2018, DOI: dx.doi.org/10.1016/j.cageo.2018.03.007
