独立设计、实现的一个三维医学图像分割软件，能够在体绘制视图中直接操作、分割三维数据，代码量10000行。

基于交互的体数据分割系统设计

1. 选题背景

体数据(Volume data)是一种常见于医学图像等领域的数据类型，具有三个维度。

如同二维图像的分割，体数据的分割也是其常用的操作。然而，增加的一个维度导致了数据量的大幅增加，进而增大了分割的难度；与此同时，三维的数据无法直接在二维的屏幕上显示，又进一步增加了难度。

现有的体数据分割方法可分为自动、半自动、人工分割三种。其中，自动分割方法有很大局限性，经常难以分割出想要的结果；相比之下，人工分割更加可靠，但需要一定的专业知识，也十分耗时。

上图为使用Seg3D软件进行人工标注分割的操作界面。图中所展示的是196个切片中的第112个，其中橙色为肿瘤，浅蓝色为健康组织。使用画刷工具为各个部位刷上颜色便完成了一个切片的分割，对于全部196张切片皆需要进行这样的操作，因而工作量极大。

我选择的研究内容是优化人工分割的交互方式，在不降低分割精度的前提下尽可能加速分割的过程。

2. 国内外研究现状

2.1 自动方法

2.1.1 基于水平集的分割

水平集方法的原理是给定一个初始轮廓，根据图像的亮度、梯度等信息，不断使该轮廓演化，最终逼近物体的真实边界。

传统的水平集方法对噪声敏感，轮廓容易停止于噪点处。Yang等人[1]在2015年提出的嵌入MRF的水平集方法改进了这一问题，对于医学图像分割也有较好的效果。

水平集方法适用于二维图像的分割，对于三维图像则需要逐层分割。因为未考虑层间的关系，这一方法对三维图像并不十分适用。

2.1.2 基于atlas的分割

基于atlas的分割是常用的医学图像分割方法，其中atlas可以理解为分割的模板。这类方法是在三维空间上进行的，且使用了已分割好的数据作为先验知识，因而能够得到较好的结果。

上图解释了基于单个atlas的分割方法的工作原理。首先选定一个参考图像作为atlas，这个参考图像需要人工地进行标注。然后输入待分割图像(目标图像)，将两图像进行配准，得到一个配准映射。对atlas的分割应用配准映射得到的新的分割即可作为目标图像的分割。

[2-5]是近年来基于atlas的分割方法中效果较好的几篇文章，分别从不同的环节对传统方法做出了优化。

这类方法的不足有两点:

• 耗时长。为了分割精确通常需要使用多个atlas，这使得分割时间成倍增加，分割时间通常长于人工分割时间。
• 可靠性差。由于分割完全依赖于配准映射，分割结果与预期常有或多或少的出入。

2.2 人工交互方法

体数据分割相比于二维图像分割的难点在于增加的一个维度。对于二维图像，用户可以使用简单的操作方法，如阈值、Floodfill(photoshop中的魔棒工具)方法来分割。但是对于体数据，

2.2.1 Seg3D软件

Seg3D是一个免费软件，可以用于体数据的分割。软件中提供了阈值、画刷、连通域、布尔运算、图层、滤波等功能，可以辅助人工分割。但是，所有的分割操作均需要在体数据的三个方向的截面(横断面、矢状面、正中面)进行，均是二维的操作。对于体数据分割，没有三维的操作显然不能满足需求。

2.2.2 Teddy: A Sketching Interface for 3D Freeform Design[6]

Teddy是Igarashi等人于1999年开发的基于sketch的三维物体建模工具，截至目前已有超过1600次的引用量。论文中将二维操作转化为三维操作的思想十分有借鉴意义。

2.2.3 WYSIWYG Volume Visualization[7]

WYSIWYG(what you see is what you get)体可视化是北京大学可视化实验室提出的方法。

该方法实现了通过交互来为体数据上色，并体绘制的技术。

2.2.4 Interactive Volume Segmentation with Threshold Field Painting[8]

这是本次作业主要参考的文章。

这篇文章同样是Igarashi等人发表的，发表时间是2016年10月。2.2.3中介绍的WYSIWYG方法针对的是体数据，这篇文章与之类似，但处理的是等值面的数据。

文章的思路非常简单，指出使用阈值法做体数据分割面临的一个大问题是不存在能够恰好分割数据的合适的阈值，于是提出在不同空间位置应用不同的阈值进行分割。文章的主要贡献是设计了构建这样的阈值场的方法。

3. 实现过程

3.1软件结构设计

为了日后与实验室项目整合，此次代码采用了与实验室项目相同的环境，即:

• Visual Studio 2012
• Qt 4.8.0
• VTK 7.0.0

软件主界面由工具区和显示区构成。工具区包含了各种操作按钮。显示区包含了若干个查看器(Viewer类)，可以方便地设置各种布局的组合。

3.2 体数据读取

医学体数据通常是用DICOM格式(.dcm)存储的，每个切片存成一个文件，读入时读取一个文件夹。DICOM标准非常混乱，不易操作，而且许多开源软件也不支持。因此，我选择处理应用范围更广的NIfTI格式(.nii, .nii.gz)。

NIfTI网站提供了开源的工具包，编译好即可使用其中的读取函数。从NIfTI图像中可读到图像的尺寸、层间距信息，图像以short一维数组形式存储，数据范围大约为0~5000。

3.3 切片显示

程序的一个功能就是查看图像的切片。每显示一张切片都需要从体数据中将特定坐标的体素数值遍历一遍，但实现发现对于不太大的图像(小于$512\times 512$)，速度尚可。

一开始我使用的是OpenGL，实现了约两天时间(下图)，debug非常困难，因此改用封装好的三维引擎来显示。先是选择了Ogre，在迁移过程中为了解决各种问题耗费了近一周的时间。最终还是选择使用VTK来显示，迁移又花了一些时间。

三维数据相比于二维数据具有更多的朝向，因此处理起来也比二维数据难得多。为了正确摆放切片，我将数据中坐标小于50的体素全部变成黑色，便可方便地确定其方向，如下图的三维视图中，两个黑色区域没有对齐，便可知道有一个切片的方向放错了。

三维数据的切片间距处理也是一个难点，需要对显示的物体的scale做相应调整，图像在三个维度上的spacing通常是不同的。

3.4 三维视图的旋转、平移、缩放

查看三维的医学数据时，通常不希望出现数据的上下方向发生倾斜，例如进行旋转操作时头顶方向朝向了右侧。这就需要对VTK默认的Trackball Controller进行改进。Three.js中的Orbit Controller恰好符合这样的要求。

于是我利用VTK的vtkCamera类实现了Orbit Controller，可用鼠标、键盘操作模型，同时保证相机的竖直方向在世界坐标系中也始终保持竖直。

鼠标的交互还考虑到了视图的缩放问题，经过调整参数，进行平移和旋转操作时，鼠标指针与模型间的相对位移很小。

3.5 等值面生成

论文中的显示及交互都是基于等值面的，VTK中也提供了用于生成等值面的Marching Cubes算法:

然而，该方法效率过低，对于上图中$200^3$尺寸的数据都要生成近一分钟的时间，而多数医学图像比这个还要大。导师说，Marching Cubes算法可以使用GPU加速，然后我在GitHub上找到了一个使用OpenCL实现的Marching Cubes代码。

但是将该代码迁移到工程中并不容易。代码最后一次更新是五年前，在此期间OpenCL甚至都发生了变化。代码编译始终不通过，最终发现是作者使用了cl::vector而不是std::vector，改成后者则解决了问题。

但是他的程序运行一段时间就会崩。OpenCL的部分无法debug，后来使用GPU-Z查看发现程序占用的显存会越来越大，到了极限(4G)就崩了。因此考虑应该是代码在重新生成等值面时没有清除显存，手动清除即解决了问题。修正之后，等值面生成速度非常快，$512\times 512\times 381$的数据也能达到实时。

他的代码的渲染使用的是OpenGL，而我需要使用VTK来显示。他的代码实现非常粗糙，我首先重构了代码，整理成了面向对象的。Marching Cubes算法生成的等值面是以顶点数组的形式存储在显存中的，每个顶点占用6个float，表示空间坐标及法向量，使用OpenGL可以用glDrawArrays可以直接绘制。但是VTK应该是没有办法直接绘制显存中的顶点数组，于是我将数组复制到了内存中，这个步骤可以在瞬间完成(尽管顶点数可多达千万)。然而VTK初始化网格模型的过程很慢，比直接用OpenGL显示慢了数十倍。

整个等值面生成工作一共花了两周时间。

3.6 体绘制

整合好等值面生成的代码后，需要对不同区域着上不同颜色，对于数据量动辄多达千万的体数据这是非常困难的。

然后突然意识到并不需要使用等值面，用体绘制完全可以实现论文的方法。

3.6.1 模型显示

体绘制通常是用来高质量地渲染三维数据的。因为医学图像中，不同组织通常有特定的亮度值，使用预设的传递函数将亮度值映射到特定的颜色、透明度便实现了逼真的渲染效果。

为了用体绘制模拟等值面的显示，我构造了这样的传递函数:

$$f_{opacity}(I)=\begin{equation} \begin{cases} 1, I\ge threshold\\0, \rm otherwise\end{cases}\end{equation}$$

它使体绘制只显示大于阈值的部分。

体绘制是使用显卡计算的，所以速度很快，可以达到实时绘制。

3.6.2 在体绘制中显示不同颜色

3.6.2.1 初步尝试

VTK的体绘制不支持为有标记数据显示多种颜色，因此可构造一个特殊的传递函数来实现这一功能:

$$f_{opacity}(I)=\begin{equation} \begin{cases} 1, I\in [T(A),T(A)+\delta]\cup [T(B),T(B)+\delta]\\0, \rm otherwise\end{cases}\end{equation}$$

$$f_{color}(I)=\begin{equation} \begin{cases} Color(A), I\in (-\infty, T(B))\\Color(B),I\in [T(B),+\infty)\end{cases}\end{equation}$$

其中，$T(A)$、$T(B)$表示映射到两个颜色的阈值，$\delta$表示区间的宽度(只需要保证两区域不重叠即可)。如下图，亮度处于区域A中的体素会呈黄色，亮度处于区域B中的体素会呈蓝色。在本项目中，数据一开始都位于黄色区间，此时低于阈值的区域的opacity为0，表现为透明; 高于阈值的区域为不透明区域(见3.6.1节中的绘制效果)。

接下来，将分割区域的全部体素加上10000(一个足够大的数)，让其映射到区域B，那这部分体素便(理应)显示出蓝色。

但实际的绘制结果并不是这样。在蓝色的外部，与透明区域交界处，出现了一层金色的膜。

分析得知，这是体绘制的插值引起的。

体数据是离散表示的，各个体素之间存在间距，为了在体绘制中显示更精细的结果，通常要在两个相邻的体素间进行插值。对于通常的、自然的数据，这样做是没有问题的，但是对于我上面使用的特殊的传递函数，问题便出现了。

现在要显示的是区域B，它与外面的透明区域直接相连。根据上面插值的原理，"空气"与区域B中间的数值将会穿过区域A(见上图)，而区域A是不透明的，这导致了上述薄膜的出现。因此应该改变方法。

3.6.2.2 尝试修改传递函数

经过思考，我想到可以利用负数区域。现有数据均是大于0的，所以，将区域B移动到比数据最小值还小的位置，从空气到区域B之间便不会穿过区域A，也就没有先前的虚面了:

注意金色区域与银色区域中间有缝隙，这也是插值导致的，因为插值经过了一个透明区域。

现在新的问题来了，标记区域的表面出现了"台阶"，不像先前那样光滑(见上图)。

3.6.2.3 去除台阶效果

上图中，金色区域是标记出的体素，我将它们全变成了负值，便映射成了金色。但是它们周围，没有被标记为分割区域的部分是没有改变的，仍是较大的正值。未经修改的地方看起来光滑是因为数据是连续变化的，这样使用插值即可得到体素间本不存在的数据点; 而金色区域与空气间的差非常大，插值不起作用，所以出现了台阶效果。

解决方法是将分割区域连同其外部附近的区域都变为负值，而负数区域的阈值保持先前的阈值，即得到了光滑的表面(下图)。

3.6.2.4 显示更多的颜色

上述方法虽然能完美地显示两种颜色，但如果想显示三种或者更多的颜色则是做不到的。如果设置三段的传递函数则会出现下图的情况:

这里同样出现了插值导致的虚面，这是从空气(0左右)到蓝色区域(-20000左右)中间穿过了黄色区域(-10000左右)时发生的。

从原理上考虑，该现象无法避免。所幸，VTK提供了多种插值模式(nearest、linear、cubic)，选择近邻插值虽然牺牲了显示效果，却避免了上面的问题(下图)。这样，当只有金色、银色的时候，使用立方插值，有蓝色的时候使用近邻插值，就可以了。

3.6.3 体绘制的拾取

VTK中包含的vtkCellPicker类可直接对体绘制进行拾取，方法是从视点发出一条射线，在穿过体数据的时候逐步积累透明度，当透明度超过阈值便判定当前点为拾取点。

不幸的是VTK的拾取使用的是线性插值，所以会错误地拾取到上述的虚球面上。

解决方法是继承vtkCellPicker类(在代码中记为VolumePicker类)，重写其中的透明度累积函数，改成近邻差值便可正确拾取。

3.7 数据标记工具

软件提供了多种数据标记工具。使用这些工具将首先在体绘制结果中标记出蓝色区域，然后用户可以选择接下来的操作(按钮位于界面左侧，均为红色):

3.7.1 球形标记工具

使用球形标记工具时，鼠标右键点击体绘制模型可进行拾取，在拾取区域将出现一个球体。此时按下Enter可以标记该区域，按下Delete可取消蓝色标记。也可以在按住Shift的同时右键拖动鼠标来连续标记，或按住Ctrl来连续取消标记。使用鼠标滚轮可以改变拾取球的大小。

该工具的实现只需要找出以拾取坐标为中心的拾取球内所包含的体素，但由于图像具有spacing，图像空间坐标与模型所在的世界坐标是不一样的。为了实现这一功能，需要先求出拾取球的包围盒在世界坐标系下的坐标，将包围盒转换到图像空间坐标，使用三层循环遍历之，对于每个体素，将体素的坐标(图像空间)再转换回世界坐标，判断与圆心的距离即可。

若打开PU按钮(pickup slices)，则在三维视图中的拾取可在二维视图中同步更新; 同样地，在二维视图中点击也可改变三维中拾取球的位置。用户也可以在二维视图里进行标记操作。

对球形标记工具的一个优化是加入了拾取物体内部的方法。拾取时，ray casting的结果是物体的表面，若在此处生成一个球体，那么球会有大约一半的体积在物体外部，而我们希望的是球体尽可能将物体包裹起来，所以需要将球体向物体内部移动。我将拾取球沿着拾取点的法向量的反方向移动$\frac{R}{1.5}$($R$是拾取球半径)，这样做的好处是，当我想使拾取球的球心位于一个球形物体的球心时，只需要调整拾取球的半径为球形物体半径的1.5倍，再点击物体表面即可。

3.7.2 floodfill

球形标记工具适用于简单的场景，但是对于复杂的情况，例如血管等，使用该工具是较为费时的。

floodfill工具可以快速地将当前拾取点的连通域全部标记为蓝色，连通与否是根据当前的全局阈值确定的。

该工具也使用了球形工具中拾取物体内部的优化，否则用户点击的将是物体表面，得到的结果可能不对。

3.7.3 region growing

与floodfill工具类似，region growing工具也要先拾取坐标点。该工具会找出与起始点连通、亮度差异小于一个阈值，且梯度小于另一个阈值的所有体素。

该方法可以较好地分割血管。

3.7.4 dilation & erosion

使用floodfill、region growing方法分割的结果是具有"突变"边界的，也就会产生前文所述的"台阶"现象，如上图，金黄色的血管表面较为粗糙。前面已经提到，这是因为分割结果与外部存在亮度的跳变，为了使结果看起来更加平滑，可以多分割一些区域，将分割物体的边缘体素包含进来。

实现这一过程可以使用形态学处理中的腐蚀、膨胀操作。类似于二维图像的腐蚀、膨胀，我使用了这样的模板:

4. 实现结果及分析

4.1 程序的全部功能

4.1.1 查看器

查看器是用于显示三维或二维模型、图像的控件。

程序中包含两种查看布局，分别是"左1右3"(默认)、"左2右2"的形式。点击程序左下方layout按钮可以切换两种布局模式。

各个查看器是独立的，且可以通过鼠标拖动调整大小。在查看器左下角的下拉菜单中可以选择当前查看的内容，可选项包括"Volume"(三维数据)、"Sagittal"(矢状面)、"Coronal"(冠状面)、"Axial"(横截面)。

查看器左下方显示了一些文字信息，包括:

• (3D) 当前拾取点在图像空间的坐标，以及此处图像亮度值。
• (3D) 当前等值面阈值。
• (3D) 当前窗位。
• (3D) 当前窗宽。
• (2D) 当前切片所在位置。

被选中的查看器的边框会变为橙色，键盘操作仅在此时有效。

二维查看器的相关操作有:

• 鼠标左键点击或拖动: 改变十字线所指位置，该操作会更新除当前查看器外其他两个查看方向的切片位置。
• 鼠标滚轮滚动: 切换当前查看器的切片位置，每滚动一下切换三张切片; 如果按住Ctrl则切换一张(减速); 如果按住Shift则切换七张(加速)。
• 鼠标右键拖动: 改变窗位(window level)、窗宽(window center)，水平方向是窗位，垂直方向是窗宽。按住Ctrl可减速; 按住Shift可加速。因为图像数据的范围约为0~5000，而显示器只能显示256中不同的灰度，所以需要做一个映射: $$Value=\frac{Intensity(x)-(level-width/2)}{width}\times 255$$。

二维查看器下方的按钮有:

• CR(Cross lines): 显示/隐藏十字线。
• TR(Transparency): 启用/禁用透明模式。若启用，则图像的亮度通道也将赋值给透明度通道，表现为黑色区域变得透明、越亮的区域越不透明。由于VTK的deep peeling方法实现得不好，该功能存在一些问题。
• TS(Threshold): 启用/禁用分割预览。若启用，则切片图像会稍变暗，同时，当前阈值下的显示部分将显示绿色(对应于三维视图中的银色部分)，三维视图中的金色、蓝色部分在二维视图中也对应地显示为黄色、蓝色。

三维查看器的相关操作有(以下操作均在View模式下):

• 鼠标左键拖动、W/S/A/D键: 旋转模型。这里使用了Orbit Controller，旋转时竖直方向不会发生倾斜，且上下旋转是有边界限制的。按住Ctrl可减速; 按住Shift可加速。
• 鼠标滚轮滚动、E/Q键: 放大/缩小模型。按住Ctrl可减速; 按住Shift可加速。
• 鼠标右键拖动、↑↓←→键: 移动模型。按住Ctrl可减速; 按住Shift可加速。
• 按住Alt键: 暂时切换当前操作模式。例如，现在是View模式，按住Alt则切换到Operation模式。
• 空格键: 切换操作模式。
• R键: 重置相机位置。

三维查看器下方的按钮有:

• SL(Slice): 显示/隐藏切片。若显示，则在三维视图中会将二维视图的切片显示在其真实的位置处。
• VR(Volume Rendering): 显示/隐藏体绘制。
• MC(Marching Cubes): 显示/隐藏等值面。注意，仅金色部分会被显示。
• PU(Pickup slices): 同步二维、三维查看器的切片位置。在三维视图中的拾取也可在二维视图中同步更新; 同样地，在二维视图中点击也可改变三维中拾取球的位置。用户也可以在二维视图里进行标记操作。
• SG(Segmentation): 单独显示/取消分割部分。若启用则界面中会只显示蓝色、金色区域。

4.1.2 分割功能

4.1.2.1 概述

程序进入时默认的状态是View模式，在各个查看器中按下空格，或者点击左侧面板中的operation按钮则可以切换到Operation模式。仅在此模式下，可以修改图像数据(例如进行分割操作)。

在三维查看器体绘制中，

• 银色区域表示图像中亮度值大于当前全局阈值的部分;
• 金色区域表示被分割出的部分;
• 蓝色区域表示当前选中部分。

对于被选中部分(蓝色)，可以进一步地进行下列操作:

4.1.2.2 操作方法

与分割相关的左侧面板按钮有:

• sph(sphere): 球形标记工具。
• fld(floodfill): floodfill工具。
• rgn(region growing): 区域生长工具。
• lck(lock): 标记区域锁定工具。在启用状态下，floodfill与区域生长工具的作用范围将被局限在当前蓝色区域内，即，遇到非蓝色区域将会使两个算法停止生长。使用两种工具标记后，将取消当前全部蓝色标记，并用新的蓝色标记更新。
• dil(dilate): 膨胀工具。点击此按钮将直接对蓝色区域运行一次膨胀运算。
• ero(erode): 腐蚀工具。点击此按钮将直接对蓝色区域运行一次腐蚀运算。

相关操作有:

• 鼠标左键拖动: 旋转模型。按住Ctrl可减速; 按住Shift可加速。
• 鼠标滚轮滚动: 改变拾取球的大小; 若同时按住Shift则会改变当前的全局阈值; 若同时按住Shift和Ctrl则会缓慢改变全局阈值。
• 鼠标右键点击或拖动: 改变拾取球位置。
• 软件下方的滑动条: 改变全局阈值。
• E/Q键: 改变拾取球大小。
• R键: 将相机关注点移动至拾取球球心。
• ←→键: 改变全局阈值。
• sph工具下:
  • Shift+鼠标右键拖动: 将标记球经过的区域标记为蓝色;
  • Ctrl+鼠标右键拖动: 取消标记球经过的区域的蓝色标记;
  • Enter: 将当前标记球内区域标记为蓝色;
  • Delete: 取消当前标记球内区域的蓝色标记。
• fld或rgb工具下:
  • Enter: 从当前标记球球心开始运行floodfill或者region growing算法，算法算出的区域会被标记为蓝色;
  • Delete: 从当前标记球球心开始运行floodfill或者region growing算法，算法算出的区域会被取消蓝色标记。

除此此外，还有4.1.2.1节中提到的四种操作方式。

4.1.3 数据读取及保存

左侧按钮面板上方的import、save按钮对应于数据的读取及保存功能。

软件支持NIfTI格式图像的读取，其扩展名为.nii或.nii.gz。

进行完分割操作后，可保存当前分割的图像为.nii格式，其中存储了分割的状态信息。也可以选择导出网格文件(.obj格式)，该操作会生成金色分割区域的等值面并保存，保存结果可用其他网格查看软件打开。

4.2 结果分析

测试数据包括CT与核磁图像，均可正常打开、进行各种交互操作。

对于脑部核磁图像，本软件无法较好地进行分割，这是因为脑核磁分割的目标通常是脑的某一功能区，例如海马体。在核磁中，海马体与周围的灰度差异不明显，较难使用阈值法分割。

对于足部的CT图像，软件可以很容易地分割出其中的某一块骨骼。分割过程为:

• 调整合适阈值，观察三维、二维查看器，使骨骼间完全分离;
• 使用floodfill工具，拾取待分割骨骼上的点，观察二维查看器，确保拾取点在骨骼内部;
• 运行floodfill，观察运行结果;
• 若floodfill结果中包含其他骨骼，则使用球形标记工具，切除骨骼间的关联，重新运行floodfill;
• 使用膨胀工具，使分割边界变得平滑;
• 调低阈值，使骨骼中间未被分割区域在阈值以上;
• 使用球形标记工具将中间未被分割区域填充上;
• 点击左侧面板select按钮，确认分割;
• 点击MC按钮，生成等值面。

对于头颈部CT图像，软件可以分割出其中的血管。粗血管的分割过程为:

• 调整合适阈值，使血管可以清晰地显示;
• 使用区域生长工具，选择一条较粗的血管，调整拾取球半径为血管半径1.5倍左右，点击血管表面;
• 运行区域生长，观察运行结果;
• 使用若干次膨胀工具，使分割边界变得平滑;
• 点击左侧面板select按钮，确认分割。

细血管的分割过程为:

• 选取球形标记工具，右键点击血管，调整拾取球大小;
• 按住Shift，使用鼠标右键沿着血管拖动，标记血管;

• 若分割区域包含其他不想要的部分，则:
  • 若不想要的部分与血管连通，则使用球形工具拾取连通点，按Delete键去除关联;
  • 否则，使用floodfill工具，拾取血管，运行floodfill;
• 点击左侧面板select按钮，确认分割。

脑内部血管观察方法:

• 选取球形标记工具，右键点击模型出现拾取球，调整拾取球大小大约为脑内部空间大小;
• 在二维视图上选择脑内中心位置;
• 按Enter键标记球内区域;
• 点击三维查看器下方的SG按钮，观察分割结果;
• 点击左侧面板select按钮，确认分割。

5. 总结与展望

5.1 总结

参考论文[8]的体数据分割效果非常好，但其中展示的系统仅是一个demo，选取的图像尺寸很小，且除分割功能外没有其他功能。

本项目实现了论文中大部分功能，以及下列创新点:

• 改用体绘制技术来更高效地显示更大规模的数据，支持显示多种颜色标记，且渲染效果比论文中使用的等值面更光滑;
• 添加区域锁定工具，限制了floodfill与region growing的范围，该工具可以与球形工具结合使用，增强了分割功能;
• 实现Orbit Controller的三维旋转，限制旋转方向;
• 添加腐蚀、膨胀工具;
• 添加二维图像查看器，可显示原始图像及彩色的分割标记;
• 三维、二维的拾取可以联动进行，方便用户查看当前分割的位置。

整个工程的代码共约10000行(其中有500行的开源代码 )，85个代码文件，实现时间为40天。

5.2 展望(不足)

在人工分割工作中，论文[8]中提到的圈选去除工具(圈出屏幕一个区域并删除其中的物体)是非常重要的工具，实现该工具需要计算鼠标圈出的切割体与体数据的交集。但由于时间有限，我未能将其实现。

此外，无论是原论文，还是我实现的工程，都无法很好地处理藏于物体内部的结构，将来的版本中，可实现为将体绘制暂时地"剖开"(裁剪)进行处理，VTK提供了这样的函数。

各种分割操作都需要对图像数据进行赋值操作，若数据量大则会有卡顿感，影响了交互体验。有些操作可以尝试使用GPU加速。

受GUI限制，当前的腐蚀膨胀操作半径都是1体素，执行多次操作之后区域边界将会出现棱角。

目前VTK不支持有标记数据的体绘制，所以我使用了替代方法来为体绘制结果上色。而如果VTK原生地支持彩色图像(每个体素均有颜色值而不单是亮度)便可更好地渲染。

等值面生成使用的是简单的Marching Cubes算法，显示效果比体绘制差，特别是对于较细的血管，该算法生成的等值面质量很低。对此可以尝试使用插值的方法，先对数据进行超采样，再生成等值面。

论文方法与本软件的适用范围有限，较适合骨骼、血管的分割，而不适合诸如脑功能区、肿瘤等阈值区分不明显的区域的分割，且不支持阈值上界(当前只有下界)，仍待改进。

参考文献

[1] Yang, Xi, et al. "An efficient MRF embedded level set method for image segmentation." IEEE transactions on image processing 24.1 (2015): 9-21.

[2] Bai, Wenjia, et al. "A probabilistic patch-based label fusion model for multi-atlas segmentation with registration refinement: application to cardiac MR images." IEEE transactions on medical imaging 32.7 (2013): 1302-1315.

[3] Wang, Hongzhi, et al. "Regression-based label fusion for multi-atlas segmentation." Computer Vision and Pattern Recognition (CVPR), 2011 IEEE Conference on. IEEE, 2011.

[4] Wang, Hongzhi, et al. "Multi-atlas segmentation with joint label fusion." IEEE transactions on pattern analysis and machine intelligence 35.3 (2013): 611-623.

[5] Sanroma, Gerard, et al. "Learning-based atlas selection for multiple-atlas segmentation." Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. 2014.

[6] Igarashi, Takeo, Satoshi Matsuoka, and Hidehiko Tanaka. "Teddy: a sketching interface for 3D freeform design." Acm siggraph 2007 courses. ACM, 2007.

[7] Guo, Hanqi, Ningyu Mao, and Xiaoru Yuan. "Wysiwyg (what you see is what you get) volume visualization." IEEE Transactions on Visualization and Computer Graphics 17.12 (2011): 2106-2114.

[8] Igarashi, Takeo, et al. "Interactive Volume Segmentation with Threshold Field Painting." Proceedings of the 29th Annual Symposium on User Interface Software and Technology. ACM, 2016.

实现了二维图像的配准工具，支持平移、旋转、仿射、透视变换，可使用暴力法、梯度下降法、模拟退火法优化求解。

GitHub Pagelink

1. 概述

本次实验中，我基于OpenCV，实现了一个二维图像配准工具，全部代码均为自行实现，OpenCV用于计算图像变换与相似度。

该工具能够将一幅图像进行变换，并与另一幅图像相匹配。支持包括平移、旋转(含平移、缩放)、仿射与透视共四种变换，使用L1、L2、无穷范数作为优化的目标函数，实现了暴力算法、梯度下降法、模拟退火算法来求解该优化问题。

2. 应用问题

如果两幅图像，它们是在同一场景、不同角度下拍摄的，那么，存在一种图像变换，使得其中一幅图像经过变换后，能与另一图像大部分重合。

上述图像变换被称为配准(registration)$T$，被变换的图像被称为参考图像$I_M$，另一图像被称为目标图像$I_F$。优化的目标是使变换后的参考图像$T(I_M)$与目标图像$I_F$的差异尽可能低。

最简单的图像变换是平移变换，需要确定两个参数: $\Delta x$和$\Delta y$; 旋转变换通常与缩放、平移共同进行，需要确定四个参数: $\Delta x$、$\Delta y$、$\theta$、$scale$; 仿射变换将矩形图像线性映射至一个平行四边形，需要确定三个坐标点，共六个参数，三个坐标点分别表示原图左上、右上、左下角变换至新图像的坐标位置; 透视变换与仿射变换相似，不同的是原图像的四个顶点可变换至任意的四边形，所以需要确定四个坐标点，共八个参数。此外，也有更为精细的图像变换方法，但相比于上述简单变换，其参数较多，难以优化，故本次实验不予考虑。

对于图像相似度，需针对使用场景选择合适的度量方法。本实验中，实现的方法有L1(1范数)、L2(2范数)、无穷范数三种。

总的来说，问题可以总结为如下步骤:

• 输入参考图像、目标图像;
• 选择合适的变换，确定参数范围;
• 设置初始参数，在这个参数下变换参考图像，并计算与目标图像的差异;
• 调整参数，使上述差异达到最小值;
• 输出最优参数作为配准变换。

3. 数据来源

本实验使用在实验室拍摄的四张照片作为测试数据。

4. 实现算法

4.1 软件界面

本实验搭建了一个二维图像配准工具。"主界面"可进行配准的各项设置，并展示参考图像与目标图像:

"结果界面"会在运行配准之后弹出，分左中右三栏，分别为变换后的参考图像、变换后的参考图像与目标图像的对比(会交替显示两图像)、目标函数优化曲线:

4.2 优化算法

在上述基础上实现了三种优化算法，分别是暴力算法、梯度下降法、模拟退火算法。执行算法前，首先会在Registration::initialize()函数中针对不同变换方式，初始化变换的各项参数，例如旋转变换(四个参数):

params.resize(4); // cx cy angle scale
limits.resize(4);
steps.resize(4);
limits[0] = std::make_pair(-c / 2, c / 2);
limits[1] = std::make_pair(-r / 2, r / 2);
limits[2] = std::make_pair(-180, 180);
limits[3] = std::make_pair(0.5, 2);
for (int i = 0; i < params.size(); i++) {
  params[i] = limits[i].first;
}
steps[0] = 1;
steps[1] = 1;
steps[2] = 1;
steps[3] = 0.01;

其中params表示参数数组，limits表示各个参数的取值范围，steps表示各个参数邻域的步长。

暴力算法使用深度优先搜索遍历整个可行解空间，可求解出全局的最优解:

void Registration::optimizeNaive(int pos) {
    if (pos >= params.size()) {
        completeIteration();
        return;
    }
    for (float i = limits[pos].first; i <= limits[pos].second; i += steps[pos]) {
        params[pos] = i;
        optimizeNaive(pos + 1);
    }
}

梯度下降法在每轮迭代中，寻找当前最优解的邻域中最优的可行解，并以此代替最优解，直到邻域中没有比当前解更好的解为止。我实现的函数中，某个可行解的邻域定义为将第$i$个参数增加或减少steps[i]得到的可行解的集合($i$为某个随机参数的序号)。因梯度下降法与初始值选取密切相关，选择不当会导致陷入局部最优解，我在函数中增加了一个外循环，每轮循环随机选择一组初始值，并开始一轮新的梯度下降:

void Registration::optimizeGD() {
    for (int i = 0; i < 10000; i++) {
        double tmp_loss = 1e30;
        double cur_loss = 1e30;
        for (int j = 0; j < params.size(); j++) {
            // give random start values
            params[j] = random(limits[j].first, limits[j].second);
        }
        while (true) {
            std::vector<float> old = params;
            std::vector<float> best = params;
            for (int pos = 0; pos < params.size(); pos++) {
                params[pos] = old[pos] + steps[pos];
                optimizeGDHelper(best, cur_loss);
                params[pos] = old[pos] - steps[pos];
                optimizeGDHelper(best, cur_loss);
                params[pos] = old[pos];
            }
            if (tmp_loss == cur_loss) {
                break;
            }
            tmp_loss = cur_loss;
            params = best;
            completeIteration();
        }
        completeIteration();
    }
}
void Registration::optimizeGDHelper(std::vector<float>& best, double& cur_loss) {
    applyTransform();
    double s = getSimilarity(trans_img, tar_img, similarity_type);
    if (s < cur_loss) {
        cur_loss = s;
        best = params;
    }
}

模拟退火算法与梯度下降法相比，具有跳出局部最优解的能力。当温度较低时，算法不容易跳出局部最优解，从而退化为梯度下降法。所以，如果在温度下降到很低时没能跳入最优解附近，算法同样会收敛于局部最优解。因此，与上述梯度下降法的代码类似，我也使用了一层外循环，用于生成多个随机初始参数:

void Registration::optimizeSA() {
    for (int i = 0; i < 10000; i++) {
        for (int j = 0; j < params.size(); j++) {
            // give random start values
            params[j] = random(limits[j].first, limits[j].second);
        }
        applyTransform();
        double fi = getSimilarity(trans_img, tar_img, similarity_type);
        double temp = 10;
        double r = 0.99;
        int k = 1000;
        while (k > 0) {
            // select a random neighbor
            int len = params.size();
            int r1 = rand() % len;
            double r2 = random(1, 10) * ((rand() % 2) * 2 - 1);
            params[r1] += r2 * steps[r1];
            applyTransform();
            double fj = getSimilarity(trans_img, tar_img, similarity_type);
            double delta = fj - fi;
            if (random(0, 1) < exp(-delta / temp)) {
                // jump to this neighbor
                completeIteration(iter % 10 == 0);
                fi = fj;
            } else {
                params[r1] -= r2 * steps[r1];
            }
            temp *= r;
            k--;
        }
        completeIteration();
    }
}

4.3 实现细节

1. 计算图像变换、图像相似度需要遍历图像每一像素，这一过程是较为耗时的。对此，代码中将图像先缩小至原有尺寸的$1/10$，即，像素数变为原来的$1/100$，再进行变换、相似度计算。最终再恢复至原有图像尺寸。
2. 图像变换后，其周围可能出现空白区域(下图)，在计算图像相似度时，这一区域也是需要计算的(否则会倾向于生成全是边界的图像)，所以需要填补这一区域。我使用的方法是计算全图像的平均颜色，并以平均颜色填补。
   

5. 实验结果与分析

5.1 目标函数值可视化

下图展示了4号图片到1号图片的平移变换的优化目标函数值。平移变换有两个参数，所以可以映射到二维空间。使用上述暴力算法计算出参考图像平移至所有位置，计算其目标函数值，便得到如下的可视化结果。

从图中可以看出目标函数最优值(最小值)位于图像中间附近，有局部最优值，但该函数较为光滑，因此可认为性质较好。

5.2 结果

5.2.1 暴力算法

暴力算法效率较低，所以只进行了平移变换、旋转变换的实验。其中，平移变换约耗时3秒，旋转变换耗时数分钟。

该算法可以确保找到最优解，可作为其他算法的验证手段。

下图是枚举平移变换参数时，相似度函数的变化曲线:

5.2.2 梯度下降法

下图是使用梯度下降法优化仿射变换配准问题的结果，在迭代至5000轮左右时，已经找到较好的结果。这里，每一轮迭代指的是找到一个新的最优值，或开启一次新的梯度下降。

右侧的蓝色曲线表示梯度下降过程的目标函数值变化，可以看出，在每轮梯度下降过程内，目标函数值能够快速下降。

5.2.3 模拟退火算法

下图使用了模拟退火算法来计算5.2.2中相同的优化问题以便对比。该方法约在4000轮左右取得一个较好的解，更好的解出现在50000轮左右。

相比于梯度下降法，模拟退火在每轮迭代中不需要计算所有邻域解，而梯度下降计算数量为参数个数$\times 2$，所以模拟退火的迭代速度高于梯度下降。

为了使两次模拟退火之间可以区分，我在两次模拟退火的能量函数之间插入了一个峰值。观察蓝色曲线可以发现，每次运行模拟退火时，前期能量函数波动较大(因为跳到更差的解的几率很大)，在后期逐渐趋于稳定，类似梯度下降的结果。

5.3 参数调整

三种算法中，只有模拟退火算法有需要调整的参数。包括迭代次数k、初始温度temp与温度衰减系数r。

k控制何时停止迭代，根据曲线取值即可。

起初不知道如何设置初始温度与衰减系数，将温度设置得很高，衰减也很快，然后发现优化效果不如梯度下降法。分析各个参数的作用与模拟退火算法的优势之后，我得到结论:

• 初始温度应与目标函数值在同一量级;
• 衰减系数应调整到使得温度在进行了一半的迭代轮数之后下降到足够低，进而使得优化只向更优的可行解跳跃。

最终的曲线和期望一致，一开始有较大的波动，随后逐渐稳定直至收敛。

5.4 分析与结论

本次实验我得到了如下结论:

1. 经过对比试验，对于二维照片配准问题，使用2范数作为相似性度量函数的效果较好，其次是1范数。
2. 参数空间的邻域结构非常重要。例如仿射变换矩阵尺寸为$2\times3$，有六个参数，若直接优化这六个参数，则优化结果较差，很难收敛到最优解。原因可能是这六个参数的现实意义不够明确(例如，增大第一个参数会使图像怎样变化)，且参数之间存在较大关联性(例如，同时增大第一个、第五个参数，作用是放大图像，而只改变一个则图像会变得不自然)。仿射变换的另一种表示方式是，原图像的左上、右上、左下三个顶点分别变换到了新图像的哪些位置，同样需要六个参数描述，但这些参数的可解释性更强，目标函数也更容易优化至最优。实验发现，选择合适的参数能够让目标函数的性质更好，更容易收敛到全局最优解。
3. 模拟退火算法在图像配准问题上的优势不明显。对于特别不光滑的目标函数，梯度下降算法是无能无力的，因为选择一个随机初始值下降至全局最优值的概率非常低。而模拟退火算法有助于跳出这些波动，相比梯度下降法更容易得到最优解。但是对于图像配准问题，其目标函数性质较好，局部最优解较少。于是，多次使用梯度下降算法也能够找到全局最优解，模拟退火算法的优势难以体现。

实验室的公司委托项目，是一个面向临床医生的三维医疗影像系统软件，已逐步在多家医院推广使用。

Official Sitelink