CN104299238A - 一种基于医学图像的器官组织轮廓线提取方法 - Google Patents

Abstract

本发明提供一种基于医学图像的器官组织轮廓线提取方法，首先通过对图像进行光滑处理，去掉图像中低振幅的点，加强显示显著的边，同时保留原图像的结构，然后在光滑处理基础上的图像进行Canny边缘提取，最后在提取出来的一些散乱点的基础上，应用二维散乱点自动生成算法形成符合美学观点的图像。本发明主要包括：L0梯度最小化进行光滑处理；边缘检测；二维散乱点自动生成算法。通过此算法可以生成封闭的，不相交的轮廓线，特别对尖锐角处理效果非常好。

Description

一种基于医学图像的器官组织轮廓线提取方法

技术领域

本发明涉及一种在医学图像中的器官组织轮廓中进行轮廓线提取的方法。

背景技术

最早提出的经典边缘检测方法大多从图像高频分量的增强出发，微分运算就成了边缘提取的主要方法。一阶微分算子采用梯度模的近似形式，包括RobertS算子、sobel算子、Prewitt算子、Laplacan算子以及Krish算子等，这些方法多是以将要处理像素为中心的邻域作为进行灰度分析的基础，实现对图像边缘的提取，并已经取得了较好的效果。但是这类方法也存在获得的边缘不连续、不封闭、像素宽、噪声干扰较严重等缺点。因此，这势必要求在对医学图像进行边缘提取时，设计出一种能够形成美观的边缘连续的封闭的噪声干扰比较小的方法。

由于微分运算在增强边缘的同时也会增强噪声，所以在微分运算之前需要用滤波降噪以减少噪声的影响，Marr和Hildreth提出LOG算法，利用高斯滤波器对图像进行平滑，使用各向同性的拉普拉斯算子进行求导，判定导数的零交叉点为边缘点。Johncanny提出Canny算子，并给出了边缘检测的三条准则。Canny算子也是一阶算子，实质是用一个准高斯函数做平滑运算，然后用带方向的一阶微分算子来定位导数最大值，其贡献在于给出了一种寻找最优算子的思路，即确立了最优准则。

在此背景下，本发明设计出一种形成美观的连续的封闭的噪声干扰比较小的边缘的方法，为医学图像的处理提供了帮助。

发明内容

本发明解决的技术问题是：降低了现有的医学图像中轮廓线提取的噪声干扰，同时绘制出美观的连续的轮廓线。

本发明的技术方案为：一种基于医学图像的器官组织轮廓线提取方法，包括以下三个步骤：

步骤(1)、L0梯度最小化进行光滑处理：L0梯度最小化算法通过限制邻域像素的强度变化的离散数量，在保留图像整体信息并且不影响整体精度的同时，通过增加陡度的过渡提高最突出的边缘，通过对输入图像进行处理进行L0梯度最小化光滑处理，增强整个图像的光滑度。

步骤(2)、边缘检测：图像的边缘检测，就是用离散化梯度逼近函数根据二维灰度矩阵梯度向量来寻找L0梯度最小化之后的图像的图像灰度矩阵的灰度跃变位置，然后将这些灰度跃变位置连接起来就形成了边缘，在图像的边缘图像中提取出了梯度变化比较大的点。

步骤(3)、二维散乱点自动生成算法：通过将边缘图像中的点集进行轮廓分组，然后对每一个轮廓的点集进行二维散乱点生成算法。通过点集的坐标生成没有叶子节点的欧几里得最小生成树，以此为基础生成不闭合的，轮廓长度最小化的满足格式塔原理的轮廓。最后，将每个轮廓的点集的处理结果显示在svg文件中得到整个图像的处理结果。

本发明的原理在于：

(1)在进行L0梯度最小化光滑处理时，我们用计算稀疏梯度数量的方法，通过限制强度变化点的周围的像素点的离散的数量来保持L0范数中的信息稀疏的追求。通过增强的边如果他们的结构清楚即使在低像素的图片或者细边缘也能很好的保持原图的边的特性。

(2)在进行边缘处理时，通过计算灰度矩阵的梯度向量，依靠离散化梯度逼近函数来找到图像中灰度矩阵的灰度跃变位置，然后将找到的这些位置用线连起来就构成了边缘。

(3)二维散乱点生成算法就是连接二维点生成美观的轮廓线，这些二维散乱点只给出了坐标信息。通过把对轮廓线的限制转化为一个遵循格式塔原理的问题，我们期望的轮廓Bmin是一个不相交的封闭的并且最小化边界长度的轮廓。算法的基础是没有叶子节点的欧几里得最小生成树BCmin。BCmin和Bmin能通过调整拓扑约束的参数近似表示。BCmin的近似的算法BC0能用一个贪婪算法快速的计算出来。BC0通过膨胀操作将度数大于2的点去除，在经过雕刻操作将内部度数等于0的点剔除掉就可以满足Bmin的拓扑和最小化需求，通过这两个步骤得到了最后结果Bout。通过大量的实验证明，Bout和Bmin是非常相似的。

本发明与现有技术相比的优点在于：

1、本发明使用的图像预处理技术，为接下来的轮廓线提取提供了很好的输入。

2、效率高：本发明中计算Bmin实际上是一个NP问题，然而计算Bout是线性时间此算法对那些有急剧变化的拐角处处理结果也非常好。

3、提取出的轮廓线对噪声干扰比较小。

附图说明

图1医学图像中的器官组织轮廓线提取方法的总体流程；

图2a：输入医学图像 b：L0梯度最小化光滑处理结果；

图3a：图像去噪结果 b：canny边缘检测结果；

图4二维点分组结果；

图5a:点集 b:BC0 c:BCinf d:Bmin；

图6点的分类示意图；

图7a:BC0 b:detail c:inflated；

图8a:Initial b:Exposed v0 c:Exposedv3；

图9a：点集 b：二维散乱点算法生成轮廓线；

图10输入图像的轮廓线提取结果。

具体实施方式

图1给出了医学图像中的器官组织轮廓中进行轮廓线提取的方法的总体处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种医学图像中的器官组织轮廓中进行轮廓线提取的方法，主要步骤介绍如下：

1、L0梯度最小化进行光滑处理

在进行L0梯度最小化光滑处理时，我们用计算稀疏梯度数量的方法，通过限制强度变化点的周围的像素点的离散的数量来保持L0范数中的信息稀疏的追求。通过增强的边如果他们的结构清楚即使在低像素的图片或者细边缘也能很好的保持原图的边的特性。在问题求解过程中，通过引入辅助变量交替最优化的策略求解，引入变量扩展原来的变量并且反复的更新。在解决这个问题的时候，由于是离散的数据，又引入了新的子问题，每个子问题都有解决方案，由于解决这个问题是很棘手的，所以我们的解决方法是一个很好的保持原图像的特性并且很好的展示突出的结构的近似方法。

在一维图像中，平滑在全局方式实现，通过限制非零梯度的数量提高最高对比度边缘。在一维信号中，g是输入的离散的信号，f是光滑处理的结果。公式(1)离散的计算了幅度的变化。公式(1)中p和p+1是两个相邻像素的索引，|fp-fp+1|是p的前向差分形式的梯度，#{}是计算操作符，用来输出满足|fp-fp+1|0的点的个数，也就是梯度的L0范数。c(f)不是计算梯度值，因此在改变对比度的时候结果不会受到影响。

c(f)＝#{p||fp-fp+1|≠0}   公式(1)

c(f)单独使用不会发挥作用，需要结合一定的约束条件，约束条件用来保证结果f和输入信号g在结构上相似，如公式(2)所示。公式(2)中c(f)＝k代表光滑结果f中存在着k个非零梯度值。图1.1e图是k＝6的结果，结果的信号图压缩了细节并且锐化了主要的边，并且和原图信息保持一致，k的取值越大，越限制与原图的近似程度，然而能刻画最好的显著的拐角。公式(2)中的(fp-gp)2是为了保证与原图保持最好的相似性，这样低幅度的结构可以用可控的统计的方式首先被移除，这个方法的好处是不管k设为何值，由于不采用局部滤波和平均值的操作，结果都不会出现模糊的现象。

$\underset{f}{\min }\underset{p}{\mathrm{\Σ}}(\mathrm{fp}-\mathrm{gp})2s.t.c\left(f\right)=k$    公式(2)

实际中，公式(2)中k的取值在不同像素的二维图像中可以从几十到几千，为了控制k的取值，采用了一种在结构压缩和保持相似度之间的平衡的方法。如公式(3)所示。λ是一个控制c(f)的权值，实际上是一个光滑参数，λ越大，结果图的边越少，通过实验证明，k与1/λ之间的是k关于1/λ单调递增。

$\underset{f}{\min }\underset{p}{\mathrm{\Σ}}(\mathrm{fp}-\mathrm{gp})2+\mathrm{\λ}\·c\left(f\right)$    公式(3)

在二维图像中，I代表输入图像，S代表结果图像，梯度是每个像素p与其邻像素在x和y方向的差异，如公式(4)所示，公式(4)计算了的个数，通过在一维上的分析，可以得出计算S的计算公式(5)。在彩色图像中是在rgb中的梯度的总和，∑(S-I)2限制图像结果的相似度。

$C\left(S\right)=\#\left\{p\right||\∂\mathrm{xSp}|+|\∂\mathrm{ySp}|\≠0\}$    公式(4)

$\underset{s}{\min }\left\{\underset{p}{\mathrm{\Σ}}(\mathrm{Sp}-\mathrm{Ip})2+\mathrm{\λ}\·C\left(S\right)\right\}$    公式(5)

公式(5)的求解，我们用通过引入辅助变量交替最优化的策略求解，通过引入变量扩展原来的变量并且反复的更新。在解决这个问题的时候，由于是离散的数据，又引入了新的子问题，每个子问题都有解决方案，由于解决这个问题是很棘手的，所以我们的解决方法是一个很好的保持原图像的特性并且很好的展示突出的结构的近似方法，我们引入了两个辅助变量hp和vp，hp相当于vp相当于公式(5)可以写成公式(6)，公式(6)中C(h，v)＝#{p||hp|+|vp|≠0}，β是一个控制变量(h，v)和他们梯度的相似度的自动调节参数。当足够大时，公式(6)近似等于公式(5)。公式(6)通过交替最小化(h，v)和S。在每一步中的变量都固定为上一步迭代后的值。

$\underset{S,h,v}{\min }\{\underset{p}{\mathrm{\Σ}}(\mathrm{Sp}-\mathrm{Ip})2+\mathrm{\λ}\·C(h,v)+\mathrm{\β}((\∂\mathrm{xSp}-\mathrm{hp})2+(\∂\mathrm{ySp}-\mathrm{vp})2)\}$    公式(6)

子问题1：计算S

计算S的子问题相当于最小化公式(7)，公式(7)通过省略公式(6)中不涉及S的项，此公式是二次方程式，因此有最小值，在快速傅里叶变换之后对角化微分算子来快速解决这个问题，如公式(8)，公式(8)中F是傅里叶变换操作符，F()*是共轭复数操作符，F(1)是狄克拉函数的傅里叶变换，相比于公式(7)中要涉及大量的转置运算，公式(8)的快速傅里叶变换的运算要快很多。

$\underset{p}{\mathrm{\Σ}}(\mathrm{Sp}-\mathrm{Ip})2+\mathrm{\β}((\∂\mathrm{xSp}+\mathrm{hp})2+(\∂\mathrm{ySp}-\mathrm{vp})2)\}$    公式(7)

$S=F^{-1}\left(\frac{F\left(I\right)+\mathrm{\β}(F(\∂x)*F\left(h\right)+F(\∂y)*F\left(v\right))}{F\left(1\right)+\mathrm{\β}(F(\∂x)*F(\∂x)+F(\∂y)*F(\∂y))}\right)$    公式(8)

子问题2：计算(h，v)

计算(h，v)的子问题相当于公式(9)，其中C(h，v)返回|h|+|v|不等于0的元素的个数，这个复杂的问题实际上可以用hp和vp单独估计的空间分解的方法快速解决。这样公式(9)可以写成公式(10)。公式(10)中H(|hp|+|vp|)是一个二分函数，如果|hp|+|vp|≠0就返回1，否则返回0.公式(10)中的每一个像素p的公式可表示为公式(11)，通过数学证明可以证明在下面条件下Ep取得最小值Ep*，如公式(12).通过计算每个像素的最小值Ep*的和∑pEp*计算公式(10)。

$\underset{h,v}{\min }\{\underset{p}{\mathrm{\Σ}}(\∂\mathrm{xSp}-\mathrm{hp})2+(\∂\mathrm{ySp}-\mathrm{vp})2)+\frac{\mathrm{\λ}}{\mathrm{\β}}C(h,v)\}$    公式(9)

$\underset{p}{\mathrm{\Σ}}\underset{\mathrm{hp},\mathrm{vp}}{\min }\{(\mathrm{hp}-\∂\mathrm{xSp})2+(\mathrm{vp}-\∂\mathrm{ySp})2+\frac{\mathrm{\λ}}{\mathrm{\β}}H\left(\right|\mathrm{hp}|+|\mathrm{vp}\left|\right)\}$    公式(10)

$\mathrm{Ep}=\{(\mathrm{hp}-\∂\mathrm{xSp})2+(\mathrm{vp}-\∂\mathrm{ySp})2+\frac{\mathrm{\λ}}{\mathrm{\β}}H\left(\right|\mathrm{hp}|+|\mathrm{vp}\left|\right)\}$    公式(11)

   公式(12)

为了提高算法的效率，我们选用Matlab来实现这一过程，由于我们整体***是用C++代码编写，因此需要在C++中调用Matlab计算的结果。图2展示了处理结果。

2、边缘检测

边缘检测的方法已经有很多，本文中采用的边缘检测的方法是canny边缘检测。通过边缘检测可以将尖锐的梯度变化大的角点提取出来。然而，在应用边缘检测之前，我们需要运用去噪处理，将图像中一些信噪比低的点去掉，从而使边缘检测可以得到很好的处理结果。

高斯滤波是将输入图像的每一个像素点与高斯内核做卷积操作，并且将卷积的和当作输出的像素值。通过高斯函数的图像的中间大两边小的性质，我们就可以有效的将噪点去除掉。我们使用了openCV中的cvSmooth函数进行去噪处理。通过将函数的光滑处理的类型设置为CV_GAUSSIAN，就可以进行高斯滤波处理了。

通过计算灰度矩阵的梯度向量，依靠离散化梯度逼近函数来找到图像中灰度矩阵的灰度跃变位置，然后将找到的这些位置用线连起来就构成了边缘。在实际情况中，理想的灰度改变不多见，其中大量的噪声干扰，滤波成为了我们工作的第一步。边缘检测分为三个步骤：滤波处理，增强像素点邻域强度变化大的点，剔除掉一些梯度变化大但不是我们想要的点。由去噪处理后的结果进行处理，我们使用了openCV中的cvCanny函数进行边缘检测的处理。cvCanny函数只能处理单通道图像，所以在处理之前需要把图像变为单通道图像，我们采用的方法是在加载图像的时候cvLoadImage函数时候将最后一个参数设置为0，就可以直接加载单通道图像。图3展示了处理结果。

3、二维散乱点自动生成算法

a)二维点分组的计算

由于二维散乱点生成轮廓线的算法生成的是密闭的轮廓，首先要将上述处理后的点集进行分组处理，将图像中的点集进行分组处理，首先将图像通过控制阀值将图像变换为二值图像，在二值图像中找到各个轮廓，找到值为1的像素点然后在其邻域中进行查找，如果是值为1的像素点，则将点加入这个轮廓中，否则进行下一个轮廓的查找。反复进行这个查找过程，直到所有的点都被遍历。这样既可以把所有的点集分成若干个轮廓。

由处理后的二值图像进行二维点的分组处理，我们使用openCV中的cvFindContours函数进行处理，输入是上述步骤中的二值图像，由此函数处理之后我们可以得到各个分组轮廓的指针，通过对各个轮廓的每个点进行访问，可以得到二维散乱点算法的输入。通过此函数的处理，将canny处理后的图像进行了轮廓的分组处理。图4展示了处理结果。

b)构建点集P的Delaunay三角剖分图DG(P)

三角剖分是指，假设V是二维实数域上额有限点集，边e是由点集中的点作为端点构成的封闭线段，E为e的集合。那么该点集V的一个三角剖分T(V，E)是一个平面图G，该平面图满足：

1、除了端点，平面图中的边不包含点集中的任何点。

2、没有相交边。

3、平面图中所有的面都是三角面，且所有三角面的合集是散点集V的凸包。

本文中用到的三角剖分是Delaunay三角剖分，它是一种特殊的三角剖分。当前仅当T中的每个三角形的外接圆的内部不包含V中任何的点，T是V的一个Delaunay三角剖分。

对于给定的点集P构建Delaunay三角剖分图，此处使用的是CGAL库中的Delaunay_triangulation_2类来进行三角剖分的构建。

c)构建理想轮廓线BCmin的一个近似轮廓BC0

由点集的Delaunay三角剖分图，构建BC0，采用的算法是一个贪婪算法，而且时间复杂度是线性复杂度，BC0通过保证每个点的度数大于2很好的遵循了格式塔原则中的封闭的原则，但是却不能保证好的连续性，由于和封闭的要求矛盾，但是在下面会有一些列的处理步骤，保证其连续性。

首先，将DG(P)中的边按照边的长度从小到大放在队列PQ中，当BC0不是一个连接的图形或者BC0中的点的度数有小于2的时候，将PQ中的第一个边从队列中移除，如果这个边连接了BC0中的两个连接的图形，或者这条边和BC0中的叶顶点连接着，那么久把这条边***到BC0中，直到BC0成为一个连接的图形，并且BC0中的点的度数都大于等于2。

d)膨胀操作

在膨胀操作的时候需要构建BC0的封闭轮廓Be，在BC0中所有满足下面的任一条件的边组成了Be。T(Be)是Be中所有的在DG(P)的三角形。

1、这条边是点集P的凸包中的边；

2、这条边是DG(P)中的三角形t的一条边，t有至少一条边在凸包上；

3、从DG中有一条边在凸包上的三角形开始，能通过遍历相邻三角形的边到达，并且不通过BC0中的任何边。

我们定义Be中点Vi的一个伞是点Vi的相邻的两条边按照一个给定的方向所夹的结构。对于点集P中的Vi能通过伞的数目u进行分类，如果u等于0，那么Vi是内点，如果u等于1，那么Vi是Be中的流点，否则，Vi不是合格的点。如图6中，虚线的边在DG中，实线在BC0中，这里也是Be，黑点是流点，灰点(v0是内点，v1，v2，v3是不合格点)，v3有两个伞。

我们定义膨胀候选三角形是所有Be内部的连接不合格的点的三角形ti。下面通过图7展示膨胀操作，将图7(a)红色部分放大，在图7(b)中放大的部分只有一个不合格的点v0，DG如虚线所示，Be如实线所示，可知t0和t1都是膨胀操作的候选三角形，由于增加t0中不再Be中的边，并且减去t0中已经在Be中的边即Be的长度变化||t0||最小，t0被添加到T(Be)中，结果为图7(c)这时候v0现在是流点，所以t1不再是候选三角形。通过膨胀操作可以将Be中所有的大于2的点去除掉了。

e)雕刻操作

在膨胀操作之后的结果BCinfl的封闭边界Be是流形轮廓，但是在Be中还会有一些内点。实际上有一条边在Be上的三角形中，在Be上的边的对面一定是内点，这样可以保证Be仍然是流形轮廓并且它的所有内点都可以暴露在Be上。

我们定义Be的雕刻候选三角形是在Be内部的三角形，并且三角形的一条边在Be上，它对面的点在Be的内部。下面通过图8展示雕刻操作，粗线连接的是Be，细线是DG中的线，灰色的点是内点，雕刻操作的候选三角形用黄色表示，在图8(a)中，有10个雕刻候选三角形，其中||t0||是最小的，在图8(b)中，通过把t0从T(Be)中移除，v0***到Be中，在v0点，t1，t2就成为了候选三角形，之后选择移除t3，使v1暴露出来，在图8(c)中，t4，t5之后被移除，从而把v2，v3暴露出来，从而形成Bout。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (4)

1.一种基于医学图像的器官组织轮廓线提取方法，其特征在于包括以下三个步骤：

步骤(1)、输入医学图像，用L0梯度最小化进行光滑处理：L0梯度最小化算法通过限制邻域像素的强度变化的离散数量，在保留图像整体信息并且不影响整体精度的同时，通过增加陡度的过渡提高最突出的边缘，通过对输入图像进行处理进行L0梯度最小化光滑处理，增强整个图像的光滑度；

步骤(2)、边缘检测：图像的边缘检测，就是用离散化梯度逼近函数根据二维灰度矩阵梯度向量来寻找L0梯度最小化之后的图像的图像灰度矩阵的灰度跃变位置，然后将这些灰度跃变位置连接起来就形成了边缘，在图像的边缘图像中提取出了梯度变化比较大的点；

步骤(3)、二维散乱点自动生成算法：通过将边缘图像中的点集进行轮廓分组，然后对每一个轮廓的点集进行二维散乱点生成算法；通过点集的坐标生成没有叶子节点的欧几里得最小生成树，以此为基础生成不闭合的，轮廓长度最小化的满足格式塔原理的轮廓；最后，将每个轮廓的点集的处理结果显示在svg文件中得到整个图像的处理结果。

2.根据权利要求1所述的医学图像中的器官组织轮廓线提取方法，其特征在于：步骤(1)中所述的L0梯度最小化进行光滑处理方法，使用计算稀疏梯度数量的方法，通过限制强度变化点的周围的像素点的离散的数量来保持L0范数中的信息稀疏的追求，结合图像稀疏性的特点，利用l0-范数作为正则项约束，保证了恢复图像的稀疏性要求；结构清楚的边通过增强操作，即使是低像素的图片的边或者细边缘也能很好的保持原图的边的特性。

3.根据权利要求1所述的医学图像中的器官组织轮廓线提取方法，其特征在于：步骤(2)中所述的边缘检测方法，该方法采用的是canny算法，在进行canny算法处理之前，首先进行滤波处理，去掉噪声，使图像处理结果光滑；采用的滤波方法是高斯滤波，在使用滤波之后使用canny边缘检测，有效的将一些噪点去除。

4.根据权利要求1所述的医学图像中的器官组织轮廓线提取方法，其特征在于：步骤(3)中所述的二维散乱点自动生成算法，由于二维散乱点生成轮廓线的算法最后生成的是全局的封闭的不相交的轮廓线，所以在此算法之前需要将canny算法的处理结果的图像的点集进行分组，然后对每个分组进行二维散乱点生成算法生成轮廓线，最后再将每个分组的处理结果显示在一个SVG文件中；所以在此算法之前需要将canny算法的处理结果的图像的点集进行分组，首先将图像通过控制阀值将图像变换为二值图像，在二值图像中找到各个轮廓，找到值为1的像素点然后在其邻域中进行查找，如果是值为1的像素点，则将点加入这个轮廓中，否则进行下一个轮廓的查找；反复进行这个查找过程，直到所以的点都被遍历；这样就可以把所有的点集分成若干个轮廓然后对每个分组进行二维散乱点生成算法生成轮廓线，最后再将每个分组的处理结果显示在一个SVG文件中。