CN104268878B - 一种超声图像中运动探针检测和定位方法 - Google Patents

Abstract

本发明涉及一种超声图像中运动探针检测和定位方法，步骤如下：1)对超声序列图像向量化和矩阵化，得到超声影像数据矩阵，并交互式初始化；2)根据步骤1)形成的二维超声影像数据矩阵进行低秩建模和稀疏分解，得到矩阵低秩模型，并进行连续奇异像素支持区的交替检测，得到整体影像序列奇异像素支持区；3)利用步骤2)中得到的奇异像素支持区，转化其每一列向量为原始超声图像单帧大小，得到单帧稀疏块掩模区，并计算该掩模中的所有8联通区域，选取较大区域作为检测探针区域并进行空洞填充；4)通过计算检测到的探针区域中心点位置方法，找到探针具体空间位置信息。

Description

一种超声图像中运动探针检测和定位方法

技术领域

本发明涉及医疗及医学图像处理领域，更具体地说，涉及一种超声图像中运动探针检测和定位方法。

背景技术

超声影像运动探针监测和定位是医学图像导引手术治疗过程中的一项关键技术，其定位准确的有效性直接关系到手术的成败。

随着超声医学影像技术的发展及对疾病进一步进行研究和治疗的需要，如何利用图像分析引导技术从超声图像中快速精确的检测运动探针区域，并定量地获得介入式探针器械的相关位置信息，是提高手术成功率的关键，并已成为在超声影像介入式治疗中一个迫切需要解决的问题。

目前，大多数介入式治疗中运动探针位置的检测和定位是借助于超声医学影像技术，并由人工判断，需要依赖医生的临床经验，需耗费大量的人力资源。因此，让计算机自行对超声图像中的针状物进行检测和定位具有广泛的实际应用价值。

近年来，随着一些新兴图像分割技术的出现和发展，医学超声图像分割得到了迅速的发展。然而，这些分割方法主要集中应于器官、血管的分割，而对于类似探针器械物体的分割，国内外的研究比较少。虽然，国外有研究机构在CT和MRI图像中做过针状物的定位和分割的工作，但需加入过多的先验知识，自适应性不强，并且对噪声较为敏感，有一定局限性。此外，这些方法直接应用在超声图像中运动探针检测和定位的准确性和鲁棒性等方面还是有许多改进的地方，离实际临床应用还有很多工作要做。

一般而言，超声影像图像序列往往存在低对比度、高信噪比、边界模糊、探针与其他灰度相近肌能组织粘连等不利因素，使得准确的探针检测和定位仍是一个挑战性难题。此外，在超声影像图像形成的过程中，当超声波长与照射物体表面粗糙度相当时，就会极易产生斑点噪声。这些斑点噪声的存在使得超声图像质量较差，尤其是掩盖和降低了图像某些细节信息，在一定程度上干扰了探针位置的信息检测，为病情诊断及定量分析带来不利影响。

发明内容

本发明的目的在于克服现有技术的不足，提供一种能够克服超声图像中密度分布不均匀、低对比度、存在噪声和弱边界等不利影响的超声图像中运动探针检测和定位方法。

本发明的技术方案如下：

一种超声图像中运动探针检测和定位方法，步骤如下：

1)对超声序列图像向量化和矩阵化，得到超声影像数据矩阵，并交互式初始化；

2)根据步骤1)形成的二维超声影像数据矩阵进行低秩建模和稀疏分解，得到矩阵低秩模型，并进行连续奇异像素支持区的交替检测，得到整体影像序列奇异像素支持区；

3)利用步骤2)中得到的奇异像素支持区，转化其每一列向量为原始超声图像单帧大小，得到单帧稀疏块掩模区，并计算该掩模中的所有8联通区域，选取较大区域作为检测探针区域并进行空洞填充；

4)通过计算检测到的探针区域中心点位置方法，找到探针具体空间位置信息。

作为优选，步骤2)的矩阵低秩模型具体通过如下模型求解：

首先建立矩阵低秩数学模型，对于含有运动探针的超声影像图像序列，将每幅图像排列为一个向量；

然后将所有超声图像序列对应的向量排列成为一个矩阵，则稳定的背景部分对应于矩阵低秩部分，运动探针区域应于矩阵稀疏部分。

作为优选，设定运动探针为连续出现的异常信息进行求解，设ij表示影像序列矩阵中第j个列向量中的第i个像素，S∈{0,1}^p×F表示奇异像素支持区，其中p为列向量维度，F为序列帧数，则S_ij＝1表示ij位置为探针区域，S_ij＝0为背景区域；将探针检测的目的转化为在构建低秩模型的情况下，寻求最佳奇异像素支持区的估计过程。

作为优选，在奇异像素支持区的检测过程中，通过对稀疏矩阵每一列分量转化为原始超声影像图像矩阵大小，并进行平滑处理，达到消除细小奇异像素点的目的。

作为优选，进一步地，利用高斯模型建模，采用自适应阈值选择法选取目标稀疏块掩模，来进行连续奇异像素支持区的检测。

作为优选，一步地，步骤3)中，单帧检测到的单帧稀疏块掩模区，是通过求解奇异像素支持区得到的每一列向量的矩阵转化形式。

作为优选，步骤1)具体为：对于含有运动探针的超声影像图像序列，将每幅尺度为m×n的图像排列为一个向量，再将所有超声图像序列对应的向量排列成为一个矩阵D＝[vec(I₁)|…|vec(I_F)]∈R^p×F，其中p＝m×n为向量维度，F为序列帧数，则稳定的背景部分对应于矩阵低秩部分，而运动探针区域应于矩阵稀疏部分；

根据统计学习理论中低秩数学模型建立规则：将含破坏的数据矩阵D∈R^p×F分解为两矩阵之和，D＝A+E，其中，矩阵A是低秩的并逼近原始数据矩阵，E是稀疏的噪声数据矩阵；

于是矩阵恢复可用以下为凸的优化问题来描述：subjectto A+E＝D；其中，目标函数为矩阵A的秩以及噪声矩阵E的零范数，λ为噪声所占的权重。

作为优选，让表示在奇异像素支持区S的映射的情况下对矩阵X的映射变换，

设为其互补变换，则

由于矩阵的核范数是矩阵秩的包络，当矩阵的0范数和1范数等价时，矩阵恢复过程可以通过以下优化模型求解出低秩部分A和稀疏支持区S：

其中，核范数σ_k(A)表示矩阵A的第k个较大奇异值，参数λ为大于0的正数；

针对此优化函数为非凸的，并且含有约束条件，根据D＝A+E，可以利用以下能量优化函数进行低秩模型构建和奇异区域支持区检测：

作为优选，采取交替求解方法进行A和S的求解，步骤如下：

当给定初始估计的奇异像素支持区，则低秩矩阵A的恢复可以转化以下矩阵补全问题：

此问题描述为从部分观测数据中进行矩阵补全求解，采取SOFT-IMPUTE算法进行迭代求解：

其中，Θ_α＝UΣ_αV′，为奇异值阈值操作算子，Σ_α＝diag[(d₁-α)₊,…,(d_r-α)₊]，diag[·]为矩阵对角算子，UΣV′为的SVD分解，d₁,…,d_r为其奇异值；t₊＝max(t,0)，t₊是(d₁-α)₊中的元素取值算子。

作为优选，稀疏支持矩阵S由稀疏矩阵E的奇异值分布决定：

给定低秩矩阵则稀疏奇异矩阵E通过获得，将稀疏矩阵E进行如下绝对值归一化操作：

对于归一化后的稀疏矩阵中每一个t列分量转化为原始超声影像图像矩阵大小1≤x≤m，1≤y≤n，并利用9×9维零均值二维高斯滤波进行平滑处理，达到消除细小奇异点的目的：其中符号“*”为卷积操作算子；

针对获得的单帧前景稀疏块矩阵Em(x,y,t)，利用高斯模型建模：并采用自适应阈值选择法选取目标稀疏块掩模，其中阈值τ＝μ+σ，则稀疏块掩模可通过以下式子得到：

其中，1≤x≤m，1≤y≤n，1≤t≤F，p＝m×n；

通过矩阵与向量之间的转化得到：

采取交替求解方法进行和的求解，直到收敛为止。

作为优选，根据向量和矩阵转换规则，针对每一个与原始影像图像矩阵大小相同的稀疏掩模矩阵进行8连通域检测，并将8连通区域作为第当前第t帧目标探针检测区域并进行空洞填充。

作为优选，根据得到的探针检测区域结果，通过计算检测到的探针区域中心点位置方法，找到探针具体位置(x^*,y^*,t)，其中

本发明的有益效果如下：

对于超声影像图像序列进行数学建模，模型易于理解，且符合实际的物理意义。

利用矩阵秩最小化及稀疏分解方面研究的最新成果，采用高效的凸优化方法使得方法实现起来简单而高效。

本发明中的运动探针检测方法能够有效克服弱边界及噪声等不利影响，尤其是因为引入了连续奇异像素支持区的检测原理，能够很好的区分开运动探针和相邻器官及软组织，有效的避免了误检测显现出现，可用于超声介入式手术治疗。

本发明要解决的技术问题是，克服现有技术的不足，提供一种能够克服灰度不均匀、低对比度、噪声弱边界等不利影响的超声图像运动探针检测和定位方法。

本发明附加的方面和优点将在下面的描述中部分的给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

图1是本发明具体实施流程图；

图2是原始超声影像输入单帧图像；

图3是超声影像序列进行低秩建模后单帧背景部分；

图4是超声影像序列进行低秩建模后单帧奇异检测区域；

图5是稀疏奇异区域进行零均值二维高斯滤波平滑效果图；

图6是利用自适应阈值选择法选取目标稀疏块掩模效果图；

图7是目标8联通区域选择模效果图；

图8是检测到探针区域效果图；

图9是定位到探针位置效果图。

具体实施方式

以下结合图及实施例对本发明进行进一步的详细说明。

图1是本发明具体实施流程，具体实施过程如下：

对于含有运动探针的超声影像图像序列，将每幅尺度为m×n图像排列为一个向量，m为图像的长度，n为图像的宽度,再将所有超声图像序列对应的向量排列成为一个矩阵D＝[vec(I₁)|…|vec(I_F)]∈R^p×F，其中p＝m×n为向量维度，F为序列帧数，则稳定的背景部分对应于矩阵低秩部分，而运动探针区域应于矩阵稀疏部分。

根据统计学习理论中低秩数学模型建立规则：将含破坏的数据矩阵D∈R^p×F分解为两矩阵之和，即D＝A+E，其中矩阵A是低秩的并逼近原始数据矩阵，E是稀疏的噪声数据矩阵，于是矩阵恢复可用以下的凸的优化问题来描述：

其中目标函数为矩阵A的秩以及噪声矩阵E的零范数(即E的非零元素的个数)，λ表明噪声所占的权重。

针对运动探针检测，本发明视其为连续出现的异常奇异信息求解。令ij表示影像序列矩阵中第j个列向量中的i个像素，S∈{0,1}^p×F表示奇异像素支持区，例如S_ij＝1表示ij位置为探针区域，S_ij＝0为背景区域。因此，探针检测的目的，可以转化为在构建低秩模型A的情况下，寻求最佳S的估计过程。

让表示在奇异像素支持区S的映射的情况下对矩阵X的映射变换，

设为其互补变换，则

由于矩阵的核范数是矩阵秩的包络，矩阵的0范数和1范数(本发明定义为所有元素绝对值的和)在一定的条件下可以等价。因此，本发明提出通过以下优化模型求解出求解低秩部分A和稀疏支持区S：

其中，核范数σ_k(A)表示矩阵A的第k个较大奇异值，参数λ为大于0的正数。

针对此优化函数为非凸的，并且含有约束条件，根据D＝A+E，本发明中提出利用以下能量优化函数进行低秩模型构建和奇异区域支持区检测：

此在，本发明采取交替求解方法进行A和S的求解。给定初始估计的支持矩阵S，则低秩矩阵A可以转化以下矩阵补全问题：

此问题描述为从部分观测数据中进行矩阵补全求解，本发明采用了R.Mazumder等人“R.Mazumder，T.Hastie，and R.Tibshirani，“Spectral Regularization Algorithmsfor Learning Large Incomplete Matrices，”Journal of machine learning research，vol.11，pp.2287–2322，2010“提出的SOFT-IMPUTE算法进行迭代求解：

其中Θ_α为奇异值阈值操作算子：Θ_α＝UΣ_αV′，Σ_α＝diag[(d₁-α)₊,…,(d_r-α)₊]，diag[·]为矩阵对角算子，UΣV′为X经典的SVD分解，d₁,…,d_r为其奇异值；t₊＝max(t,0)，t₊是(d₁-α)₊中的元素取值算子。

由于稀疏支持矩阵S是用来标识奇异检测区域，可由稀疏矩阵E的奇异值分布决定。给定低秩矩阵则稀疏奇异矩阵E可以通过获得。一般情况下，异常奇异值对应的区域相比稳定的背景会有较大偏差，针对此特性，本发明将稀疏矩阵E进行如下绝对值归一化操作：

对于归一化后的稀疏矩阵中每一个列分量转化为原始超声影像图像矩阵大小，并利用9×9维零均值二维高斯滤波进行平滑处理，达到消除细小奇异点的目的：其中符号“*”为卷积操作算子，本实例中δ＝2。

针对获得的单帧前景稀疏块矩阵利用高斯模型建模：并采用自适应阈值选择法选取目标稀疏块掩模，其中阈值τ＝μ+σ。则稀疏块掩模可通过以下式子得到：

其中1≤x≤m，1≤y≤n，1≤t≤F，p＝m×n。

通过矩阵与向量之间的转化，可以得到：

采取交替求解方法进行和的求解，直到收敛为止。

根据向量和矩阵转换规则，针求解得到的奇异像素支持区得到其每一帧t与原始影像图像矩阵大小相同的稀疏掩模矩阵进行8连通域检测，并选取较大联通区区域作为检测探针区域

通过计算检测到的探针区域中心点位置方法，找到探针具体位置(x^*,y^*,t)，其中

本实例采用选择连续时间段内的超声影像30帧图像序列(F＝30)作为输入，如图2所示(实例为第影像序列的第5帧)，每一帧超声图像尺寸大小为240×240。

将每幅输入超声图像排列为一个向量，再将所有超声图像序列对应的向量按照时序关系堆叠成为一个数据矩阵D＝[vec(I₁)|…|vec(I₃₀)]∈R^57600×30。

将形成的二维矩阵进行低秩建模和稀疏分解。

先初始化稀疏支持区矩阵S为0，利用SOFT-IMPUTE方法进行低秩求解然后通过获得稀疏奇异矩阵E并进行绝对值归一化操作；紧接着，如图4所示，将归一化后稀疏矩阵E中每一个57600×1列分量转化为原始超声影像图像矩阵大小240×240。如图5所示，针对第t帧序列，利用9×9维零均值二维高斯滤波G_(0,δ)进行平滑操作1≤x≤240，1≤y≤240，达到消除细小奇异点的目的。

利用获得的单帧平滑稀疏块矩阵对其高斯模型建模：并采用自适应阈值选择法选取目标稀疏块掩模，其中阈值τ＝μ+σ。

其中1≤x≤240，1≤y≤240，1≤t≤30。

通过进一步将240×240维单帧目标稀疏掩块掩模矩阵转化为57600×1维向量，并综合所有帧序列可以得到输入序列的奇异像素支持区：得到的

采取交替求解方法进行和的求解，直到收敛为止，如图3和6所示。

根据向量和矩阵转换规则，针对每一个与原始影像图像矩阵大小相同的稀疏掩模矩阵进行8连通域检测。如图7所示，并选取较大8连通区域作为第5帧目标探针检测区域并进行空洞填充。

通过检测到的探针区域，如图8所示，通过计算中心点位置方法，找到第5帧探针具体位置(x^*,y^*,5)。如图9所示，坐标求得探针空间位置信息为(128，97)，与实际探针位置相符，实验结果证明了本文方法的有效性。

上述实施例仅是用来说明本发明，而并非用作对本发明的限定。只要是依据本发明的技术实质，对上述实施例进行变化、变型等都将落在本发明的权利要求的范围内。

Claims (11)

1.一种超声图像中运动探针检测和定位方法，其特征在于，步骤如下：

1)对超声图像序列向量化和矩阵化，得到二维的超声影像数据矩阵，并交互式初始化；

2)根据步骤1)形成的超声影像数据矩阵进行低秩建模和稀疏分解，得到矩阵低秩模型，并进行连续奇异像素支持区的交替检测，得到超声图像序列奇异像素支持区；

3)利用步骤2)中得到的奇异像素支持区，转化其每一列向量为超声图像单帧大小，得到单帧稀疏块掩模区，并计算该单帧稀疏块掩模区中的所有8连通区域，选取最大面积的连通区域作为检测探针区域并进行空洞填充；

4)通过计算空洞填充后的检测探针区域的中心点位置方法，得到探针具体空间位置信息。

2.根据权利要求1所述的超声图像中运动探针检测和定位方法，其特征在于，步骤1)中，对于含有运动探针的超声图像序列，将每幅图像排列为一个列向量；

步骤2)的矩阵低秩模型具体通过如下模型求解：

2.1)首先建立矩阵低秩数学模型；

2.2)然后对超声影像数据矩阵利用矩阵低秩数学模型进行建模，则稳定的背景部分对应于矩阵低秩部分，运动探针区域对应于矩阵稀疏部分。

3.根据权利要求2所述的超声图像中运动探针检测和定位方法，其特征在于，步骤2.2)具体为：设定运动探针为连续出现的异常信息进行求解，设ij表示超声影像数据矩阵中第j个列向量中的第i个像素，S∈{0,1}^p×F表示连续奇异像素支持区，其中p为列向量维度，F为序列帧数，则S_ij＝1表示ij位置为探针区域，S_ij＝0表示ij位置为背景区域；将探针检测的过程转化为在构建矩阵低秩模型的情况下，寻求最佳连续奇异像素支持区的估计过程。

4.根据权利要求2所述的超声图像中运动探针检测和定位方法，其特征在于，步骤2)中，在连接奇异像素支持区的交替检测过程中，进行平滑处理，达到消除细小奇异像素点的目的。

5.根据权利要求4所述的超声图像中运动探针检测和定位方法，其特征在于，进一步地，步骤2)中，对单帧稀疏块掩模区利用高斯模型建模，采用自适应阈值选择法选取单帧稀疏块掩模区，来进行连续奇异像素支持区的交替检测。

6.根据权利要求3所述的超声图像中运动探针检测和定位方法，其特征在于，步骤2.2)具体为：对于含有运动探针的超声图像序列，将每幅尺度为m×n的图像排列为一个列向量，再将超声图像序列对应的列向量排列成为一个超声影像数据矩阵D＝[vec(I₁)|…|vec(I_F)]∈R^p×F，其中p＝m×n为向量维度，F为序列帧数，则稳定的背景部分对应于矩阵低秩部分，而运动探针区域对应于矩阵稀疏部分；

根据统计学习理论中矩阵低秩数学模型建立规则：将含破坏的数据矩阵D∈R^p×F分解为两矩阵之和，D＝A+E，其中，矩阵A是低秩的并逼近原始数据矩阵，E是稀疏的噪声数据矩阵；

于是矩阵恢复可用以下为凸的优化问题来描述：subjectto A+E＝D；其中，目标函数为矩阵A的秩以及稀疏矩阵E的零范数，λ为噪声所占的权重。

7.根据权利要求6所述的超声图像中运动探针检测和定位方法，其特征在于，寻求最佳连续奇异像素支持区的估计过程具体为：让表示在奇异像素支持区S的映射的情况下对矩阵X的映射变换，

设为其互补变换，则

由于矩阵的核范数是矩阵秩的包络，当矩阵的0范数和1范数等价时，矩阵恢复过程通过以下优化模型求解出矩阵A和奇异像素支持区S：

其中，核范数σ_k(A)表示矩阵A的第k个奇异值，参数λ为大于0的正数；

针对此优化函数为非凸的，并且含有约束条件，根据D＝A+E，利用以下能量优化函数进行矩阵A和奇异像素支持区S的求解：

8.根据权利要求7所述的超声图像中运动探针检测和定位方法，其特征在于，采取交替求解方法进行A和S的求解，步骤如下：

当给定初始估计的奇异像素支持区，则矩阵A的恢复转化以下矩阵补全问题：

此问题描述为从部分观测数据中进行矩阵补全求解，采取SOFT-IMPUTE算法进行迭代求解：

其中，Θ_α＝UΣ_αV′，为奇异值阈值操作算子，Σ_α＝diag[(d₁-α)₊,…,(d_r-α)₊]，diag[·]为矩阵对角算子，UΣV′为的SVD分解，d₁,…,d_r为矩阵Z的奇异值；t₊＝max(t,0)，t₊是(d₁-α)₊中的元素取值算子。

9.根据权利要求8所述的超声图像中运动探针检测和定位方法，其特征在于，奇异像素支持区S由稀疏矩阵E的奇异值分布决定：

给定低秩矩阵则稀疏奇异矩阵E通过获得，将稀疏矩阵E进行如下绝对值归一化操作：

<mrow> <mover> <mi>E</mi> <mo>&amp;OverBar;</mo> </mover> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mi>E</mi> <mo>|</mo> <mo>-</mo> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <mo>|</mo> <mi>E</mi> <mo>|</mo> <mo>)</mo> </mrow> </mrow> <mrow> <mi>max</mi> <mrow> <mo>(</mo> <mo>|</mo> <mi>E</mi> <mo>|</mo> <mo>)</mo> </mrow> <mo>-</mo> <mi>min</mi> <mrow> <mo>(</mo> <mo>|</mo> <mi>E</mi> <mo>|</mo> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>;</mo> </mrow>

对于归一化后的稀疏矩阵中每一个t列分量转化为超声图像单帧大小1≤x≤m，1≤y≤n，并利用9×9维零均值二维高斯滤波进行平滑处理，达到消除细小奇异点的目的：其中符号“*”为卷积操作算子；

针对获得的单帧前景稀疏块矩阵利用高斯模型建模：并采用自适应阈值选择法选取单帧稀疏块掩模区，其中阈值τ＝μ+σ，则单帧稀疏块掩模区可通过以下式子得到：

<mrow> <mover> <mi>S</mi> <mo>&amp;OverBar;</mo> </mover> <mi>m</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>1</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mtable> <mtr> <mtd> <mrow> <mi>i</mi> <mi>f</mi> </mrow> </mtd> <mtd> <mrow> <mover> <mi>E</mi> <mo>&amp;OverBar;</mo> </mover> <mi>m</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;GreaterEqual;</mo> <mi>&amp;tau;</mi> </mrow> </mtd> </mtr> </mtable> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>w</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow>

其中，1≤x≤m，1≤y≤n，1≤t≤F，p＝m×n；

通过矩阵与向量之间的转化得到：

采取交替求解方法进行和的求解，直到收敛为止。

10.根据权利要求9所述的超声图像中运动探针检测和定位方法，其特征在于，步骤3)中，根据向量和矩阵转换规则，针对每一个与超声影像数据矩阵大小相同的单帧稀疏块掩模区进行8连通区域检测，并将8连通区域作为第当前第t帧检测探针区域并进行空洞填充。

11.根据权利要求10所述的超声图像中运动探针检测和定位方法，其特征在于，步骤4)中，根据得到的探针检测区域结果，通过计算检测到的探针区域中心点位置方法，得到探针具体空间位置信息(x^*,y^*,t)，其中