基于改进的Live-Wire算法在ARPlanner中的应用
引言
患者进行放疗之前，需要对患者制定放疗计划，而精准的靶区及危及器官的勾画是精准放疗的前提[1]。在精准勾画前提下，靶区附近危及器官受照剂量才会更精确，从而提高治疗的精度，降低副反应发生率。医生在勾画靶区以及危及器官时需要反复修改，因而花费了大量时间，降低了勾画靶区的工作效率。一些危及器官虽然经过医生的反复修改，但是勾画的精度仍远远不够。
目前，靶区的勾画主要分为三种[2-4]：自动勾画、人机交互式半自动勾画、手动勾画。由于医学图像的复杂性，自动勾画难以得到准确的器官勾画结果，容易造成轮廓的偏离；手动勾画需要消耗医生大量的时间与精力，随着CT设备的不断提升，层数数量也在不断的增加，这给医生手动勾画带来了大量的工作，勾画一个患者往往需要数个小时甚至几天，需要不断的修改，由于耗时长，精准度又受主观的影响较大，极大地降低了医生的工作效率；人机交互式勾画，可以通过医生的经验配合图像处理的技术，因此本文采用Live-Wire算法对患者的靶区以及危及器官进行交互式勾画，可以极大地提高医生的工作效率与勾画的精准度。
1 方法
Live-Wire[5-7]首先是快速地定位到该图像区域的边缘上，通过与用户交互，更精确地完成整个感兴趣区域的勾画、分割。Live-Wire算法的实现通过两步，第一步是对构造成本函数的计算，第二步是对最短路径的搜索。首先将待分割的图像映射成一个加权有向图，图像中的像素点对应有向图中的结点。连接图像中的相邻像素的边对应于连接图中相邻结点的边，并且在边上定义构造函数，使得图像中的器官边缘具有较小的权重。给予非边缘区域更大的权重，然后根据图像的特征计算每个像素的Laplace过零点特征函数、梯度特征函数、梯度方向函数局部代价，然后，基于加权有向图寻路算法Dijkstra执行最短路径搜索，并根据初始种子点和目标点搜索最短路径，获得感兴趣区域的边缘。传统的Live-Wire算法的梯度幅值由水平方向和垂直方向来计算。本文对Live-Wire算法进行了改进，在梯度幅值[8-10]的计算由原算法中的水平方向和垂直方向改进为由水平方向、45°方向、垂直方向和135°方向来计算。
1.1 改进的构造成本函数
每个像素的构造成本函数[11]计算公式如(1)所示。
其中，公式(1)中函数l(p,q)表示像素p与相邻像素q之间的成本函数；函数fZ(q)表示像素q处的Laplace零交叉点特征函数；函数fG(q)表示像素q的梯度特征函数；函数fD(p,q)表示点(p,q)的梯度方向。经过多次实验，权值wZ=0.43,wG=0.43,wD=0.14。
1.1.1 Laplace过零点特征函数fZ(q)的计算
Laplace过零点特征函数[11-12]是一种二值化的边缘检测方法[13-16]。因为图像边缘的灰度变化较大，所以图像的一阶偏导数在边缘处有局部最大值或最小值，二阶偏导数在边缘处会通过“0”点。过“0”点指两个相邻的像素二阶偏导数由正变为负，或者由负变为正的点。这一正一负的两个像素点，值与“0”接近的点称为过“0”点，取值0。Laplace过零点特征函数计算得到的结果是一个像素宽的窄带，表示图像内目标的边缘。Laplace过零点特征函数中的过“0”点表示具有较强的边缘特征的点。假设IL(q)是图像I在q点处的函数值，则其公式如下：
其中，函数IL(q)是像素点q经过Laplace变换后的值，函数fZ(q)的作用是突出医学图像中器官的边缘。
1.1.2 改进的梯度特征函数fG(q)的计算
图像使用Laplace过零点特征函数计算得到的是一个二值矩阵，计算的结果只有极少的点值是0。由于计算得到的是一个二值矩阵，因此它仅仅能表示像素点是否属于目标边缘而不能说明像素点属于图像边缘概率的大小、梯度大小及边缘特征的强弱。因此在公式(1)中引入了一个表示像素点图像边缘特性强度的梯度特征函数[11,17]。有很多种方法的算子可以计算图像的梯度值，在本论文中采用的是高斯函数。由于在一副图像中每个像素都具有八个相邻的像素点，当边缘与检测方向垂直才能很好地检测到边缘点，因此在上述方法中只用到水平方向和垂直方向来计算梯度的幅值在一定条件下存在漏检情况，为了减少漏检，本文在计算梯度幅值方法上做了改进，将用水平方向Ix、45°方向 I45°、垂直方向 Iy、135°方向 I135°作为梯度幅值的计算方向，以下为四个方向的算子。
梯度G的计算公式(4)如下：
梯度值越大说明图像的边缘特征越强，那么像素点的梯度特征函数得到的权值就应该越小，
为使高梯度产生低能量，令
其中G´=G-min(G)，为了保证梯度特征最大值计算方式的统一，梯度特征函数值需要使用欧几里得距离进行处理。如果q是p的斜向的相邻点，fG(q)需要除以
；如果q是p的水平或垂直方向的相邻点，fG(q)需要除以1，如公式(6)所示。
1.1.3 梯度方向fD(p,q)的计算
在边缘提取中，保证图像边缘的平滑是必要的。在公式(1)中，梯度方向函数[17]就是用来控制目标边缘平滑的平滑约束因子。梯度方向函数给边缘方向变化较大的点赋予一个较高的权值，给边缘变化平滑的点赋予较小的权值。梯度方向值是一个由Ix和Iy决定的单位向量。假设D(p)是同p点处梯度方向垂直的单位向量，梯度方向函数的公式是：
其中
从公式(7)可以看出，如果与相邻点的梯度方向具有相同或者相似梯度方向的点具有较小的权值，而与相邻点梯度方向垂直的点具有较大的权值。其中公式(9)中的Ix(p)，Iy(p)，Ix(q)，Iy(q)分别表示像素p、像素q沿方向x、方向y的方向梯度。
由上述三个计算成本函数的公式可以方便的找到图像的边缘点。
1.2 计算最小成本路径
当成本函数确定之后，下一步就是使用加权有向图的路径寻找算法来寻找最佳路径。算法的思想：只需每次设置种子点，然后重新计算从种子点到其他结点的最短路径。也就是说，所有路径的集合构成从种子点扩散的最短路径树。在本文中，Dijkstra算法[18-21]用于找到最短路径，并且可以找到从单个源点到其他顶点的最短路径。
Dijkstra算法的输入为：有向图G和G中的源顶点S，用V来表示G中所有顶点的集合。用E来表示G中连接结点边的集合。(u,v)表示从顶点u到v的路径是连通的，边缘的权重由权重函数w:E→[0,∞]定义，因此，w(u,v)是顶点u到顶点v的非负成本（cost），两个顶点之间的距离可用于计算边的成本。寻路算法还可以找到从一个顶点到图中任何其他顶点的最短路径。如图1所示，结点A、B、C、D、E分别为平面上的5个结点，现在需要分别求出结点A到其他4个结点B、C、D、E的最短路径，优先队列Q={A、B、C、D、E}。
图1 有向图G及各边权重
寻路过程如下：
（1）设起点的距离为0，即键值为0，其他的点为∞，因为此时不知道如何到达，即键值为∞，从队列中取出最小值，取出队列中的最小值A，然后观察A，把A加到集合S中，A在队列里移出来，就不会再放回去了。
（2）接下来更新队列中其他结点的键值，观察跟A有边连接的点，从A到B有一条，它的权值是10，沿AB边走，它的花费是10，总权值是0+10=10，更新队列里B的权值为10，从A到C是另外一条边，0+3=3，更新C的键值3，然后在队列里找出最小值3即C点，说明这里的最短路径是：从A到C，从队列里删除C。
（3）然后看所有从C出来的边，有一条是到B的，权值为4+3=7，比之前的10小，所以从这条路径到B会更好，更新B的键值为7，这里还有一条从C到D的，费用为8的边加上3为11，更新D的键值为11，然后还有一条从C到E的，费用为3+2=5，更新E的键值为5，到现在所有顶点都有一条有限的路径相连了。从这个队列里找出最小值E，此时这是最短路径，到达E的方案是从A->C->E，称这是最短路径。
（4）然后看从E点出发的边，这里只有一条，从E到D，费用是5+9=14，比之前的11要大，所以不走这条路径。
（5）现在剩下B和D，A、C、E都被删除了，取出最小值B的键值为7，这里有一条路径B到C，费用是7+1=8，要比3大，接着B到D，费用为7+2=9，比11要小，此时我们找到了一条更短的路径，所以当前D的最短路径的权值为9，即这条路径是经过A->C->B->D，移除B点，现在只剩D点了。
（6）观察从D点出发的边，这里只有一条，从D到E，费用为9+7=16，比之前的5要大，所以E的最短路径已经确定了，现在删除D点，队列是空的。
其寻路结果，见表1。
表1 A点到B、C、D、E各点的最短路径
目的 最短路径 最短距离A=>A A->A 0 A=>B A->C->B 7 A=>C A->C 3 A=>D A->C->B->D 9 A=>E A->C->E 5
在本文中使用的算法是Dijkstra算法的一种变体，输入包括： 种子点，有向图。输出为： 输入图中的最小路径树，每一个结点指向它的前一个结点，沿着种子点到该结点的最小代价路径。还将为每个结点分配总成本，对应于从结点到种子的最小成本路径。在有向图中，每一个结点都有以下三种状态：初始状态、活动状态、扩展状态。当全部结点都被扩展时，本算法终止。当算法开始时，图中的全部结点进行初始化。当算法运行时，所有活动结点都保持在优先队列Q中，该优先队列根据结点到种子的当前总成本排序。本文的变体的Dijkstra算法流程，见图2。
图2 本文变体的Dijkstra算法流程图
2 结果
本文在验证算法过程中在200例不同癌种患者中随机抽取了患者的CT图像序列，实验效果表明改进后的勾画效果比原算法勾画效果更能满足临床要求，本文采取了4幅典型图像进行对比显示，见图3～6。
本文算法在原有算法的基础上对梯度幅值的计算进行了改进，从水平方向和垂直方向来计算梯度的幅值改进为从水平方向、45°方向、垂直方向、135°方向计算梯度幅值，添加医学图像在多个梯度方向都很重要，而原方法只能提取2个方向梯度，对对角线方向梯度不能有效提取，这将导致在边缘提取的时候无法实现精准提取，而本文中方法明显对对角线方向实现梯度提取，兼顾了多个方向，边缘提取的精度得到有效保障。改进之后的算法能够使边缘更精确的被检测出来，本文对改进的Live-Wire算法和原Live-Wire算法的效果进行了对比，通过对比图可以发现，改进的Live-Wire算法能够更精确的器官的边缘进行检测。本文的算法能够使靶区附近危及器官受照剂量更精确，提高治疗的精度，降低副反应发生率，提高了医生勾靶的工作效率，为临床勾画工作降低了时间成本。
图3 肺部勾画效果示例
注：a. 改进的Live-Wire勾画的肺部；b.原算法勾画的肺部。
图4 肝脏勾画效果示例
注：a. 改进的Live-Wire勾画的肝脏；b. 原算法勾画的肝脏。
图5 肾脏勾画效果示例
注：a. 改进Live-Wire勾画的肾脏；b. 原算法勾画的肾脏。
图6 脑部勾画效果示例
注：a. 改进Live-Wire勾画的脑；b. 原算法勾画的脑。
将本文改进的交互式自动分割算法Live-Wire应用到ARPlanner软件（ARPlanner软件为北京全域医疗技术集团有限公司开发的放射治疗计划轮廓勾画软件，与医用加速器配套，适用于肿瘤患者的放射治疗。支持患者CT/MRI图像输入，并作三维图像重建与显示；支持任意切面解剖图像的显示；支持PET/CT、CT/MR、增强/平扫图像之间的融合配准；支持手动勾画、半自动勾画、自动勾画靶区和器官轮廓；具有DICOM网关功能，支持设备间数据通讯）中，可以看到通过本文算法能够更精确的勾画出器官（表2）。通过与Monaco软件（Monaco是由Elekta开发生产的放射治疗计划软件，主要用于常规方法的轮廓绘制、影像处理、模拟、影像融合、计划优化等）半自动勾画功能进行对比，发现本文的算法勾画精确度要高于Monaco中勾画的精度。从而说明本文算法具有一定的临床应用价值。
表2 手动勾画与本文算法勾画以及Monaco中边缘检测勾画对比表
A R P l a n n e r手动勾画的器官M o n a c o中边缘探测勾画的器官A R P l a n n e r本文算法勾画的器官头部心脏 肝脏 盆腔 眼球 肾脏 肺
3 结论
本论文将改进的Live-Wire算法运用到实际的放射治疗计划轮廓勾画软件中，极大地提高了肿瘤医学勾画的精准度与速度，为医生勾画靶区以及危及器官提高了工作效率以及精准度，为后续做计划提供了精准的前提，从而提高了放疗的精准度。本文的主要贡献是：一是引入水平方向、45°方向、垂直方向、135°方向计算梯度幅值，添加医学图像在多个梯度方向都很重要，而原方法只能提取2个方向梯度，对对角线方向梯度不能有效提取，这将导致在边缘提取的时候无法实现精准提取，而本文中方法明显对对角线方向实现梯度提取，兼顾了多个方向，边缘提取的精度得到有效保障。二是根据Dijkstra原始算法进行变体，使得该算法在确定初始种子点以后，可以灵活自由的进行路径选择。与当前应用比较广泛的放射治疗软件Monaco中的半自动勾画功能进行了对比，发现本文算法的勾画效果要优于Monaco中的半自动勾画算法勾画的效果，有一定的临床应用价值。
[1]闫梦梦,王卫东,郎锦义.影像组学技术及其在肿瘤精准放疗中的应用[J].肿瘤预防与治疗,2018,31(5):66-70.
[2]刘晓东,李婷婷,陈海,等.放射治疗临床信息管理系统的需求分析与设计[J].中国医疗设备,2018,33(1):126-127.
[3]Antony J,Mcguinness K,Welch N,et al.An interactive segmentation tool for quantifying fat in lumbar muscles using axial lumbar-spine MRI[J].IRBM,2016,37(1):11-22.
[4]Ananth C,Sankari IU.Automatic image segmentation method based on ALG[J].Int J Innov Res Comput Commun Eng,2014,2(4):3716-3721.
[5]徐义春,刘军清,董方敏,等.基于多特征的人机交互式目标分割[J].广西大学学报（自然科学版）,2017,42(2):742-747.
[6]Li H,Dai Y,Zhang J.Improved Live-Wire algorithm for kidney image segmentation[A].2017 36th Chinese Control Conference[C].Dalian:IEEE,2017:10890-10895.
[7]Jessica D,Rosemary-Claire C.If biodiversity offsets are a dead end for conservation,what is the live wire? A response to Apostolopoulou & Adams[J].Oryx,2016,51(1):35-39.
[8]杨恢先,姜德财,谭正华,等.局部方向梯度幅值与相位差分的人脸识别算法[J].计算机工程与应用,2017,53(11):217-222.
[9]曹洪运,赵宇峰,高佳佳.基于梯度和标准差值的模糊边缘检测[J].自动化与仪表,2019,34(1):50-54.
[10]沈德海,鄂旭,张龙昌,等.一种基于梯度的细胞图像边缘检测算法[J].信息技术,2018,(3):6-9.
[11]徐晓刚,于金辉,马利庄.复杂物体轮廓提取[J].中国图象图形学报,2018,6(5):455-459.
[12]薛犇犇.半自动虚拟人切片图像分割方法研究[D].南京:南京理工大学,2011.
[13]王振华,胡伏原.基于显著性边缘的运动模糊图像恢复方法[A].第十五届中国体视学与图像分析学术会议论文集[C].南宁:中国体视学学会,2017.
[14]漆阳.CT图像边缘提取算法及质量评估分析研究[D].成都:西南交通大学,2016.
[15]范晞,费胜巍,储有兵.基于Canny算子的改进型图像边缘提取算法[J].自动化与仪表,2019,34(1):46-49.
[16]郭逸伦.基于OpenCV的边缘检测算法效率分析[J].科学技术创新,2019,(1):92-93.
[17]Mortensen EN,Barrett WA.Interactive segmentation with intelligent scissors[J].Graph Models Image Process,1998,60(5):349-384.
[18]Thomas HC,Charles EL,Ronald LR,Clifford S.算法导论（原书第3版）[M].殷建平,徐云,王刚,等,译.北京:机械工业出版社,2012:373-387.
[19]Krzysztof CC,Alexandre XF,Paulo AVM.Path-value functions for which Dijkstra’s algorithm returns optimal mapping[J].J Math Imaging and Vis,2018,60:1025-1036.
[20]Rini W,Cin C,Suhartono,et al.TransTrip: a shortest path finding application for Jakarta public transportation using Dijkstra algorithm[J].J Comput Sci,2018,14(7):939-944.
[21]岳晓娟.基于一般Dijkstra的改进算法在最短路径问题中的应用[J].农村经济与科技,2018,(4):289-290.
Application of Improved Live-Wire Algorithm in ARPlanner