CN101826209B - 一种基于Canny模型的三维医学图像分割方法 - Google Patents

Abstract

一种基于Canny模型的三维医学图像分割方法，属于图像处理技术领域。首先通过用户交互在原始三维医学图像I_mnk中截取包含用户感兴趣目标的三维区域图像I_mnk；然后对三维区域图像I_mnk进行中值滤波以去除图像噪声；接着再采用Canny方法获得中值滤波后三维区域图像I_mnk的Canny边缘信息图像C_mnk；之后根据用户在某一帧图像目标边缘附近的交互点坐标(i₀，j₀)、搜索目标边缘、提取完整的封闭的目标边缘；最后提取出三维区域图像I_mnk中所有帧的感兴趣目标边缘E_mnk作为三维医学图像的分割结果。采用本发明提出的基于Canny模型的三维医学图像分割方法，使用了较少的用户交互过程，计算量较小，能够快速有效地提取三维医学图像中感兴趣目标的边缘信息，从而完成三维医学图像的分割。

Description

一种基于Canny模型的三维医学图像分割方法

技术领域

本发明属于图像处理技术领域，涉及三维医学图像的交互式快速分割方法。

背景技术

为了准确的分辨医学图像中的正常组织结构和异常病变，需要对医学图像进行分割。传统的图像分割方法主要包括：

(1)基于边缘的分割方法：通常利用区域间不同性质(如区域内灰度不连续性)划分出各个区域之间的分界线，这类方法包括并行微分算子法(如Roberts、Sobel、Laplacian、Marr等算子)、串行边界搜索方法、基于曲面拟合的方法等；

(2)基于区域的分割方法：通常利用同一区域内的均一性识别图像中的不同区域，包括阈值法、区域生长和分裂合并、分类器和聚类、基于随机场的方法等。

但由于医学图像的成像原理和组织本身的特性差异，而且图像的形成受到诸如噪音、场偏移效应、局部体效应和组织运动，组织和组织之间，组织和器官之间，器官和器官之间等等之间的影响，快速准确的实现三维医学图像分割是非常必要的。因此，有必要针对医学应用这个领域，研究一种方法能够按照人们的主观意思自愿的去分割目标区域，而且又能快速准确的提取目标边缘信息。这样，有利于快速对感兴趣目标进行整体观察和准确分析。

发明内容

本发明提供一种基于Canny模型的三维医学图像分割方法，该方法不仅具有较快较准确的图像分割速度，而且具有人机交互功能，能够根据操作者的意愿分割得到感兴趣的目标边缘；同时，本发明便于观察到病变组织(或器官)形状等，有利于做出准确的医学分析及进一步判断。

为了方便地描述本发明内容，首先对一些术语进行定义。

定义1：感兴趣目标。待分割的对象物体，如肿瘤，肝脏等。

定义2：图像尺寸大小。二维医学切片图像的尺寸大小为M×N，M表示二维医学切片图像长度方向的像素点个数，N表示二维医学切片图像宽度方向的像素点个数。三维医学图像的尺寸大小为M×N×K，表示有K张尺寸大小为M×N的二维医学切片图像。

定义3：图像灰度。指图像中的亮度信息，没有任何颜色信息。

定义4：3×3中值滤波器。将当前像素点及其周围8个像素点的灰度值按从大到小排列，取排列在最中间的那个数值作为当前像素点的灰度值。

定义5：图像梯度。指图像灰度值变化的向量场，包括梯度大小和梯度方向。梯度大小指灰度值的最大变化率，梯度方向指灰度值增长最快的方向。

定义6：生长算法。首先选择一个初始像素点，然后向四周搜索满足一定规则、与初始像素点属于同一类别且相互连通的所有像素点。将这些像素点保存起来，即可得到所期望的结果。

定义7：边缘端点。指边缘线的断裂点。

定义8：边缘断裂方向。指边缘线在端点处的切线方向并指向断开的边缘位置。该方向与边缘端点处的梯度方向相垂直。

定义9：毛刺。指从边缘线上分叉出来的小短线，由端点和分叉点的连线构成。一般是由于噪声和细化过程造成的。

定义10：数学形态学。用具有一定形态的结构元素去量度和提取图像中对应形状以达到对图像分析和识别的目的。膨胀是其中一个基本运算，公式为：

其中A为图像集合，B为结构元素，^表示做关于原点的映射，( )_x表示平移x，⌒表示交集，φ表示空集，

为膨胀运算符。

本发明详细技术方案如下：

一种基于Canny模型的三维医学图像分割方法，如图1所示，包括以下步骤：

步骤1：截取包含感兴趣目标的三维区域图像。在尺寸为M×N×K的原始三维医学图像I_MNK中，截取一个包含感兴趣目标的、尺寸为m×n×k的三维区域图像I_mnk；其中m≤M，n≤N，k≤K。

本步骤将待处理的原始三维医学图像的尺寸大小就从原来的M×N×K缩小到了m×n×k，可加快图像处理速度，也有利于提高目标分割的准确率。

步骤2：对步骤1所得的三维区域图像I_mnk进行中值滤波。中值滤波的目的是为了减少噪声对三维区域图像I_mnk的影响；中值滤波时，可采用3×3大小或5×5大小的中值滤波器。

步骤3：计算中值滤波后的三维区域图像I_mnk所有像素点的梯度，得到三维区域图像I_mnk的梯度图像D_mnk。

将中值滤波后的三维区域图像I_mnk中任意一帧切片图像的像素点表示为(i，j)，像素点(i，j)的灰度值表示为I(i，j)，其中1≤i≤m，1≤j≤n；首先计算像素点(i，j)的横向梯度d_x(i，j)和纵向梯度d_y(i，j)，其中d_x(i，j)＝I(i+1，j)-I(i，j)，d_y(i，j)＝I(i，j+1)-I(i，j)；然后计算像素点(i，j)的梯度值大小M(i，j)和梯度方向θ(i，j)，得到三维区域图像I_mnk的梯度图像D_mnk，其中：

$M(i,j)=\sqrt{d_x^2(i,j)+d_y^2(i,j)},$ θ(i，j)＝arctan[d_x(i，j)/d_y(i，j)]

步骤4：提取梯度图像D_mnk的Canny边缘信息图像C_mnk。对梯度图像D_mnk中的每一帧梯度图像D_mn，进行如下操作：

步骤4-1：规定像素点(i，j)的梯度线，即将梯度方向属于(0，π/8]、(7π/8，9π/8]或(15π/8，2π]的像素点(i，j)归入水平线梯度像素点，将梯度方向属于(π/8，3π/8]或(9π/8，11π/8]的像素点(i，j)归入第一对角线梯度像素点，将梯度方向属于(3π/8，5π/8]或(11π/8，13π/8]的像素点(i，j)归入垂直线梯度像素点，将梯度方向属于(5π/8，7π/8]或(13π/8，15π/8]的像素点(i，j)归入第二对角线梯度像素点。

步骤4-2：采用非极大值抑制方法提取梯度值屋脊。取像素点(i，j)及位于像素点(i，j)的梯度线上相邻的两个像素点，比较三个像素点的梯度值大小，若像素点(i，j)的梯度值小于其梯度线上相邻两个像素点梯度值，则将像素点(i，j)的梯度值置为“0”，否则保留像素点(i，j)的梯度值，得到非极大值抑制后的图像N(i，j)，其宽屋脊带已细化成只有一个像素点宽，并保留了屋脊的高度值。

步骤4-3：采用双门限阈值方法实现Canny边缘信息提取。选取两个固定的阈值τ₁和τ₂，且τ₂＝2*τ₁，0＜τ₂＜M_max(i，j)，其中M_max(i，j)为该帧梯度图像D_mn中的最大梯度值；将非极大值抑制后的图像N(i，j)中所有像素点的梯度值与阈值τ₁或τ₂分别进行比较，将梯度值大于阈值τ₁或τ₂的像素点的梯度值置为“1”，否则置为“0”，分别得到两个阈值边缘图像T₁和T₂；然后以高阈值边缘图像T₂为基础，加入低阈值边缘图像T₁中与高阈值边缘图像T₂中边缘相连通的所有边缘点，得到该帧梯度图像D_mn的Canny边缘信息图像C_mn。

步骤4-4：对梯度图像D_mnk中的每一帧梯度图像D_mn进行步骤4-1至步骤4-3的操作，得到梯度图像D_mnk的Canny边缘信息图像C_mnk。

步骤5：用户交互。用户在三维区域图像I_mnk中间第t张图像I_mnt中感兴趣目标的边缘附近确定一个像素点位置(i₀，j₀)，其中0≤t≤k。

步骤6：确定三维区域图像I_mnk中间第t张图像I_mnt的感兴趣目标边缘E_mnt，具体包括以下步骤：

步骤6-1：搜索邻近边缘点位置。在相应的第t张Canny边缘信息图像C_mnt中，根据步骤5确定的像素点位置(i₀，j₀)，由近及远的向四周搜索边缘点。

步骤6-2：提取边缘线。以步骤6-1搜索到的边缘点作为起始点，通过生长算法，获取与其相连通的所有边缘点，从而提取整条边缘线。

步骤6-3：提取封闭的边缘线。若步骤6-2获得的是封闭的边缘线，则转入步骤6-4；否则，找到一条边缘线，然后沿着当前边缘线两个端点的断裂的方向分别查找另外的边缘线，将找到的边缘线的端点与当前边缘线相邻的端点用直线连接起来作为感兴趣目标边缘线；重复本步骤，直至完成封闭边缘线的提取。

步骤6-4：去除毛刺。在步骤6-3得到的封闭边缘线中，删除所有的毛刺边缘，得到三维区域图像I_mnk中间第t帧图像I_mnt的感兴趣目标边缘E_mnt。

步骤7：确定相邻帧Canny边缘信息图像中感兴趣目标边缘区域。将第t帧图像I_mnt的目标边缘E_mnt膨胀成为15～25像素宽度的边缘区域，作为第t+1帧图像I_mn(t+1)和第t-1帧图像I_mn(t-1)的感兴趣目标边缘搜索区域。

步骤8：搜索感兴趣目标边缘线。在第t+1张和t-1张Canny边缘信息图像C_mn(t+1)和C_mn(t-1)中对应的步骤7所确定的感兴趣目标边缘区域内分别搜索最长的边缘线。

步骤9：提取封闭目标边缘并去除毛刺。对于步骤8搜索到的第t+1帧Canny边缘信息图像C_nm(t+1)和第t-1帧Canny边缘图像C_mn(t-1)中的最长边缘线，依照步骤6-3和步骤6-4的处理过程，即可分别获得三维区域图像I_mnk中第t+1帧图像I_mn(t+1)的感兴趣目标边缘E_mn(t+1)和t-1帧图像I_mn(t-1)的感兴趣目标边缘E_mn(t-1)。

步骤10：重复步骤7至步骤9，直到获得整个三维区域图像I_mnk的感兴趣目标边缘E_mnk，就能从整个原始三维医学图像I_MNK中提取到感兴趣目标边缘E_mnk，从而完成三维医学图像的分割。

需要说明的是：

1、步骤1获得的包含用户感兴趣目标的三维区域图像I_mnk尺寸大小是由用户交互决定的。

2、步骤2中的中值滤波器窗口大小选择由图像噪声决定，窗口越大，去除噪声效果越好，但相应的计算量也会增加。

3、步骤3到步骤4实现了对三维区域图像I_mnk的Canny边缘的提取。由于Canny边缘能够准确的描述目标对象的准确边缘，这样就为后续的处理提供了准确保障。其中步骤4-3中所述双门限阈值τ₁和τ₂满足0＜τ₁＜τ₂＜M_max(i，j)的关系，其中M_max(i，j)为该帧梯度图像D_mn中的最大梯度值。通常情况下，若双门限阈值τ₁和τ₂取值越大(接近于M_max(i，j))，所提取到的Canny边缘就越少，可能会漏掉真实的边缘信息；若双门限阈值τ₁和τ₂取值越小(接近于0)，则所提取到的Canny边缘就越多，可能会将一些噪声误作为边缘信息而提取。因此，双门限阈值τ₁和τ₂取值应当适中，最好是τ₂＝2τ₁＝2/3M_max(i，j)。

本发明采用基于Canny边缘模型的方法，首先通过用户交互截取包含用户感兴趣目标的三维区域图像I_mnk；然后对三维区域图像I_mnk进行中值滤波以去除图像噪声；接着再采用Canny方法获得中值滤波后三维区域图像I_mnk的Canny边缘信息图像C_mnk；之后根据用户在某一帧图像目标边缘附近的交互点坐标(i₀，j₀)、搜索目标边缘、提取完整的封闭的目标边缘；最后提取出三维区域图像I_mnk中所有帧的感兴趣目标边缘E_mnk作为三维医学图像的分割结果。采用本发明提出的基于Canny模型的三维医学图像分割方法，使用了较少的用户交互过程，计算量较小，能够快速有效地提取三维医学图像中感兴趣目标的边缘信息，从而完成三维医学图像的分割。

附图说明

图1是本发明流程示意图。

具体实施方式

本发明技术方案在实现时，首先使用Matlab语言编写程序；然后使用三维MRI或CT医学图像数据进行参数设定和程序优化处理；最后使用C++语言重写程序代码及交互界面框架，以提高程序性能。

Claims (3)

1.一种基于Canny模型的三维医学图像分割方法，包括以下步骤：

步骤1：截取包含感兴趣目标的三维区域图像；在尺寸为M×N×K的原始三维医学图像I_MNK中，截取一个包含感兴趣目标的、尺寸为m×n×k的三维区域图像I_mnk；其中m≤M，n≤N，k≤K；

步骤2：对步骤1所得的三维区域图像I_mnk进行中值滤波；

步骤3：计算中值滤波后的三维区域图像I_mnk所有像素点的梯度，得到三维区域图像I_mnk的梯度图像D_mnk；

将中值滤波后的三维区域图像I_mnk中任意一帧切片图像的像素点表示为(i，j)，像素点(i，j)的灰度值表示为I(i，j)，其中1≤i≤m，1≤j≤n；首先计算像素点(i，j)的横向梯度d_x(i，j)和纵向梯度d_y(i，j)，其中d_x(i，j)＝I(i+1，j)-I(i，j)，d_y(i，j)＝I(i，j+1)-I(i，j)；然后计算像素点(i，j)的梯度值大小M(i，j)和梯度方向θ(i，j)，得到三维区域图像I_mnk的梯度图像D_mnk，其中：

$M(i,j)=\sqrt{d_x^2(i,j)+d_y^2(i,j)},$ θ(i，j)＝arctan[d_x(i，j)/d_y(i，j)]

步骤4：提取梯度图像D_mnk的Canny边缘信息图像C_mnk；对梯度图像D_mnk中的每一帧梯度图像D_mn，进行如下操作：

步骤4-1：规定像素点(i，j)的梯度线，即将梯度方向属于(0，π/8]、(7π/8，9π/8]或(15π/8，2π]的像素点(i，j)归入水平线梯度像素点，将梯度方向属于(π/8，3π/8]或(9π/8，11π/8]的像素点(i，j)归入第一对角线梯度像素点，将梯度方向属于(3π/8，5π/8]或(11π/8，13π/8]的像素点(i，j)归入垂直线梯度像素点，将梯度方向属于(5π/8，7π/8]或(13π/8，15π/8]的像素点(i，j)归入第二对角线梯度像素点；

步骤4-2：采用非极大值抑制方法提取梯度值屋脊；取像素点(i，j)及位于像素点(i，j)的梯度线上相邻的两个像素点，比较三个像素点的梯度值大小，若像素点(i，j)的梯度值小于其梯度线上相邻两个像素点梯度值，则将像素点(i，j)的梯度值置为“0”，否则保留像素点(i，j)的梯度值，得到非极大值抑制后的图像N(i，j)；

步骤4-3：采用双门限阈值方法实现Canny边缘信息提取；选取两个固定的阈值τ₁和τ₂，且0＜τ₁＜τ₂＜M_max(i，j)，其中M_max(i，j)为该帧梯度图像D_mn中的最大梯度值；将非极大值抑制后的图像N(i，j)中所有像素点的梯度值与阈值τ₁或τ₂分别进行比较，将梯度值大于阈值τ₁或τ₂的像素点的梯度值置为“1”，否则置为“0”，分别得到两个阈值边缘图像T₁和T₂；然后以高阈值边缘图像T₂为基础，加入低阈值边缘图像T₁中与高阈值边缘图像T₂中边缘相连通的所有边缘点，得到该帧梯度图像D_mn的Canny边缘信息图像C_mn；

步骤4-4：对梯度图像D_mnk中的每一帧梯度图像D_mn进行步骤4-1至步骤4-3的操作，得到梯度图像D_mnk的Canny边缘信息图像C_mnk；

步骤5：用户交互；用户在三维区域图像I_mnk中间第t张图像I_mnt中感兴趣目标的边缘附近确定一个像素点位置(i₀，j₀)，其中0≤t≤k；

步骤6：确定三维区域图像I_mnk中间第t张图像I_mnt的感兴趣目标边缘E_mnt，具体包括以下步骤：

步骤6-1：搜索邻近边缘点位置；在相应的第t张Canny边缘信息图像C_mnt中，根据步骤5确定的像素点位置(i₀，j₀)，由近及远的向四周搜索边缘点；

步骤6-2：提取边缘线；以步骤6-1搜索到的边缘点作为起始点，通过生长算法，获取与其相连通的所有边缘点，从而提取整条边缘线；

步骤6-3：提取封闭的边缘线；若步骤6-2获得的是封闭的边缘线，则转入步骤6-4；否则，找到一条边缘线，然后沿着当前边缘线两个端点的断裂的方向分别查找另外的边缘线，将找到的边缘线的端点与当前边缘线相邻的端点用直线连接起来作为感兴趣目标边缘线；重复本步骤，直至完成封闭边缘线的提取；

步骤6-4：去除毛刺；在步骤6-3得到的封闭边缘线中，删除所有的毛刺边缘，得到二维区域图像I_mnk中间第t帧图像I_mnt的感兴趣目标边缘E_mnt；

步骤7：确定相邻帧Canny边缘信息图像中感兴趣目标边缘区域；将第t帧图像I_mnt的目标边缘E_mnt膨胀成为15～25像素宽度的边缘区域，作为第t+1帧图像I_mn(t+1)和第t-1帧图像I_mn(t-1)的感兴趣目标边缘区域；

步骤8：搜索感兴趣目标边缘线；在第t+1张和t-1张Canny边缘信息图像C_mn(t+1)和C_mn(t-1)中对应的步骤7所确定的感兴趣目标边缘区域内分别搜索最长的边缘线；

步骤9：提取封闭目标边缘并去除毛刺；对于步骤8搜索到的第t+1帧Canny边缘信息图像C_mn(t+1)和第t-1帧Canny边缘图像C_mn(t-1)中的最长边缘线，依照步骤6-3和步骤6-4的处理过程，即可分别获得三维区域图像I_mnk中第t+1帧图像I_mn(t+1)的感兴趣目标边缘E_mn(t+1)和t-1帧图像I_mn(t-1)的感兴趣目标边缘E_mn(t-1)；

步骤10：重复步骤7至步骤9，直到获得整个三维区域图像I_mnk的感兴趣目标边缘E_mnk，就能从整个原始三维医学图像I_MNK中提取到感兴趣目标边缘E_mnk，从而完成三维医学图像的分割。

2.根据权利要求1所述的基于Canny模型的三维医学图像分割方法，其特征在于，步骤2对步骤1所得的三维区域图像I_mnk进行中值滤波时，采用3×3大小或5×5大小的中值滤波器。

3.根据权利要求1所述的基于Canny模型的三维医学图像分割方法，其特征在于，步骤4-3中所述阈值τ₂＝2τ₁＝2/3M_max(i，j)。

Images