一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质技术

本发明专利技术提供一种基于三维空间体拟合的医学图像层间插值方法、设备及可读存储介质，配置三维空间体素拟合及配准的医学图像层间插值算法的两个约束条件；构建三维空间体；修正空间体；构造逼近空间体，对空间体进行重采样生成待插值医学图像。对空间体进行重采样生成质量较好的待插值层间图像.通过采样公式对原图像序列反向采样局部地构造三元二次多项式空间体,将二次多项式空间体加权平均生成单位空间体,再将所有的单位空间体拼合成拟合三维空间体,最后对空间体进行重采样生成待插值图像.方法生成的层间图像的质量不论在视觉评价还是客观数值上均有所提升,且有效地提高了插值效率。

【技术实现步骤摘要】
一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质
本专利技术涉及医学图像处理
，尤其涉及一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质。
技术介绍
断层医学图像层间插值技术是医学图像超分辨重建和三维重建的关键技术之一。即,在相邻两层CT切片间通过层间插值技术生成新的切片,增加切片数量,提高序列的层间分辨率。如图1所示。层间插值技术在医学图像处理中有广泛的应用和重要的研究意义,各种新技术层出不穷。近年来的方法可以分为三种:基于灰度的插值方法、基于形状的插值方法以及基于配准的插值方法。1)基于灰度的插值方法直接利用原始序列中上下两层图像的灰度信息通过一组基函数内插层间图像。最近邻插值、线性插值、三次B样条函数插值等是此类插值方法常用的类型。在层间距较小的情况下,基于灰度的插值方法计算简单、时间复杂度低且易于实现,应用较广。然而,对于相邻切片间组织结构形变较大的CT图像序列,该方法得到的插值图像通常过于平滑且包含伪影,不能准确反映组织轮廓的自然渐变状态。随后,Goshtasby提出匹配插值算法,在一定程度上改进了边缘模糊问题,提高了插值图像的精度,但是时间复杂度高。2)基于形状的插值方法是通过上下两层原始切片的形状,直接生成待插值切片的轮廓,相较于基于灰度的插值方法它能够有效消除伪影,提高插值图像质量。但是,基于形状的方法依然具有很大的局限性。一方面,其对组织轮廓的质量有很高的要求,适合上下两层切片十分相似的情况。如果层间距过大,或连续两层CT图像存在较大的结构差异,都会严重影响该类方法的效果,导致插值结果出现较大误差。另一方面,基于形状的方法时间复杂度高于基于灰度的插值方法。3)基于配准的插值方法在医学图像领域应用较为广泛,该方法通过找到相邻切片结构之间的对应关系来匹配局部解剖结构,根据找到的形变场获得待插值图像像素点的形变信息,进行图像层间插值。基于配准的插值能在消除边界伪影的同时,考虑到人体器官组织结构的各种非刚性形变，可以产生视觉上不错的结果。目前流行的非刚性配准方法可以分为2类：基于空间变换的配准，例如Rueckert等提出基于B样条的自由形变配准方法；基于物理模型的配准，例如基于光流估计的配准方法。但是，这些基于配准的插值方法在配准过程中存在近似计算，导致原图像中多个像素可能映射到插值图像中的同一位置，使得在插值图像中存在没有原图像像素映射到的区域，也就是像素缺失，出现异常点。相较于基于形状、基于灰度的插值，其插值计算复杂度更高。
技术实现思路
为了克服上述现有技术中的不足，本专利技术提供一种基于三维空间体拟合的医学图像层间插值方法，方法包括：步骤一，配置三维空间体素拟合及配准的医学图像层间插值算法的两个约束条件；步骤二，构建三维空间体fi，j，k(x，y，z)；步骤三，修正空间体fi，j，k(u，v，w)；步骤四，构造逼近空间体f(x，y，z)，对空间体进行重采样生成待插值医学图像。进一步需要说明的是，步骤一还包括：设相邻两个断层切片Sk，Sk+1之间的距离为1，待插值层间切片Sk+d(0＜d＜1)由m×n个体素Vi，j，k+di＝1，2，...，m，j＝1，2，...，n组成，待插值体素配置为原空间体F(x，y，z)上的采样值；设每个体素Vi，j，k+d是从单位体积上采样得到，即其中，w(x，y，z)是权函数；F(x，y，z)定义在[0.5，m+0.5]×[0.5，m+0.5]×[0.5，t+0.5]Vi，j，k+d(i＝1，2，...，m；j＝1，2，...，n；k＝1，2，...，t)上，且对于每一个都满足式(1)。同时，将权函数w(x，y，z)设置为1。进一步需要说明的是，步骤二还包括：在三维空间中，空间体fi，j，k(x，y，z)，i＝2，3，...，m-1，j＝2，3，...，n-1，k＝2，3，...，t-1满足以下条件：如果式(1)中的F(x，y，z)是三元二次多项式，则fi，j，k(x，y，z)重建F(x，y，z)，即fi，j，k(x，y，z)＝F(x，y，z)；令u＝x-i，v＝y-j，w＝z-k，则fi，j，k(x，y，z)在uvw空间的[-1.5，1.5]×[-1.5，1.5]×[-1.5，1.5]三维空间内写成：若F(x，y，z)由式(2)定义，Pi，j，k由式(1)定义，则：其中，e1＝(Pi+1，j，k-Pi-1，j，k)/2，e2＝(Pi，j+1，k-Pi，j-1，k)/2，e3＝(Pi，j，k+1-Pi，j，k-1)/2，e4＝(Pi+1，j+1，k-Pi-1，j-1，k)/2，e5＝(Pi+1，j-1，k-Pi-1，j+1，k)/2，e6＝(Pi+1，j，k+1-Pi-1，j，k-1)/2，e7＝(Pi+1，j，k-1-Pi-1，j，k+1)/2，e8＝(Pi，j+1，k+1-Pi，j-1，k-1)/2，e9＝(Pi，j+1，k-1-Pi，j-1，k+1)/2，e10＝(Pi+1，j-1，k-1-Pi-1，j+1，k+1)/2，e11＝(Pi+1，j-1，k+1-Pi-1，j+1，k-1)/2，e12＝(Pi+1，j+1，k-1-Pi-1，j-1，k+1)/2，e13＝(Pi+1，j+1，k+1-Pi-1，j-1，k-1)/2。为使得fi，j，k(x，y，z)反映出医学图像的轮廓特性，用带约束的最小二乘法确定a7，a8，a9。目标函数G(a7，a8，a9)为：G(a7，a8，a9)＝w1(a7-e1)2+w2(a8-e2)2+w3(a9-e3)2+w4(a7+a8-e4)2+...+w13(a7+a8+a9-e13)2，其中，wi(i＝1，2，3，...，13)为权函数，通过联立式(2)和式(4)构建方程组求得：当中心点像素沿z方向接近线性变化时，则Δ3＝2a3≈0，w3与Δ3成反比；权函数wi(i＝1，2，3，...，13)的定义如下：其中，α和β为修正参数，取α和β为1。最小化目标函数即可求得a7，a8，a9：通过式(4)可知，式(2)推出如下形式：用带约束的最小二乘法确定a1，a2，...，a6，构建目标函数G(a1，a2，a3，。。。，a6)，其中，wp，q，l为权函数，定义规则同前；设置式(2)中的fi，j，k(x，y，z)沿z轴方向呈线性变化，那么Pi，j，k+1，Pi，j，k，Pi，j，k-1在同一条直线上，用式(7)中的w3代替，求出其他权重：进一步需要说明的是，步骤三还包括：当符合以下两种条件时，配置fi，j，k(u，v，w)为255×(fi，j，k(u，v，w)-fi，j，k(up，vp，wp))/fi，j，k本文档来自技高网...

【技术保护点】
1.一种基于三维空间体拟合的医学图像层间插值方法，其特征在于，方法包括：/n步骤一，配置三维空间体素拟合及配准的医学图像层间插值算法的两个约束条件；/n步骤二，构建三维空间体f

【技术特征摘要】
1.一种基于三维空间体拟合的医学图像层间插值方法，其特征在于，方法包括：
步骤一，配置三维空间体素拟合及配准的医学图像层间插值算法的两个约束条件；
步骤二，构建三维空间体fi，j，k(x，y，z)；
步骤三，修正空间体fi，j，k(u，v，w)；
步骤四，构造逼近空间体f(x，y，z)，对空间体进行重采样生成待插值医学图像。


2.根据权利要求1所述的方法，其特征在于，
步骤一还包括：
设相邻两个断层切片Sk，Sk+1之间的距离为1，待插值层间切片Sk+d(0＜d＜1)由m×n个体素Vi，j，k+di＝1，2，...，m，j＝1，2，...，n组成，待插值体素配置为原空间体F(x，y，z)上的采样值；
设每个体素Vi，j，k+d是从单位体积上采样得到，即



其中，w(x，y，z)是权函数；
F(x，y，z)定义在[0.5，m+0.5]×[0.5，n+0.5]×[0.5，t+0.5]Vi，j，k+d(i＝1，2，...，m；j＝1，2，...，n；k＝1，2，...，t)上，且对于每一个都满足式(1)；
同时，将权函数w(x，y，z)设置为1。


3.根据权利要求1或2所述的方法，其特征在于，
步骤二还包括：
在三维空间中，空间体fi，j，k(x，y，z)，i＝2，3，...，m-1，J＝2，3，...，n-1，k＝2，3，...，t-1满足以下条件：
如果式(1)中的F(x，y，z)是三元二次多项式，则fi，j，k(x，y，z)重建F(x，y，z)，即fi，j，k(x，y，z)＝F(x，y，z)；
令u＝x-i，v＝y-j，w＝z-k，则fi，j，k(x，y，z)在uvw空间的[-1.5，1.5]×[-1.5，1.5]×[-1.5，1.5]三维空间内写成：



若F(x，y，z)由式(2)定义，Pi，j，k由式(1)定义，则：



其中，
e1＝(Pi+1，j，k-Pi-1，j，k)/2，
e2＝(Pi，j+1，k-Pi，j-1，k)/2，
e3＝(Pi，j，k+1-Pi，j，k-1)/2，
e4＝(Pi+1，j+1，k-Pi-1，j-1，k)/2，
e5＝(Pi+1，j-1，k-Pi-1，j+1，k)/2，
e6＝(Pi+1，j，k+1-Pi-1，j，k-1)/2，
e7＝(Pi+1，j，k-1-Pi-1，j，k+1)/2，
e8＝(Pi，j+1，k+1-Pi，j-1，k-1)/2，
e9＝(Pi，j+1，k-1-Pi，j-1，k+1)/2，
e10＝(Pi+1，j-1，k-1-Pi-1，j+1，k+1)/2，
e11＝(Pi+1，j-1，k+1-Pi-1，j+1，k-1)/2，
e12＝(Pi+1，j+1，k-1-Pi-1，j-1，k+1)/2，



为使得fi，j，k(x，y，z)反映出医学图像的轮廓特性，用带约束的最小二乘法确定a7，a8，a9；目标函数G(a7，a8，a9)为：
G(a7，a8，a9)＝w1(a7-e1)2+w2(a8-e2)2+w3(a9-e3)2+w4(a7+a8-e4)2+...+w13(a7+a8+a9-e13)2，
其中，wi(i＝1，2，3，...，13)为权函数，通过联立式(2)和式(4)构建方程组求得：



当中心点像素沿z方向接近线性变化时，则Δ3＝2a3≈0，w3与Δ3成反比；
权函数wi(i＝1，2，3，...，13)的定义如下：



其中，α和β为修正参数，取α和β为1；最小化目标函数即可求得a7，a8，a9：



通过式(4)可知，式(2)推出如下形式：



用带约束的最小二乘法确定a1，a2，...，a6，构建目标函数G(a1，a2，a3，...，a6)，



其中，wp，q，l为权函数，定义规则同前；



设置式(2)中的fi，j，k(x，y，z)沿z轴方向呈线性变化，那么Pi，j，k+1，Pi，j，k，Pi，j，k-1在同一条直线上，用式(7)中的w3代替，求出其他权重：
w-1，0，0＝w1，0，0＝w1，w0，-1，0＝w0，1，0＝w2，
w0，-1，0＝w0，1，0＝w2，w1，i，0＝w-1，-1，0＝w4，
w1，-1，0...

【专利技术属性】
技术研发人员：刘慧，林毓秀，李钰，张彩明，
申请(专利权)人：山东财经大学，
类型：发明
国别省市：山东;37