CN108198235A - 一种三维超声重建方法、装置、设备及存储介质 - Google Patents

Abstract

本发明适用计算机技术领域，提供了一种三维超声重建方法、装置、设备及存储介质，该方法包括：在多角度超声扫描到的堆图像中将运动误差最小的堆图像设置为模板堆图像，将剩余堆图像与模板堆图像进行全局配准，根据剩余堆图像、模板堆图像构建三维体数据，通过基于核回归函数的体重建方式对该三维体数据进行更新，判断更新后的三维体数据是否收敛，是则将剩余堆图像、模板堆图像与三维体数据进行局部配准，否则跳转至对该三维体数据进行更新的步骤，判断局部配准后的三维体数据是否收敛，是则输出局部配准后的三维体数据，否则跳转至构建三维体数据的步骤，从而实现了基于多角度扫查的三维超声重建，有效地提高了三维超声图像的重建质量。

Description

一种三维超声重建方法、装置、设备及存储介质

技术领域

本发明属于医疗图像技术领域，尤其涉及一种三维超声重建方法、装置、设备及存储介质。

背景技术

三维超声***可以很直观的向医生提供感兴趣器官或组织在体内的空间位置以及三维形态，在临床上更准确地反映出人体器官病变的形状、体积、轮廓、以及器官间的邻接。目前，三维超声***主要包括两大类，一类是基于二维面阵探头(也称为三维超声专用容积探头)的容积成像方法，该方法将换能晶片和驱动装置包容于探头内，机械装置驱动晶片做等距的扇扫或环形扫描，另一类是利用传统的二维超声设备结合特定的空间定位信息来获取一序列的二维超声图像，然后利用三维超声重建方法生成三维超声体数据，主要包括基于机械臂的三维超声***(Motorized 3D ultrasound)、无***的三维超声***(Sensorless 3D ultrasound)和自由式三维超声***(Freehand 3D ultrasound)。从三维超声的临床应用角度来看，Freehand三维超声***更符合医生操作习惯和手术室华宁，是目标应用前景比较广的一种方案。

目前，三维超声重建算法可以分为三类：基于体素的方法(Voxel Based Method，VBM)、基于像素的方法(Pixel Based Method，PBM)和基于函数的方法(Function BasedMethod，FBM)。这三类重建算法都未考虑到两类关键性问题，一是医生在对单一角度扫查到的一个二维超声序列图像进行重建后的结果进行观察时，可能会由于骨组织和其它器官的遮挡，无法得到很好的诊断定论，即单一角度扫查的局限性，二是运动伪影对超声数据的影响，运动伪影的来源主要包含因不可遮挡的生理因素所带来的运动伪影(例如呼吸运动，虽然可在全身麻醉的条件下进行超声扫描来避免呼吸运动所造成的伪影，但是使用麻醉剂对患者有一定的伤害)、扫描受试者在行动上不配合扫查实验(例如：婴儿的不自主运动和胎动)、以及超声探头的挤压或重力场等自然因素造成的软组织等形变。

发明内容

本发明的目的在于提供一种三维超声重建方法、装置、图像处理设备及存储介质，旨在解决由于现有技术的三维超声重建未考虑到单一角度扫查的局限性和运动伪影对超声数据的影响，导致三维超声图像的重建质量不佳的问题。

一方面，本发明提供了一种三维超声重建方法，所述方法包括下述步骤：

当接收多角度超声扫描到的堆图像时，将所有堆图像中运动误差最小的堆图像设置为模板堆图像，并将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准；

根据所述剩余堆图像和所述模板堆图像，通过预设的基于点扩展函数的体重建方式构建三维体数据；

通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新，判断更新后的所述三维体数据是否收敛；

当所述三维体数据收敛时，将所述剩余堆图像、所述模板堆图像分别与所述三维体数据局部配准，否则跳转至通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新的步骤；

判断局部配准后的所述三维体数据是否收敛，是则输出所述三维体数据，否则跳转至通过预设的基于点扩展函数的体重建方式构建三维体数据的步骤。

另一方面，本发明提供了一种三维超声重建装置，所述装置包括：

全局配准单元，用于当接收多角度超声扫描到的堆图像时，将所有堆图像中运动误差最小的堆图像设置为模板堆图像，并将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准；

三维重建单元，用于根据所述剩余堆图像和所述模板堆图像，通过预设的基于点扩展函数的体重建方式构建三维体数据；

三维体更新单元，用于通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新，判断更新后的所述三维体数据是否收敛；

局部配准单元，用于当所述三维体数据收敛时，将所述剩余堆图像、所述模板堆图像分别与所述三维体数据局部配准，否则由所述三维体更新单元执行通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新的操作；以及

体数据输出单元，用于判断局部配准后的所述三维体数据是否收敛，是则输出所述三维体数据，否则由所述三维重建单元执行通过预设的基于点扩展函数的体重建方式构建三维体数据的操作。

另一方面，本发明还提供了一种图像处理设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上述一种三维超声重建方法所述的步骤。

另一方面，本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如上述一种三维超声重建方法所述的步骤。

本发明在多角度超声扫描到的堆图像中选取运动误差最小的堆图像作为模板堆图像，将剩余堆图像与模板堆图像全局配准，根据配准后的剩余堆图像和模板堆图像，通过基于点扩展函数的体重建方式构建三维体数据，再通过核回归函数的体重建方式对三维体数据进行更新，判断更新后的三维体数据是否收敛，是则将剩余堆图像、模板堆图像分别与三维体数据局部配准，否则继续对三维体数据进行更新，判断局部配准后的三维体数据是否收敛，输出三维体数据，否则继续跳转至通过基于点扩展函数的体重建方式构建三维体数据的步骤，从而实现了多角度扫描的三维超声重建，解决了单一角度扫描给医生带来的观察限制，有效地去除了三维超声图像中的伪影，提高了三维超声图像的重建质量。

附图说明

图1是本发明实施例一提供的一种三维超声重建方法的实现流程图；

图2是本发明实施例二提供的一种三维超声重建装置的结构示意图；

图3是本发明实施例二提供的一种三维超声重建装置的优选结构示意图；以及

图4是本发明实施例三提供的图像处理设备的结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

以下结合具体实施例对本发明的具体实现进行详细描述：

实施例一：

图1示出了本发明实施例一提供的一种三维超声重建方法的实现流程，为了便于说明，仅示出了与本发明实施例相关的部分，详述如下：

在步骤S101中，当接收多角度超声扫描到的堆图像时，将所有堆图像中运动误差最小的堆图像设置为模板堆图像，并将所有堆图像中的剩余堆图像与模板堆图像全局配准。

本发明实施例适用于三维超声***。在不同角度进行超声扫描，获得不同角度扫描得到的二维超声图像序列(即二维切片序列)，同一角度的二维超声图像序列可称为一个堆图像。可采用基于运动量估计的方法在所有堆图像中选取一个运动误差最小的堆图像作为模板堆图像，再将所有堆图像中的剩余堆图像与模板堆图像进行全局配准，以使每个剩余堆图像与模板堆图像对齐，降低超声图像重建过程中的伪影。作为示例地，不同角度可包括轴向、冠状和矢状方向。

在本发明实施例中，在所有堆图像中选取一个运动误差最小的堆图像作为模板堆图像的过程可通过下述步骤实现：

(1)将每个二维超声图像向量化，得到每个堆图像对应的矩阵。

在本发明实施例中，每个堆图像对应的矩阵可表示为A＝[vec(I₁)；...；vec(I_k)]∈R^m*k，其中，vec(I_j)表示将堆图像中的第j(j∈1，...，k)个二维超声图像I_j作为一个列向量。

(2)根据每个堆图像对应的矩阵生成每个堆图像对应的观察矩阵，对每个堆图像的观察矩阵进行奇异值分解。

在本发明实施例中，堆图像的观察矩阵为D＝A+E，E为单位矩阵，D是满秩的且D∈R^m*k，所以D的奇异值的递减顺序为s₁≥s₂≥...≥s_k≥0，可对D进行奇异值分解得到D＝USV^T(U∈R^m*k，S∈R^k*k，V∈R^k*k)，S且是对角线上元素为奇异值的对角阵。

(3)根据观察矩阵的奇异值分解结果，生成与观察矩阵相似的矩阵，根据观察矩阵和与观察矩阵相似的矩阵，确定所有堆图像中的模板堆图像。

在本发明实施例中，可分别从U_m*k、S_k*k、V_k*k矩阵的左上角取r列子矩阵U′_m*r、S′_r*r、U′_r*k，可得D′＝U′S′V′^T，D和D′是非常相似的，可采用弗罗贝尼乌斯范数(Frobenius norm)来衡量D和D′之间的相似程度：

其中，||D||为D的弗罗贝尼乌斯范数，||D′||为D′的弗罗贝尼乌斯范数，||D-D′||用来表示D和D′之间的相似程度。所以，D和D′之间的相对相似差异可以表示为：

在满足δ_r小于预设阈值的条件下，每个堆图像都可以得出对应的使得δ_r最小的r值，然后在所有堆图像中选择具有最小的r值的堆图像作为模板堆图像。

优选地，为了更加直观和方便，使用一个用于衡量运动估计量的变量ω来简化上述过程，ω＝r·δ_r，求出每个堆图像分别对应的ω的最小值，将具有最小ω的堆图像设置为模板堆图像。

在本发明实施例中，在将剩余堆图像与模板堆图像全局配准时，可先采用刚性配准方式初步计算出变换矩阵，即在三维空间中利用平移参数和旋转参数计算变换矩阵，为了便于与后续的局部配准区分，将该变换矩阵称为全局变换矩阵，再通过全局变换矩阵对剩余堆图像进行配准，具体地，通过以下公式对剩余堆图像S_i进行全局配准：

S′_i＝T_global·S_i，i≠T，其中，T_global为全局变换矩阵，S_T表示模板堆图像，当i≠T时，S_i表示第i个剩余堆图像，S′_i表示配准后的剩余堆图像。在配准后，计算剩余堆图像与模板堆图像之间的测度，判断该测度是否满足预设的第一测度阈值，是则完成堆图像的配准，否则通过预设的优化方法对全局变换矩阵进行优化，通过优化后的全局变换矩阵对堆图像进行配准。其中，测度用来定量地衡量两个堆图像之间的匹配效果，即描述剩余堆图像与模板堆图像之间的差异性或者相似性，例如，差异性测度由均方差差值(Mean squareddifference)、熵之间的差(Entropy of difference)等等，相似性测度有归一化相互关系数(Normalised correlation coefficient)、互信息(Mutual information)等等，超声图像配准即通过差异性测度函数不断迭代求最小值的过程、或者通过相似性测度函数不断迭代求最大值的过程，优化方式可采用梯度下降法、遗传算法或模拟退火法等全局优化算法。

在步骤S102中，根据剩余堆图像和模板堆图像，通过预设的基于点扩展函数的体重建方式构建三维体数据。

在本发明实施例中，构建三维体数据时，需要将剩余堆图像、模板堆图像中二维超声图像上的像素转换为三维体数据上的体素，而二维超声图像上的像素与三维体数据上的体素不会很好的对齐，二维超声图像上的一个像素点会影响到三维体数据上多个体素点的值，为了更准确地重建出三维体数据，可采用基于点扩展函数的体重建方法，每次重建过程即利用二维超声图像上的像素更新三维体数据上体素值的过程。在剩余堆图像与模板堆图像全局配准后，采用基于点扩展函数的体重建方法对三维体数据进行初次重构，后续需要对初次重构的三维体数据进行更新、局部配准等操作，若经过这些操作得到的三维体数据不收敛，则再次对三维体数据进行重构(即对三维体数据进行更新)、更新、局部配准等操作，直至最后的三维体数据收敛。

在本发明实施例中，从像素点p_s到体素点p_r的空间转换过程可表示为：

其中，W_s为将堆图像从图像坐标系变换到世界坐标系的变换矩阵，W_r为将三维重建体从目标坐标系变换到世界坐标系的变换矩阵，T_total为配准矩阵，在初次重建时由于尚未进行局部配准所以T_total＝T_global，在后续多次重建过程中T_total＝T_global·T_local，T_local为用于局部配准的局部转换矩阵，令则p_s＝F·p_r。因此，可以根据该空间转换过程将剩余堆图像、模板堆图像中二维超声图像上的像素转换到三维体数据上，初步得到一个三维体数据，之后便执行步骤S103。

在后续第n次采用基于点扩展函数的体重建方法对三维体数据进行重构(或更新)时，体素值的更新公式可表示为：

为三维体上原有体素点经过空间转换得到的像素点，PSF为点扩展函数。其中，点扩展函数可选用三维高斯点扩展函数：

在步骤S103中，通过预设的基于核回归函数的体重建方式对三维体数据进行更新。

在本发明实施例中，由于三维超声的采样数据通常是稀疏的，三维体数据在重建或更新后，依然会存在一些未分配像素的空白区域，因此采用基于核回归函数的体重建方式对三维体数据上的每个体素点进行更新。

在本发明实施例中，采用基于核回归函数的体重建方式的数学观察模型为：

Y_i＝r(X_i)+ε_i，其中，i＝1，...，M，M为三维体数据上所有体素点的个数，r(·)为核回归函数，X_i＝(X_i0，X_i1，X_i2)为三维体数据X_i上体素点的三维坐标，r(X_i)得到的是体素点X_i的理想体素值，Y_i为是体素点X_i的实际体素值，ε_i为均差为0、方差为σ²的高斯噪声。依次将三维体数据上的体素点设置为待更新体素点，当当前待更新体素点X距离体素点X_i很近时，根据N阶泰勒公式可求得体素点X_i的理想体素值为：

根据最小化二乘法、上述N级泰勒公式和数学观察模型，可得优化问题：

其中，m为待求体素点X邻域(例如3*3*3大小的邻域，在此不对邻域的大小进行限制)内体素点X_i的个数，h为预设的平滑参数，K(·)为预设的核函数，核函数K(·)通常为指数函数或高斯函数，且需要满足∫tk(t)dt＝0、∫t²K(t)dt＝c的条件。上述优化问题可通过矩阵进行表示：

其中，Y为所有体素点X_i的体素值组成的向量集，Y＝[Y₁，Y₂，...，Y_m]^T，为所有β_i组成的向量集，，W＝diag[K(X₀-X)，K(X₁-X)，...，K(X_m-X)]为对角矩阵且对角上的元素为K(·)、其它元素为零，根据最小二乘法对上述通过矩阵表示的优化问题进行求解，可得：

其中， 为体素点X的体素估计值，即最终用来对待更新体素点X进行更新的体素值。

在步骤S104中，判断更新后的三维体数据是否收敛。

在本发明实施例中，可预先设置通过基于核回归函数的体重建方式对三维体数据上所有体素点进行更新的最大更新次数，在对三维体素点上所有体素进行一次更新后，通过判断当前更新次数是否超过该最大更新次数，来确定更新后的三维体数据是否收敛，是则执行步骤S105，否则将当前更新次数加一，并跳转至步骤S103继续对三维体数据进行更新。此外，也可通过更新前后两个三维体数据之间的误差是否超过预设的更新误差阈值，来确定更新后的三维体数据是否收敛。

在步骤S105中，将剩余堆图像、模板堆图像分别与三维体数据局部配准。

在本发明实施例中，将剩余堆图像、模板堆图像中的二维超声图像分别与三维体数据局部配准，以调整每个二维超声图像与三维体数据之间的空间对应关系。可先采用非刚性配准方式计算出局部变换矩阵T_local，通过局部变换矩阵对二维超声图像进行配准：

其中，j＝1，...，k，I_j表示剩余堆图像或模板堆图像中第j个二维超声图像，I′_j为对I_j进行配准后得到的二维超声图像，k表示剩余堆图像或模板堆图像中二维超声图像的数量。在配准后计算二维超声图像与三维体数据之间的测度，判断该测度是否满足预设的第二测度阈值，是则完成二维超声图像的配准，否则通过预设的优化方法对局部变换矩阵进行优化，继续通过优化后的全局变换矩阵对二维超声图像进行配准。其中，测度可为差异性测度或相似性测度，优化方式可采用梯度下降法、遗传算法或模拟退火法等全局优化算法。

在步骤S106中，判断局部配准后的三维体数据是否收敛。

在本发明实施例中，可预先设置通过基于点扩展函数的体重建方式对三维体数据构建的最大构建次数，通过判断当前构建次数是否超过该最大构建次数，来确定三维体数据是否收敛，是则执行步骤S107，否则将当前构建次数加一，并跳转至步骤S102继续通过预设的基于点扩展函数的体重建方式构建三维体数据。此外，还可通过判断前后两次构建得到的三维体数据之间的误差是否小于预设重建误差阈值，来确定三维重建体数据是否收敛。

在步骤S107中，输出三维体数据。

在本发明实施例中，在多角度超声扫查到的堆图像中选取运动误差最小的堆图像作为模板堆图像，将剩余堆图像与模板堆图像进行全局配准，根据全局配准后的剩余堆图像和模板堆图像构建三维重建数据，采用核回归方式对三维体数据进行多次更新直至三维体数据收敛，将剩余堆图像、模板堆图像中的二维超声图像与三维体数据进行局部配准，检测局部配准后的三维体数据是否收敛，是则输出三维体数据，否则继续构建三维体数据，从而实现了多角度扫描的三维超声重建，解决了单一角度扫描给医生带来的观察限制，有效地去除了三维超声图像中的伪影，提高了三维超声图像的重建质量。

实施例二：

图2示出了本发明实施例二提供的一种三维超声重建装置的结构，为了便于说明，仅示出了与本发明实施例相关的部分，其中包括：

全局配准单元21，用于当接收多角度超声扫描到的堆图像时，将所有堆图像中运动误差最小的堆图像设置为模板堆图像，并将所有堆图像中的剩余堆图像与模板堆图像全局配准。

在本发明实施例中，在不同角度进行超声扫描，获得不同角度扫描得到的堆图像。可采用基于运动量估计的方法在所有堆图像中选取一个运动误差最小的堆图像作为模板堆图像，再将所有堆图像中的剩余堆图像与模板堆图像进行全局配准，以使每个剩余堆图像与模板堆图像对齐，降低超声图像重建过程中的伪影。

在本发明实施例中，在所有堆图像中选取一个运动误差最小的堆图像作为模板堆图像的过程可通过下述步骤实现：

(1)将每个二维超声图像向量化，得到每个堆图像对应的矩阵。

在本发明实施例中，每个堆图像对应的矩阵可表示为A＝[vec(I₁)；...；vec(I_k)]∈R^m*k，其中，vec(I_j)表示将堆图像中的第j(j∈1，...，k)个二维超声图像I_j作为一个列向量。

(2)根据每个堆图像对应的矩阵生成每个堆图像对应的观察矩阵，对每个堆图像的观察矩阵进行奇异值分解。

在本发明实施例中，堆图像的观察矩阵为D＝A+E，E为单位矩阵，D是满秩的且D∈R^m*k，所以D的奇异值的递减顺序为s₁≥s₂≥...≥s_k≥0，可对D进行奇异值分解得到D＝USV^T(U∈R^m*k，S∈R^k*k，V∈R^k*k)，S且是对角线上元素为奇异值的对角阵。

(3)根据观察矩阵的奇异值分解结果，生成与观察矩阵相似的矩阵，根据观察矩阵和与观察矩阵相似的矩阵，确定所有堆图像中的模板堆图像。

在本发明实施例中，可分别从U_m*k、S_k*k、V_k*k矩阵的左上角取r列子矩阵U′_m*r、S′_r*r、U′_r*k，可得D′＝U′S′V′^T，D和D′是非常相似的，可采用弗罗贝尼乌斯范数(Frobenius norm)来衡量D和D′之间的相似程度：

其中，||D||为D的弗罗贝尼乌斯范数，||D′||为D′的弗罗贝尼乌斯范数，||D-D′||用来表示D和D′之间的相似程度。所以，D和D′之间的相对相似差异可以表示为：

在满足δ_r小于预设阈值的条件下，每个堆图像都可以得出对应的使得δ_r最小的r值，然后在所有堆图像中选择具有最小的r值的堆图像作为模板堆图像。

优选地，为了更加直观和方便，使用一个用于衡量运动估计量的变量ω来简化上述过程，ω＝r·δ_r，求出每个堆图像分别对应的ω的最小值，将具有最小ω的堆图像设置为模板堆图像。

在本发明实施例中，在将剩余堆图像与模板堆图像全局配准时，可先采用刚性配准方式初步计算出全局变换矩阵，再通过全局变换矩阵对剩余堆图像进行配准，具体地，通过以下公式对剩余堆图像S_i进行全局配准：

S′_i＝T_global·S_i，i≠T，其中，T_global为全局变换矩阵，S_T表示模板堆图像，当i≠T时，S_i表示第i个剩余堆图像，S′_i表示配准后的剩余堆图像。在配准后，计算剩余堆图像与模板堆图像之间的测度，判断该测度是否满足预设的第一测度阈值，是则完成堆图像的配准，否则通过预设的优化方法对全局变换矩阵进行优化，通过优化后的全局变换矩阵对堆图像进行配准。

三维重建单元22，用于根据剩余堆图像和模板堆图像，通过预设的基于点扩展函数的体重建方式构建三维体数据。

在本发明实施例中，构建三维体数据时，需要将剩余堆图像、模板堆图像中二维超声图像上的像素转换为三维体数据上的体素，而二维超声图像上的像素点与三维体数据上的体素不会很好的对齐，二维超声图像上的一个像素点会影响到三维体数据上多个体素点的值，为了更准确地重建出三维体数据，可采用基于点扩展函数的体重建方法，每次重建过程即利用二维超声图像上的像素更新三维体数据上体素值的过程。在剩余堆图像与模板堆图像全局配准后，采用基于点扩展函数的体重建方法对三维体数据进行初次重构，后续需要对初次重构的三维体数据进行更新、局部配准等操作，若经过这些操作得到的三维体数据不收敛，则再次对三维体数据进行重构(即对三维体数据进行更新)、更新、局部配准等操作，直至最后的三维体数据收敛。

在本发明实施例中，从像素点p_s到体素点p_r的空间转换过程可表示为：

其中，W_s为将堆图像从图像坐标系变换到世界坐标系的变换矩阵，W_r为将三维重建体从目标坐标系变换到世界坐标系的变换矩阵，T_total为配准矩阵，在初次重建时由于尚未进行局部配准所以T_total＝T_global，在后续多次重建过程中T_total＝T_global·T_local，T_local为用于局部配准的局部转换矩阵，令则p_s＝F·p_r。因此，可以根据该空间转换过程将剩余堆图像、模板堆图像中二维超声图像上的像素转换到三维体数据上，初步得到一个三维体数据，之后便由三维体更新单元23执行对三维体数据进行更新的操作。

在后续第n次采用基于点扩展函数的体重建方法对三维体数据进行重构(或更新)时，体素值的更新公式可表示为：

为三维体上原有体素点经过空间转换得到的像素点，PSF为点扩展函数。其中，点扩展函数可选用三维高斯点扩展函数：

三维体更新单元23，用于通过预设的基于核回归函数的体重建方式对三维体数据进行更新，判断更新后的三维体数据是否收敛。

在本发明实施例中，由于三维超声的采样数据通常是稀疏的，三维体数据在重建或更新后，依然会存在一些未分配像素的空白区域，因此采用基于核回归函数的体重建方式对三维体数据上的每个体素点进行更新。

在本发明实施例中，采用基于核回归函数的体重建方式的数学观察模型为：

Y_i＝r(X_i)+ε_i，其中，i＝1，...，M，M为三维体数据上所有体素点的个数，r(·)为核回归函数X_i＝(X_i0，X_i1，X_i2)为三维体数据上体素点的三维坐标，r(X_i)得到的是体素点X_i的理想体素值，Y_i为是体素点X_i的实际体素值，ε_i为均差为0、方差为σ²的高斯噪声。依次将三维体数据上的体素点设置为待更新体素点，当当前待更新体素点X距离体素点X_i很近时，根据N阶泰勒公式可求得体素点X_i的理想体素值为：

根据最小化二乘法、上述N级泰勒公式和数学观察模型，可得优化问题：

其中，m为待求体素点X的邻域内体素点X_i的个数，h为预设的平滑参数，K(·)为预设的核函数，核函数K(·)通常为指数函数或高斯函数，且需要满足∫tk(t)dt＝0、∫t²K(t)dt＝c的条件。上述优化问题可通过矩阵进行表示：

其中，Y为所有体素点X_i的体素值组成的向量集，Y＝[Y₁，Y₂，...，Y_m]^T，为所有β_i组成的向量集，，W＝diag[K(X₀-X)，K(X₁-X)，...，K(X_m-X)]为对角矩阵且对角上的元素为K(·)、其它元素为零，根据最小二乘法对上述通过矩阵表示的优化问题进行求解，可得：

其中， 为体素点X的体素估计值，即最终用来对待更新体素点X进行更新的体素值。

在本发明实施例中，可预先设置通过基于核回归函数的体重建方式对三维体数据上所有体素点进行更新的最大更新次数，在对三维体素点上所有体素进行一次更新后，通过判断当前更新次数是否超过该最大更新次数，来确定更新后的三维体数据是否收敛。此外，也可通过更新前后两个三维体数据之间的误差是否超过预设的更新误差阈值，来确定更新后的三维体数据是否收敛。

局部配准单元24，用于当三维体数据收敛时，将剩余堆图像、模板堆图像分别与三维体数据局部配准，否则由三维体更新单元23执行通过预设的基于核回归函数的体重建方式对三维体数据进行更新的操作。

在本发明实施例中，当更新后的三维体数据收敛时，将剩余堆图像、模板堆图像中的二维超声图像分别与三维体数据局部配准，以调整每个二维超声图像与三维体数据之间的空间对应关系，否则将当前更新次数加一，并由三维体更新单元23执行通过预设的基于核回归函数的体重建方式对三维体数据进行更新的操作。

在本发明实施例中，在将剩余堆图像、模板堆图像中的二维超声图像分别与三维体数据局部配准时，可先采用非刚性配准方式计算出局部变换矩阵T_local，通过局部变换矩阵对二维超声图像进行配准：

其中，j＝1，...，k，I_j表示剩余堆图像或模板堆图像中第j个二维超声图像，I′_j为对I_j进行配准后得到的二维超声图像，k表示剩余堆图像或模板堆图像中二维超声图像的数量。在配准后计算二维超声图像与三维体数据之间的测度，判断该测度是否满足预设的第二测度阈值，是则完成二维超声图像的配准，否则通过预设的优化方法对局部变换矩阵进行优化，继续通过优化后的全局变换矩阵对二维超声图像进行配准。

体数据输出单元25，用于判断局部配准后的三维体数据是否收敛，是则输出三维体数据，否则由三维重建单元22执行通过预设的基于点扩展函数的体重建方式构建三维体数据的操作。

在本发明实施例中，可预先设置通过基于点扩展函数的体重建方式对三维体数据构建的最大构建次数，通过判断当前构建次数是否超过该最大构建次数，来确定三维体数据是否收敛，是则输出三维体数据，否则将当前构建次数加一，并由三维重建单元22执行通过预设的基于点扩展函数的体重建方式构建三维体数据的操作。此外，还可通过判断前后两次构建得到的三维体数据之间的误差是否小于预设重建误差阈值，来确定三维重建体数据是否收敛。

优选地，如图3所示，全局配准单元21包括：

模板确定单元311，用于通过预设的弗罗贝尼乌斯范数计算每个堆图像对应的运动估计量，将所有堆图像中运动估计量最小的堆图像设置为模板堆图像；

全局配准子单元312，用于根据预设的全局变换矩阵将所有堆图像中的剩余堆图像与模板堆图像全局配准，计算全局配准后剩余堆图像与模板堆图像之间的测度；以及

全局配准判断单元313，用于判断剩余堆图像与模板堆图像之间的测度是否满足预设第一测度阈值，是则结束剩余堆图像与模板堆图像的全局配准，否则对全局变换矩阵进行优化，由全局配准子单元312执行将所有堆图像中的剩余堆图像与模板堆图像全局配准的操作。

优选地，局部配准单元24包括：

局部配准子单元341，用于根据预设的局部变换矩阵，将剩余堆图像、模板堆图像中的二维超声图像与三维体数据局部配准，计算二维超声图像与三维体数据之间的测度；以及

局部配准判断单元342，用于判断二维超声图像与三维体数据之间的测度是否满足预设第二测度阈值，是则结束二维超声图像与三维体数据的局部配准，否则对局部变换矩阵进行优化，并由局部配准子单元341执行将剩余堆图像、模板堆图像中的二维超声图像与三维体数据局部配准的操作。

在本发明实施例中，在多角度超声扫查到的堆图像中选取运动误差最小的堆图像作为模板堆图像，将剩余堆图像与模板堆图像进行全局配准，根据全局配准后的剩余堆图像和模板堆图像构建三维重建数据，采用核回归方式对三维体数据进行多次更新直至三维体数据收敛，将剩余堆图像、模板堆图像中的二维超声图像与三维体数据进行局部配准，检测局部配准后的三维体数据是否收敛，是则输出三维体数据，否则继续构建三维体数据，从而实现了多角度扫描的三维超声重建，解决了单一角度扫描给医生带来的观察限制，有效地去除了三维超声图像中的伪影，提高了三维超声图像的重建质量。

在本发明实施例中，一种三维超声重建装置的各单元可由相应的硬件或软件单元实现，各单元可以为独立的软、硬件单元，也可以集成为一个软、硬件单元，在此不用以限制本发明。

实施例三：

图4示出了本发明实施例三提供的图像处理设备的结构，为了便于说明，仅示出了与本发明实施例相关的部分。

本发明实施例的图像处理设备4包括处理器40、存储器41以及存储在存储器41中并可在处理器40上运行的计算机程序42。该处理器40执行计算机程序42时实现上述方法实施例中的步骤，例如图1所示的步骤S101至S107。或者，处理器40执行计算机程序42时实现上述装置实施例中各单元的功能，例如图2所示单元21至25的功能。

在本发明实施例中，在多角度超声扫查到的堆图像中选取运动误差最小的堆图像作为模板堆图像，将剩余堆图像与模板堆图像进行全局配准，根据全局配准后的剩余堆图像和模板堆图像构建三维重建数据，采用核回归方式对三维体数据进行多次更新直至三维体数据收敛，将剩余堆图像、模板堆图像中的二维超声图像与三维体数据进行局部配准，检测局部配准后的三维体数据是否收敛，是则输出三维体数据，否则继续构建三维体数据，从而实现了多角度扫描的三维超声重建，解决了单一角度扫描给医生带来的观察限制，有效地去除了三维超声图像中的伪影，提高了三维超声图像的重建质量。

实施例四：

在本发明实施例中，提供了一种计算机可读存储介质，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现上述方法实施例中的步骤，例如，图1所示的步骤S101至S107。或者，该计算机程序被处理器执行时实现上述装置实施例中各单元的功能，例如图2所示单元21至25的功能。

在本发明实施例中，在多角度超声扫查到的堆图像中选取运动误差最小的堆图像作为模板堆图像，将剩余堆图像与模板堆图像进行全局配准，根据全局配准后的剩余堆图像和模板堆图像构建三维重建数据，采用核回归方式对三维体数据进行多次更新直至三维体数据收敛，将剩余堆图像、模板堆图像中的二维超声图像与三维体数据进行局部配准，检测局部配准后的三维体数据是否收敛，是则输出三维体数据，否则继续构建三维体数据，从而实现了多角度扫描的三维超声重建，解决了单一角度扫描给医生带来的观察限制，有效地去除了三维超声图像中的伪影，提高了三维超声图像的重建质量。

本发明实施例的计算机可读存储介质可以包括能够携带计算机程序代码的任何实体或装置、记录介质，例如，ROM/RAM、磁盘、光盘、闪存等存储器。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种三维超声重建方法，其特征在于，所述方法包括下述步骤：

当接收多角度超声扫描到的堆图像时，将所有堆图像中运动误差最小的堆图像设置为模板堆图像，并将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准；

根据所述剩余堆图像和所述模板堆图像，通过预设的基于点扩展函数的体重建方式构建三维体数据；

通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新，判断更新后的所述三维体数据是否收敛；

当所述三维体数据收敛时，将所述剩余堆图像、所述模板堆图像分别与所述三维体数据局部配准，否则跳转至通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新的步骤；

判断局部配准后的所述三维体数据是否收敛，是则输出所述三维体数据，否则跳转至通过预设的基于点扩展函数的体重建方式构建三维体数据的步骤。

2.如权利要求1所述的方法，其特征在于，将所有堆图像中运动误差最小的堆图像设置为模板堆图像，并将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准的步骤，包括：

通过预设的弗罗贝尼乌斯范数计算所述每个堆图像对应的运动估计量，将所述所有堆图像中所述运动估计量最小的所述堆图像设置为模板堆图像；

根据预设的全局变换矩阵将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准，计算全局配准后所述剩余堆图像与所述模板堆图像之间的测度；

判断所述剩余堆图像与所述模板堆图像之间的测度是否满足预设第一测度阈值，是则结束所述剩余堆图像与所述模板堆图像的全局配准，否则对所述全局变换矩阵进行优化，跳转至根据预设的全局变换矩阵将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准的步骤。

3.如权利要求1所述的方法，其特征在于，将所述剩余堆图像、所述模板堆图像分别与所述三维体数据局部配准的步骤，包括：

根据预设的局部变换矩阵，将所述剩余堆图像、所述模板堆图像中的二维超声图像与所述三维体数据局部配准，计算所述二维超声图像与所述三维体数据之间的测度；

判断所述二维超声图像与所述三维体数据之间的测度是否满足预设第二测度阈值，是则结束所述二维超声图像与所述三维体数据的局部配准，否则对所述局部变换矩阵进行优化，并跳转至根据预设的局部变换矩阵，将所述剩余堆图像、所述模板堆图像中的二维超声图像与所述三维体数据局部配准的步骤。

4.如权利要求1所述的方法，其特征在于，通过预设的基于点扩展函数的体重建方式构建三维体数据的步骤，包括：

根据所述剩余堆图像、所述模板堆图像中二维超声图像上的像素值、以及所述基于点扩展函数的体重建方式，计算并更新所述三维体数据上的体素值，所述体素值的计算公式为：

其中，所述V为所述三维体数据，所述为第n迭代更新时计算得到的所述三维体数据在体素点p_r处的体素值，所述为所述二维超声图像上在像素点p_s处的像素值，所述为所述体素p_r ⁿ空间转换后得到的像素点，所述PSF为所述点扩展函数。

5.如权利要求1所述的方法，其特征在于，通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新的步骤，包括：

依次将所述三维体数据上的每个体素点设置为待更新体素点，并根据所述待更新体素点邻域内的体素点、所述邻域内所有体素点的体素值、以及所述基于核回归函数的体重建方式，计算所述待更新体素点对应的体素估计值，将所述体素估计值设置为所述待更新体素点的体素值，所述体素估计值的计算公式为：

其中，所述为所述待更新体素点X的体素估计值，所述X_i为所述X邻域内的体素点，所述i＝1,...,m，所述m为所述X邻域内体素点的个数，所述Y为所述所有X_i的体素值构成的向量集，即Y＝[Y₁,Y₂,...,Y_m]^T，所述W＝diag[K(X₀-X),K(X₁-X),...,K(X_m-X)]，所述K(·)为预设的核函数，所述

6.一种三维超声重建装置，其特征在于，所述装置包括：

全局配准单元，用于当接收多角度超声扫描到的堆图像时，将所有堆图像中运动误差最小的堆图像设置为模板堆图像，并将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准；

三维重建单元，用于根据所述剩余堆图像和所述模板堆图像，通过预设的基于点扩展函数的体重建方式构建三维体数据；

三维体更新单元，用于通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新，判断更新后的所述三维体数据是否收敛；

局部配准单元，用于当所述三维体数据收敛时，将所述剩余堆图像、所述模板堆图像分别与所述三维体数据局部配准，否则由所述三维体更新单元执行通过预设的基于核回归函数的体重建方式对所述三维体数据进行更新的操作；以及

体数据输出单元，用于判断局部配准后的所述三维体数据是否收敛，是则输出所述三维体数据，否则由所述三维重建单元执行通过预设的基于点扩展函数的体重建方式构建三维体数据操作。

7.如权利要求6所述的装置，其特征在于，所述全局配准单元包括：

模板确定单元，用于通过预设的弗罗贝尼乌斯范数计算所述每个堆图像对应的运动估计量，将所述所有堆图像中所述运动估计量最小的所述堆图像设置为模板堆图像；

全局配准子单元，用于根据预设的全局变换矩阵将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准，计算全局配准后所述剩余堆图像与所述模板堆图像之间的测度；以及

全局配准判断单元，用于判断所述剩余堆图像与所述模板堆图像之间的测度是否满足预设第一测度阈值，是则结束所述剩余堆图像与所述模板堆图像的全局配准，否则对所述全局变换矩阵进行优化，由所述全局配准子单元执行将所述所有堆图像中的剩余堆图像与所述模板堆图像全局配准的操作。

8.如权利要求6所述的装置，其特征在于，所述局部配准单元包括：

局部配准子单元，用于根据预设的局部变换矩阵，将所述剩余堆图像、所述模板堆图像中的二维超声图像与所述三维体数据局部配准，计算所述二维超声图像与所述三维体数据之间的测度；以及

局部配准判断单元，用于判断所述二维超声图像与所述三维体数据之间的测度是否满足预设第二测度阈值，是则结束所述二维超声图像与所述三维体数据的局部配准，否则对所述局部变换矩阵进行优化，并由所述局部配准子单元执行将所述剩余堆图像、所述模板堆图像中的二维超声图像与所述三维体数据局部配准的操作。

9.一种图像处理设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至5任一项所述方法的步骤。

10.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1至5任一项所述方法的步骤。