CN104835169A - 一种前列腺图像融合方法 - Google Patents

Abstract

本发明公开了一种前列腺图像融合方法，所述方法包括图像采集与分割、图像标志点提取、计算映射函数f三个预处理步骤，最后根据f定位病灶区域在超声图像中的位置。本发明通过图像融合技术，能够有效提高活检系统精确度以及跟踪病灶区域的准确性。通过将核磁共振图像(MRI)和超声(US)图像的融合，将MRI图像中反映出的可疑病灶位置，在融合之后能够计算出在US图像中的位置，利用US的实时性来监控扎针位置；由此能够有效提高扎针的准确性，同时减少前列腺不必要的损伤。

Description

一种前列腺图像融合方法

技术领域

本发明涉及图像融合技术，具体涉及一种前列腺图像融合方法。

背景技术

传统的活检系统是由超声(US)引导，由于US所反映出的组织信息有限，只能够呈现出器官的大致轮廓，无法显示出器官内部的具体结构，不能够反映病灶区域；传统活检系统采用的是对前列腺区域均匀分布扎针方法，从而检查前列腺组织是否病变，例如：12针活检方法，该方法存在较大的误差以及缺陷，未能够将病灶区域与扎针位置关联起来，也存在定位精确性的争议。

随着成像技术以及图像融合技术的发展，核磁共振图像能够清晰反映组织内部结构以及病变区域，对于医疗诊断具有重要的参考意义；本发明将核磁共振(MRI)图像作为参考，利用图像融合算法将核磁共振图像与超声图像融合，算法的精度以及实效性影响着定位扎针的精确性，这是需要解决的困难。

现有技术中主要是通过超声(US)引导，利用均匀扎针方法来检测病灶区域；未能够精确定位病灶，存在精确度的争议性；并且均匀扎针的方法覆盖了前列腺的整个区域，导致有些未病变的前列腺组织在活检过程中受到不必要的损伤。

发明内容

针对上述问题，本发明旨在提供一种能够准确定位病灶区域且能有效减少对前列腺组织的不必要损伤的前列腺图像融合方法。

为实现该技术目的，本发明的方案是：一种前列腺图像融合方法，

所述方法包括图像采集与分割、图像标志点提取、计算映射函数f三个预处理步骤，最后根据f定位病灶区域在超声图像中的位置，具体步骤如下：

第一步，图像采集与分割，首先通过核磁共振设备以及超声设备，获取前列腺的MRI图像和US图像；然后对前列腺MRI和US图像进行分割；最后由分割后的图像形成前列腺3D图像，从而得到同一个病人前列腺在两种不图像的3D图像；

第二步，图像标志点提取，通过边缘检测提取所分割图像表面的像素点，MRI以及US图像分别有标志点；

第三步，图像标志点提取，由于US图像未能够提供前列腺内部组织结构，所以需要提取能够反映出前列腺整体轮廓的三维表面的点集作为标志点，标志点的提取或者舍去的要求：

a、取出单片已分割的前列腺图像数据，通过边缘检测提取出单片前列腺图像轮廓；

b、设置提出标志点的步长(step)，遍历前列腺轮廓图像，每隔step提取出该点作为标志点，直至单片标志点提取完成；

c、同样设置提取标志点片之间的间隔步长(step)，执行a、b直至所有片都处理完毕；

第四步，假设提取的标志点为X和V，需求解出两者之间的关联函数f(x)＝z，构造能量函数：

$E(M,f)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{a=1}{\overset{K}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}{\left|\right|x_i-f\left(v_a\right)\left|\right|}^2+\mathrm{\λ}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{a=1}{\overset{K}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}\log m_{\mathrm{ai}}-\mathrm{\ζ}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{a=1}{\overset{K}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}---\left(\text{1}\right)$

其中：(a＝1,2,......,K；i＝1,2,......,N；T为退火策略中的温度参数)

$y_a=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}{x}_i.$

$m_{\mathrm{ai}}=\frac{1}{T}\exp (-\frac{{(x_i-f\left(v_a\right))}^T(x_i-f\left(v_a\right))}{2T})---\left(2\right)$

$m_{K+1,i}=\frac{1}{T_0}\exp (-\frac{{(x_i-v_{K+1})}^T(x_i-v_{K+1})}{2T_0})---\left(3\right)$

$m_{a,N+1}=\frac{1}{T_0}\exp (-\frac{{(x_{N+1}-f\left(v_a\right))}^T(x_{N+1}-f\left(v_a\right))}{2T_0})---\left(4\right)$

初始化(1)中的参数，计算(2)(3)(4)，更新参数，直至能量函数收敛。

在能量函数中存在T控制参数，作为模拟退火参数，再函数收敛过程中不断减小，直至函数收敛。

第五步，根据f定位病灶区域在超声图像中的位置，由医生在MRI图像中勾画出的病变位置，通过形变函数f计算出在US图像中的位置；并将实际结果传送给定位扎针系统，从而准确定位病灶区域。

作为优选，所述第一步中，前列腺MRI和US图像的分割可以手动分割或者自动分割，并且由医生修改以及确认，并且由医生判定病变区域以及标记出。

本发明的有益效果，本发明通过图像融合技术，能够有效提高活检系统精确度以及跟踪病灶区域的准确性。通过将核磁共振图像(MRI)和超声(US)图像的融合，将MRI图像中反映出的可疑病灶位置，在融合之后能够计算出在US图像中的位置，利用US的实时性来监控扎针位置；由此能够有效提高扎针的准确性，同时减少前列腺不必要的损伤。

附图说明

图1为本发明的流程图。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细说明。

如图1所示，本发明所述的具体实施例为一种能够准确定位病灶区域且能有效减少对前列腺组织的不必要损伤的前列腺图像融合方法及装置。

为实现该技术目的，本发明的方案是：一种前列腺图像融合方法，

所述方法包括图像采集与分割、图像标志点提取、计算映射函数f三个预处理步骤，最后根据f定位病灶区域在超声图像中的位置，具体步骤如下：

第一步，图像采集与分割，首先通过核磁共振设备以及超声设备，获取前列腺的MRI图像和US图像；然后对前列腺MRI和US图像进行分割；最后由分割后的图像形成前列腺3D图像，从而得到同一个病人前列腺在两种不图像的3D图像；

第二步，图像标志点提取，通过边缘检测提取所分割图像表面的像素点，MRI以及US图像分别有标志点；

第三步，图像标志点提取，由于US图像未能够提供前列腺内部组织结构，所以需要提取能够反映出前列腺整体轮廓的三维表面的点集作为标志点，标志点的提取或者舍去的要求：

a、取出单片已分割的前列腺图像数据，通过边缘检测提取出单片前列腺图像轮廓；

b、设置提出标志点的步长(step)，遍历前列腺轮廓图像，每隔step提取出该点作为标志点，直至单片标志点提取完成；

c、同样设置提取标志点片之间的间隔步长(step)，执行a、b直至所有片都处理完毕；

第四步，假设提取的标志点为X和V，需求解出两者之间的关联函数f(x)＝z，构造能量函数：

$E(M,f)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{a=1}{\overset{K}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}{\left|\right|x_i-f\left(v_a\right)\left|\right|}^2+\mathrm{\λ}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{a=1}{\overset{K}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}\log m_{\mathrm{ai}}-\mathrm{\ζ}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{a=1}{\overset{K}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}---\left(\text{a}\right)$

其中：(a＝1,2,......,K；i＝1,2,......,N；T为退火策略中的温度参数)

$y_a=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{m}_{\mathrm{ai}}{x}_i.$

$m_{\mathrm{ai}}=\frac{1}{T}\exp (-\frac{{(x_i-f\left(v_a\right))}^T(x_i-f\left(v_a\right))}{2T})---\left(b\right)$

$m_{K+1,i}=\frac{1}{T_0}\exp (-\frac{{(x_i-v_{K+1})}^T(x_i-v_{K+1})}{2T_0})---\left(c\right)$

$m_{a,N+1}=\frac{1}{T_0}\exp (-\frac{{(x_{N+1}-f\left(v_a\right))}^T(x_{N+1}-f\left(v_a\right))}{2T_0})---\left(d\right)$

初始化(a)中的参数，计算(b)(c)(d)，更新参数，直至能量函数收敛。

在能量函数中存在T控制参数，作为模拟退火参数，再函数收敛过程中不断减小，直至函数收敛。

第五步，根据f定位病灶区域在超声图像中的位置，由医生在MRI图像中勾画出的病变位置，通过形变函数f计算出在US图像中的位置；并将实际结果传送给定位扎针系统，从而准确定位病灶区域。

为了便于修正误差，所述第一步中，前列腺MRI和US图像的分割可以手动分割或者自动分割，并且由医生修改以及确认，并且由医生判定病变区域以及标记出。

具体推算过程如下：

TPS_RPM算法计算过程是双向优化过程，既能够得到点对应关系<V，X>，同时能够计算出优化后的形变函数f。

假设两个点集V和X(在R²或者R³空间中)，V＝{v_a,a＝1,2,...,K}以及X＝{x_i,i＝1,2,...,N}；上述两点集分别属于前列腺不同形态的表面点集，可能存在一一对应的点集；用f表示由上述点集表征的不同形态之间的形变函数，由f可得到:u_a＝f(v_a)，其中由形变函数映射得到的点集U等于V或者是V的子集，可能存在一对多的情况。本算法为的就是能够求解出优化的f以及标志点之间的一一对应关系，踢出不存在对应关系的点集。

同时引进||Lf||²作为形变函数f的平滑度约束条件，在指定形变模型之后会有特定的形式表述。需要最小化下列能量函数：

$\underset{Z,f}{\min }E(M,f)=\underset{Z,f}{\min }\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ai}}{\left|\right|x_i-f\left(v_a\right)\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|\mathrm{Lf}\left|\right|}^2-\mathrm{\&zeta;}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ai}}---\left(1\right)$

其中： ${\mathrm{\&Sigma;}}_{i=1}^{N+1}{z}_{\mathrm{ai}}=1(a=\{\mathrm{1,2},...,K\left\}\right),{\mathrm{\&Sigma;}}_{a=1}^{K+1}{z}_{\mathrm{ai}}=1(i=\{\mathrm{1,2},...,N\left\}\right),z_{\mathrm{ai}}\&Element;\{0,1\}$ 即z_ai只能取0和1值；Z是一个二值对应矩阵，如果v_a与x_i存在对应关系则z_ai＝1，否则z_ai＝0。

如果按(1)进行优化存在两个优化的问题：点集对应关系线性分配不连续问题及形变函数最小二乘连续问题。所以为了解决(1)中存在的上述问题，应用了两个方法：软分配以及退火策略。

软分配：将二值对应矩阵Z转换为具有连续值的矩阵M，其中M元素值在[0,1]之间，对应矩阵M此时允许V与X两个点集之间存在模糊以及局部的匹配关系。如(2)式

退火策略：通过软分配建立了具有模糊性的对应关系，同时需要利用退火策略来控制这种模糊性；在(1)中添加用于控制对应关系的模糊性，T是作为一个温度参数，在优化的过程中通过不断减小T来控制，此时能量函数的最小化过程类似于物理学中的退火过程；每个温度的最优化结果都会作为下个温度阶段的初始值，整个过程是探索性的过程。

利上述两个方法可以将(1)变换为如(2)式：

$E(M,f)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}{\left|\right|x_i-f\left(v_a\right)\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}\log m_{\mathrm{ai}}-\mathrm{\&zeta;}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}---\left(\text{2}\right)$

其中：

M＝{m_ai}，M是一个NxK的矩阵，用于表示两个点集之间的对应关系，其值在[0,1]之间，同时满足 ${\mathrm{\&Sigma;}}_{i=1}^{N+1}{m}_{\mathrm{ai}}=1(a=\{\mathrm{1,2},...,K\left\}\right),$ ${\mathrm{\&Sigma;}}_{a=1}^{K+1}{m}_{\mathrm{ai}}=1(i=\{\mathrm{1,2},...,N\left\}\right),m_{\mathrm{ai}}\&Element;\left[\mathrm{0,1}\right];$ 第一项是表示当v_a与x_i存在对应关系，m_ai的值更为大。

第二项是用于限制形变函数，L只是标示符，此项可以用能够保证f扭曲能量的函数代替。

第三项是利用了退火策略而添加的一项，其中T是逐渐降低，每一温度的最小值作为下个温度阶段的初始值，直至满足收敛条件，T可以根据所需要收敛的程度以及计算的精确度来设定。

第四项是为了防止过多的点不能够建立对应关系。

其中λ以及ζ是控制各项之间的权重值。

需将(2)式收敛于某个最小值，过程如下：

1、更新点集之间的对应关系，即更新M＝{m_ai}(a＝{1,2,...,K},i＝{1,2...,N})

$m_{\mathrm{ai}}=\frac{1}{T}\exp (-\frac{{(x_i-f\left(v_a\right))}^T(x_i-f\left(v_a\right))}{2T})---\left(3\right)$

不存在对应关系的点判定为局外点，当a＝K+1或者i＝N+1：

$m_{k+1,i}=\frac{1}{T_0}\exp (-\frac{{(x_i-v_{k+1})}^T(x_i-v_{k+1})}{2T_0})---\left(4\right)$

$m_{a,N+1}=\frac{1}{T_0}\exp (-\frac{{(x_{N+1}-f\left(v_a\right))}^T(x_{N+1}-f\left(v_a\right))}{2T_0})---\left(5\right)$

并且对矩阵的每一行或者列利用归一化算法，做归一化计算直至满足收敛条件。

2、更新形变函数，如果只考虑(2)式中的形变函数f，并且不考虑局外点的影响，那么需要计算下列函数的最小值：

$\underset{f}{\min }E\left(f\right)=\underset{f}{\min }\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{\left|\right|y_a-f\left(v_a\right)\left|\right|}^2+\mathrm{\&lambda;T}{\left|\right|\mathrm{Lf}\left|\right|}^2---\left(6\right)$

其中： $y_a=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}{x}_i$

将形变过程设定为薄板样条函数：

f(v_a，d，w)＝v_a·d+φ(v_a)·w           (7)

将(2)式中的L项的用于限制TPS的扭曲能量函数表示：

$E_{\mathrm{TPS}}\left(f\right)=\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{\left|\right|y_a-f\left(v_a\right)\left|\right|}^2+\frac{\mathrm{\&lambda;}}{V}\&Integral;\&Integral;\&Integral;\left[\begin{array}{c}{\left(\frac{{\&PartialD;}^{2}f}{{\&PartialD;x}^2}\right)}^2+{\left(\frac{{\&PartialD;}^{2}f}{{\&PartialD;y}^2}\right)}^2+{\left(\frac{{\&PartialD;}^{2}f}{{\&PartialD;z}^2}\right)}^2+\\ 2{\left(\frac{{\&PartialD;}^{2}f}{\&PartialD;x\&PartialD;y}\right)}^2+2{\left(\frac{{\&PartialD;}^{2}f}{\&PartialD;x\&PartialD;z}\right)}^2+2{\left(\frac{{\&PartialD;}^{2}f}{\&PartialD;y\&PartialD;z}\right)}^2\end{array}\right]\mathrm{dxdydz}---\left(8\right)$

经过推导可以将(5)化为如下：

E_TPS(d，w)＝||Y-Vd-Φw||²+λ₁trace(w^TΦw)+λ₂trace[d-I]^T[d-I]   (9)

从而TPS-RPM算法如下过程：

step1:初始化T,λ₁,λ₂,M,d,w,以及整个函数收敛条件：相邻温度之间函数值的差。

step2:根据x_i对应点v_a之间通过f计算之后的误差，以及当前温度参数T(如(3-5)式)计算出M中的m_ai值，其中T₀表示初始温度参数。

step3:根据(6)式计算出形变函数参数(d,w)，利用上述参数计算出新的形变函数。

step4:逐渐减小T,λ₁,λ₂，如果满足收敛条件跳至step5，否则跳至step2重新计算。

step5：满足收敛条件，计算出优化后的f以及M。

具体使用时，先通过核磁共振设备获取前列腺的MRI图像，以MRI图像作为融合的参考；然后通过由于超声设备获取US图像，由于US图像中无法显示前列腺解剖结构特征，所以无法利用其内部特征来进行图像融合，但是能够呈现出其的整个形状结构特征，所以使用US图像作为轮廓参考；然后通过本系统的TPS-RPM算法(基于薄板样条的具有鲁棒性的点匹配算法)将MRI和US融合。该方法是基于前列腺形状表面特征的图像配准方法，该方法需要分别取出能够表征前列腺表面的足够多的点集。

通过提取出的表面点集，无需考虑器官内部组织结构，利用最优化方法求解出MRI图像形变至US图像的函数f，利用求解出的形变函数f，将医生在MRI图像中勾画出的可疑病灶区域位置，映射至US图像中，将病灶在US中的扎针点的坐标传递给自动定位装置，实现自动定位，由医生自主扎针。

通过该图像融合技术，能够有效提高活检系统精确度以及跟踪病灶区域的准确性。而且还能够通过将核磁共振图像(MRI)和超声(US)图像的融合，将MRI图像中反映出的可疑病灶位置，在融合之后能够计算出在US图像中的位置，利用US的实时性来监控扎针位置；由此能够有效提高扎针的准确性，同时减少前列腺不必要的损伤。

以上所述，仅为本发明的较佳实施例，并不用以限制本发明，凡是依据本发明的技术实质对以上实施例所作的任何细微修改、等同替换和改进，均应包含在本发明技术方案的保护范围之内。

Claims (2)

1.一种前列腺图像融合方法，其特征在于：

包括图像采集与分割、图像标志点提取、计算映射函数f三个预处理步骤，最后根据f定位病灶区域在超声图像中的位置，具体步骤如下：

第一步，图像采集，通过核磁共振设备以及超声设备，获取前列腺的MRI图像和US图像，以核磁共振提供MRI图像作为精标准，然后在扎针过程由超声设备提供US图像，所述US图像是由超声棒在移动中成像，其中切片直接的距离以及US图像应能够反映出前列腺的轮廓以及三维形态；

第二步，图像分割，首先根据采集前列腺样本作为前列腺分割的预判断，然后分割中通过学习算法预先选择一个最为相近的样本，最后调节图像使其能够涵盖前列腺组织的区域，由此构造出病人前列腺的实际三维图像模型；

第三步，图像标志点提取，由于US图像未能够提供前列腺内部组织结构，所以需要提取能够反映出前列腺整体轮廓的三维表面的点集作为标志点，标志点的提取或者舍去的要求：

a、取出单片已分割的前列腺图像数据，通过边缘检测提取出单片前列腺图像轮廓；

b、设置提出标志点的步长(step)，遍历前列腺轮廓图像，每隔step提取出该点作为标志点，直至单片标志点提取完成；

c、同样设置提取标志点片之间的间隔步长(step)，执行a、b直至所有片都处理完毕；

第四步，假设提取的标志点为X和V，需求解出两者之间的关联函数f(x)＝z，构造能量函数：

$E(M,f)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}{\left|\right|x_i-f\left(v_a\right)\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}\log m_{\mathrm{ai}}-\mathrm{\&zeta;}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{a=1}{\overset{K}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}---\left(1\right)$ 其中：(a＝1,2,......,K；i＝1,2,......,N；T为退火策略中的温度参数)

$y_a=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ai}}{x}_i.$

$m_{\mathrm{ai}}=\frac{1}{T}\exp (-\frac{{(x_i-f\left({\mathrm{\&upsi;}}_a\right))}^T(x_i-f\left({\mathrm{\&upsi;}}_a\right))}{2T})---\left(2\right)$

$m_{K+1,i}=\frac{1}{T_0}\exp (-\frac{{(x_i-{\mathrm{\&upsi;}}_{K+1})}^T(x_i-{\mathrm{\&upsi;}}_{K+1})}{2T_0})---\left(3\right)$

$m_{a,N+1}=\frac{1}{T_0}\exp (-\frac{{(x_{N+1}-f\left({\mathrm{\&upsi;}}_a\right))}^T(x_{N+1}-f\left({\mathrm{\&upsi;}}_a\right))}{2T_0})---\left(4\right)$

初始化(1)中的参数，计算(2)(3)(4)，更新参数，直至能量函数收敛。

在能量函数中存在T控制参数，作为模拟退火参数，再函数收敛过程中不断减小，直至函数收敛。

第五步，根据f定位病灶区域在超声图像中的位置，由医生在MRI图像中勾画出的病变位置，通过形变函数f计算出在US图像中的位置；并将实际结果传送给定位扎针系统，从而准确定位病灶区域。

2.根据权利要求1所述的前列腺图像融合方法，其特征在于：所述第三步中，前列腺MRI和US图像的分割可以手动分割或者自动分割，并且由医生修改以及确认，并且由医生判定病变区域以及标记出。