CN111476888B - 一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质 - Google Patents
一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质 Download PDFInfo
- Publication number
- CN111476888B CN111476888B CN202010272244.0A CN202010272244A CN111476888B CN 111476888 B CN111476888 B CN 111476888B CN 202010272244 A CN202010272244 A CN 202010272244A CN 111476888 B CN111476888 B CN 111476888B
- Authority
- CN
- China
- Prior art keywords
- medical image
- space body
- image
- interpolation
- dimensional space
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 117
- 239000011229 interlayer Substances 0.000 title claims abstract description 65
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 25
- 238000005070 sampling Methods 0.000 claims abstract description 14
- 238000012952 Resampling Methods 0.000 claims abstract description 8
- 230000006870 function Effects 0.000 claims description 26
- 238000004590 computer program Methods 0.000 claims description 9
- 239000010410 layer Substances 0.000 description 18
- 230000003287 optical effect Effects 0.000 description 10
- 230000000007 visual effect Effects 0.000 description 10
- 238000004364 calculation method Methods 0.000 description 9
- 230000000694 effects Effects 0.000 description 9
- 230000008569 process Effects 0.000 description 9
- 238000012545 processing Methods 0.000 description 8
- 238000011156 evaluation Methods 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 6
- 238000010276 construction Methods 0.000 description 4
- 238000011158 quantitative evaluation Methods 0.000 description 4
- 238000002591 computed tomography Methods 0.000 description 3
- 238000013139 quantization Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 210000003484 anatomy Anatomy 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000000644 propagated effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 244000179970 Monarda didyma Species 0.000 description 1
- 235000010672 Monarda didyma Nutrition 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 238000013441 quality evaluation Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 208000011580 syndromic disease Diseases 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T15/00—3D [Three Dimensional] image rendering
- G06T15/005—General purpose rendering architectures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Analysis (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
Abstract
本发明提供一种基于三维空间体拟合的医学图像层间插值方法、设备及可读存储介质，配置三维空间体素拟合及配准的医学图像层间插值算法的两个约束条件；构建三维空间体；修正空间体；构造逼近空间体，对空间体进行重采样生成待插值医学图像。对空间体进行重采样生成质量较好的待插值层间图像.通过采样公式对原图像序列反向采样局部地构造三元二次多项式空间体,将二次多项式空间体加权平均生成单位空间体,再将所有的单位空间体拼合成拟合三维空间体,最后对空间体进行重采样生成待插值图像.方法生成的层间图像的质量不论在视觉评价还是客观数值上均有所提升,且有效地提高了插值效率。
Description
背景技术
断层医学图像层间插值技术是医学图像超分辨重建和三维重建的关键技术之一。即,在相邻两层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(uc，vc，wc)，使
采样得到的体素Vi，j，k满足条件0≤Pi，j，k≤255，
fi，j，k(u，v，w)在三维空间，Ωi，j，k＝[i-1，j-1，k-1]×[i+1，j+1，k+1]上，由fi，j，k(u，v，w)加权平均生成的空间体f(u，v，w)也满足0≤f(u，v，w)≤255；
①fi，j，k(u，v，w)在(xp，yp，zp)点处达到最小值fi，j，k(up，vp，wp)＜0；
②fi，j，k(u，v，w)在(xc，yc，zc)点处达到最大值fi，j，k(uc，vc，wc)＞255。
进一步需要说明的是，步骤四还包括：
在每个区域
[i，i+1]×[j，j+1]×[k，k+1]，i＝1，2，...，m-1，j＝1，2，...，n-1，
j＝1，2，...，n-1，k＝1，2，...，t-1
上，构造一个单位空间体Bi，j，k(x，y，z)，将所有的Bi，j，k(x，y，z)拼接在一起构成f(x，y，z)；对于构造单位空间体Bi，j，k(x，y，z)，
当i＝2，3，...，m-1，j＝2，3，...，n-1，k＝2，...，t-1时，区域[i，i+1]×[j，j+1]×[k，k+1]上的空间体Bi，j，k(x，y，z)
由fi，j，k(x，y，z)，fi+1，j，k(x，y，z)，fi，j+1，k(x，y，z)，fi，j，k+1(x，y，z)，fi+1，j+1，k(x，y，z)，fi+1，j，k+1(x，y，z)，fi，j+1，k+1(x，y，z)fi+1，j+1，k+1(x，y，z)加权平均生成，即
Bi，j，k(x，y，z)＝wi，j，k(x，y，z)fi，j，k(x，y，z)+…+wi+1，j+1，k+1(x，y，z)fi+1，j+1，k+1(x，y，z)+wi，j+1，k+1(x，y，z)fi+1，j+1，k+1(x，y，z)，
其中
是权值，u＝x-i，v＝y-j，w＝z-k；
对于在边界区域构造空间体
B1，1，1(x，y，z)由f2，2，2(x，y，z)定义，当i＝1，j＝1，k＝2，3，...，t-2时，B1，1，k(x，y，z)由两个空间体定义，即
B1，1，k(x，y，z)＝f2，2，k(x，y，z)(1-w)+f2，2，k+！(x，y，z)w。
当i＝1，j＝2，3，...，n-2，k＝2，3，...，t-2时，B1，j，k(x，y，z)由四个空间体定义，即
B1，j，k(x，y，z)＝f2，j，k(x，y，z)(1-v)(1-w)+f2，j+1，k(x，y，z)v(1-w)+f2，j，k+1(x，y，z)(1-v)w+f2，j+1，k+1(x，y，z)vw.。
本发明还提供一种实现基于三维空间体拟合的医学图像层间插值方法的设备，其特征在于，包括：
存储器，用于存储计算机程序及基于三维空间体拟合的医学图像层间插值方法；
处理器，用于执行所述计算机程序及基于三维空间体拟合的医学图像层间插值方法，以实现基于三维空间体拟合的医学图像层间插值方法的步骤。
本发明还提供一种具有基于三维空间体拟合的医学图像层间插值方法的可读存储介质，可读存储介质上存储有计算机程序，所述计算机程序被处理器执行以实现基于三维空间体拟合的医学图像层间插值方法的步骤。
从以上技术方案可以看出，本发明具有以下优点：
本发明在三维空间中利用连续三幅CT图像反向重建出逼近原场景的拟合空间体，并对空间体进行重采样生成质量较好的待插值层间图像。通过采样公式对原图像序列反向采样局部地构造三元二次多项式空间体，将二次多项式空间体加权平均生成单位空间体，再将所有的单位空间体拼合成拟合三维空间体，最后对空间体进行重采样生成待插值图像。将本发明方法与线性插值和基于形状的插值等插值方法在3组CT图像序列上进行比较，相较于其他现有算法,方法生成的层间图像的质量不论在视觉评价还是客观数值上均有所提升,且有效地提高了插值效率。本发明的方法减少计算时间、提高计算效率的同时,插值出的层间图像在观测质量和量化精度上均获得了不错的效果.
附图说明
为了更清楚地说明本发明的技术方案，下面将对描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为断层图像层间插值示意图；
图2为基于三维空间体拟合的医学图像层间插值方法流程图；
图3为幅图像构造三维空间体图；
图4为幅图像构造三维空间体图；
图5为算法框架图；
图6为体素沿13个方向变化图；
图7为选取的3组测试图像序列图；
图8为四种不同插值方法对数据集Ⅰ的插值结果图；
图9为四种不同插值方法对数据集Ⅱ的插值结果图；
图10为四种不同插值方法对数据集Ⅲ的插值结果图；
图11为原始图像与不同插值结果局部对比图。
具体实施方式
本领域普通技术人员可以意识到，结合本发明中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。
三维计算机断层扫描(Three-Dimensional Computed Tomography,3D-CT)医学成像方法通常将获取的断层医学影像数据表示为一组切片序列。CT图像是层面图像,常用的是横断面。层间距是指相邻两层切片中心之间的距离。受成像设备、辐射剂量等因素的影响,国内用于临床的CT层间距一般为2。5mm～5mm。CT切片间的层间距大于切片内像素间距离,层间分辨率低于层内分辨率,导致CT切片序列中的体素呈各向异性。而在可视化、医学图像处理中,要求体素的维数是各向同性或者是接近各向同性的,例如,多幅图像超分辨重建、三维重建、图像分割等。因此,需要对断层图像进行插值,增加切片数量,以便CT图像各层之间过渡更加平滑，减弱各向异性。
本发明通过一种新颖的层间插值策略提高断层医学图像的层间分辨率,以获得高分辨率的体数据。经过插值后的切片序列层间分辨率等于或接近层内分辨率,消除或减弱体素的各向异性,能够有效改善后续图像超分辨重建和图像三维重建结果,辅助临床分析和诊断,对帮助医生提高临床诊断准确率具有重要的现实意义。
本发明将平面图像放缩中层内像素插值的理论用于三维空间内层间体素插值。通常地,二维图像是对场景表面的采样,通过求解原图像采样的逆过程反向重建出原场景,再重采样获得高质量的放缩图像。本发明将该理论扩展到三维空间,以获得精确的非刚体医学CT图像的估计配准信息；
本发明提出一种基于三维空间体素拟合及配准的医学图像层间插值方法。在三维空间中,利用三幅相邻层图像构造拟合空间体,并对构造的空间体进行重采样生成待插值图像,以在取得的好的视觉效果的同时大大缩短时间复杂度。
为解决技术问题存在插值精度低或插值效率低等缺点,本发明提出一种基于三维空间体素拟合及配准的医学图像层间插值方法,在减少计算时间、提高计算效率的同时,插值出的层间图像在观测质量和量化精度上均获得了不错的效果。
具体步骤包括如图2和图3所示，图像放缩过程中,一幅图像是对一个场景表面的采样,通过求出原图像采样的逆过程反向重建出原场景,就可获得高质量的放缩图像。本发明假设三维空间中已知的三幅图像可以构成一个三维空间体,图3和图4所示，空间体可由多个分块空间体表示。通过反向采样过程构造已知图像序列的原三维空间体,再对该空间区域重采样生成新的体素点,以提高插值图像的精度。具体地,由于原图像数据的有限性,很难根据有限的采样数据重建出原空间体,因此,本发明方法构造原空间体的逼近空间体。此外,医学图像具有非刚体特征,轮廓边缘形变是非线性的,于是,本发明采用带约束的空间体构造方法,能够构造出视觉效果良好的插值图像。
在发明中，步骤一，配置三维空间体素拟合及配准的医学图像层间插值算法的两个约束条件；
具体的，如图5所示，基于三维空间体素拟合及配准的医学图像层间插值算法依赖于两个约束条件:
(1)待插值CT图像和原始CT图像之间的像素是连续变化的,它们含有相似的解剖结构。
(2)待插值CT图像和原始两幅CT图像的相似度与它距离这两幅图像的距离成反比关系。
假设相邻两个断层切片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)是权函数,在不同应用场景中可以设置不同的权函数来获得更精确的体素值。
因为体素值在一般情况下是整数形式,所以,式(1)近似成立且误差
小于0。5。若已知F(x,y,z),则由式(1)可以很容易地求出相应的Vi,j,k+d。因此,如何构造合理的F(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。
步骤二，构建三维空间体fi,j,k(x,y,z)；
由于F(x,y,z)可以是任意多项式函数,仅依靠原始图像序列很难准确构造出F(x,y,z),只能构造对原三维空间体的逼近空间体fi,j,k(x,y,z)。根据数值逼近理论,任何连续的函数
都可以在连续三维空间内的一点(xa,ya,za)处展开成多项式泰勒级数,该泰勒级数是在点(xa,ya,za)的空间邻域内对
的逼近。假设F(x,y,z)在每个体素点的空间邻域内可由泰勒级数的前三项逼近,即该空间体具有二次多项式逼近精度。
在三维空间中,空间体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]三维空间内可写成:
其中,a1,a2,…,a10都是未知量。
先确定a7,a8,a9,再求剩余未知量。在这里为方便说明,以Pi,j,k为例。
推论1.若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.
证明
联立式(1)和式(2)求方程组得
插值图像轮廓处的质量影响图像的直观视觉效果,因此,在生成待插值图像时,要使得三维空间体fi,j,k(x,y,z)的各个曲面都尽可能反映出图像的边缘特性。在三维数据空间中,每一个体素与周围相邻体素(上下、左右、前后、对角)具有某种相关性。如图6所示,在中心体素与其邻域体素形成的13个方向中,如果图像像素沿某个方向变化比较缓慢,则该方向为图像边缘的可能性比较大,如果图像像素沿某个方向呈线性变化,则希望fi,j,k(x,y,z)沿该方向为线性函数。式(3)中有13个等式,三个未知量a7,a8,a9,为使得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)构建方程组求得:
如果式(2)中的fi,j,k(x,y,z)沿z轴方向呈线性变化,则e3应当对a9起决定作用,也就是要使w3对应一个相当大的值。当中心点像素沿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在同一条直线上,Pi,j,k-1以及Pi,j,k+1应当对fi,j,k(x,y,z)沿z轴方向变化起决定性作用,也就是要使w0,0,-1和w0,0,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,1,0＝w-1,-1,0＝w4,
w1,-1,0＝w-1,1,0＝w5,w1,0,1＝w-1,0,-1＝w6,
w1,0,-1＝w-1,0,1＝w7,w0,1,1＝w0,-1,-1＝w8,
w0,1,-1＝w0,-1,1＝w9,w1,-1,-1＝w-1,1,1＝w10,
w1,-1,1＝w-1,1,-1＝w11,w1,1,-1＝w-1,-1,1＝w12,
w1,1,1＝w-1,-1,-1＝w13.。
本发明涉及的步骤三，修正空间体fi,j,k(u,v,w)；
具体的，为确保采样得到的体素Vi,j,k满足条件0≤Pi,j,k≤255,fi,j,k(u,v,w)在三维空间Ωi,j,k,Ωi,j,k＝[i-1,j-1,k-1]×[i+1,j+1,k+1]上,需要满足0≤fi,j,k(u,v,w)≤255,使得由fi,j,k(u,v,w)加权平均生成的空间体f(u,v,w)也满足0≤f(u,v,w)≤255。
因此,当符合以下两种情况时,
修改fi,j,k(u,v,w)为255×(fi,j,k(u,v,w)-fi,j,k(up,vp,wp))/fi,j,k(uc,vc,wc)。
①fi,j,k(u,v,w)在(xp,yp,zp)点处达到最小值fi,j,k(up,vp,wp)<0；
②fi,j,k(u,v,w)在(xc,yc,zc)点处达到最大值fi,j,k(uc,vc,wc)＞255。
本发明涉及的步骤四，构造逼近空间体f(x,y,z)，对空间体进行重采样生成待插值医学图像。
本发明用空间体fi,j,k(x,y,z),i＝2,3,…,m-1,j＝2,3,…,n-1,k＝2,3,...,t-1确定逼近空间体f(x,y,z)。具体步骤是:
在每个区域
[i,i+1]×[j,j+1]×[k,k+1],i＝1,2,...,m-1,j＝1,2,...,n-1,j＝1,2,...,n-1,k＝1,2,...,t-1上,构造一个单位空间体Bi,j,k(x,y,z),将所有的Bi,j,k(x,y,z)拼接在一起构成f(x,y,z)。对于构造单位空间体Bi,j,k(x,y,z),本发明分两种情况讨论。
①当i＝2,3,…,m-1,j＝2,3,...,n-1,k＝2,…,t-1时,区域[i,i+1]×[j,j+1]×[k,k+1]上的空间体Bi,j,k(x,y,z)由fi,j,k(x,y,z),fi+1,j,k(x,y,z),fi,j+1,k(x,y,z),fi,j,k+1(x,y,z),fi+1,j+1,k(x,y,z),fi+1,j,k+1(x,y,z),fi,j+1,k+1(x,y,z)fi+1,j+1,k+1(x,y,z)加权平均生成,即
Bi,j,k(x,y,z)＝wi,j,k(x,y,z)fi,j,k(x,y,z)+…+wi+1,j+1,k+1(x,y,z)fi+1,j+1,k+1(x,y,z)+wi,j+1,k+1(x,y,z)fi+1,j+1,k+1(x,y,z),
其中
是权值,u＝x-i,v＝y-j,w＝z-k。
②对于在边界区域构造空间体
如图4所示,B1,1,1(x,y,z)由f2,2,2(x,y,z)定义,当i＝1,j＝1,k＝2,3,...,t-2时,B1,1,k(x,y,z)由两个空间体定义,即
B1,1,k(x,y,z)＝f2,2,k(x,y,z)(1-w)+f2,2,k+！(x,y,z)w.
当i＝1,j＝2,3,...,n-2,k＝2,3,...,t-2时,B1,j,k(x,y,z)由四个空间体定义,即
B1,j,k(x,y,z)＝f2,j,k(x,y,z)(1-v)(1-w)+f2,j+1,k(x,y,z)v(1-w)+f2,j,k+1(x,y,z)(1-v)w+f2,j+1,k+1(x,y,z)vw.
由于对称性,其他情况同样遵循上述规律。
为验证本发明方法的有效性,分别用三组符合DICOM标准的不同层间距CT数据进行仿真实验,包括人体胸部(冠状面、横断面)和头部断层图像,这三种CT图像类别在医学影像处理研究领域中十分常见。每个数据类别都由的独立的图像序列组成,CT图像序列均由山东省千佛山医院提供。图5提供了示例图像,表1列出了相关的成像参数。
插值得到的层间图像利用误差图像进行视觉评价和量化评价缺乏真实图像依据。因此,在实验过程中,分别从Ⅰ、Ⅱ、Ⅲ组图像序列中随机选取连续4幅断层图像,先用前1幅和后2幅断层图像来插值中间切片,第2幅图像模拟中间断层图像用来评定插值效果；然后用新插值出的切片与中间断层图像进行对比。本发明算法由C++编写实现,运行环境为Windows 10,处理器为Intel(R)Core(TM)i7-4790 CPU@3。60GHz,内存为16GB。
表1实验数据信息
在数值误差定量评估中,本发明选择均方根误差(Root Mean Square Error,RMSE),DSSIM(Structural Dissimilarity)2种评价指标进行定量评估,这2种方法广泛应用于图像质量评价中。
RMSE是真实图像IGT图像与估计插值图像Iip的均方根误差,
其中,M×N为图像大小,均方根误差可以用来衡量插值图像与模拟真实图像之间的像素误差。RMSE值越小,说明该实验插值生成的新图像具有更高的精确度,更接近真实CT图像。RMSE针对异常值比较敏感,即,有一个预测值与模拟值相差较大，那么RMSE值就会较大。
DSSIM是常见的衡量图像间相似程度的指标,可以更好地估计局部差异,它是基于结构相似性(structural similarity,SSIM)得出的一个距离度量。
其中,μx和μy分别代表模拟真实图像和插值图像的灰度均值,σx和σy为标准差,σxy是x和y的协方差。C1＝(k1L)2,C2＝(k2L)2,且k1＝0.01,k2＝0.03,L＝255。DSSIM值越小,表示图像失真程度越小,更接近真实图像。
为了综合评价所提出的插值方法,基于三组不同部位的CT图像数据集进行了相关仿真实验测试。主要从客观数据、视觉效果以及算法运行时间3个方面来评价本发明算法的插值效果。
本发明分别使用线性插值方法、基于形状的插值方法、基于光流插值方法和本发明算法进行插值处理和比较。之所以选择上述4种方法进行对比,是因为线性插值是基于灰度插值的常用算法,基于光流的插值方法在非刚性配准插值中表现良好,它们均是被诸多发明献引用且比较过的算法。
首先选取3组原始医学图像序列中的指定层,并将该层的原图像模拟无失真层间图像,再利用上述4种插值算法生成该指定层,与插值图像进行比较。由于进行仿真实验,本发明所用的Ⅲ组实验数据实际层间距在原层间距基础上均扩大一倍。因此用中间隔一幅原始断层图像的相邻图像插值出的中间结果与原始图像有一定的差别是符合客观实际的。
表2给出了4种不同插值算法的量化数据,每一行数据从上到下分别是RMSE和DSSIM。结合表2数据可以看出,相对于比较算法,本发明提出的算法拥有最低平均RMSE和DSSIM值,在客观数值评价中具有较大的优势。本发明算法的RMSE比线性插值平均降低了约37.6％,比基于形状的插值平均降低了22.2％,比基于光流的插值平均降低11.2％。通过对三组不同数据集的量化结果进行分析,本发明方法在第1组和第2组图像序列上表现好于第3组图像序列,可以看出本发明方法能够较好地插值层间距较小的图像,而对于较大的层间距则不太明显。实验表明本发明方法在层间插值的表现要优于线性插值,略好于或相当于基于形状的插值和基于光流的插值,说明该实验插值生成的新图像具有较高的精确度,接近真实CT图像。
在视觉评价方面,图7～图11分别展示了上述4种方法在三组不同数据集上的定性评估结果,红色框中为图像主要变化区域。通过比较原始切片与算法插值结果发现,本发明方法插值图像在主观视觉上与原始切片非常接近,边缘模糊现象显著改善。与基于形状的插值方法和基于光流配准的插值方法相比,提出的方法可达到优于或相似的插值效果。
为了更加清晰地比较4种插值结果的细微区别,在图11显示出了图8～图10红色矩形框内的局部放大图。通过对图像细节的放大,能够从视觉效果上更直观地评估各种插值方法生成图像的质量。从图中可以看出,线性插值因为只是将待插值切片的体素中心作为插值点,并进行简单地线性计算,无法准确反映非刚性形变过程,有明显的边缘模糊,层间插值结果较差。而本发明方法是对相邻切片采样构造三维空间体,再重采样生成待插值体素点,因此边缘模糊现象较线性插值方法明显改善,与原切片差异明显减小。基于形状的插值方法伪影明显改善,但在形变细节较多的情况如图9所示插值图像质量不高。基于光流配准的插值方法虽然能够反映非刚性形变过程，有效消除边缘伪影,但是会出现细小孔洞现象,影像插值图像质量。综合来说本发明方法可以明显改善边缘重影且对小形变的处理能力较强。
表2不同插值算法RMSE和DSSIM值对比
在实际的应用过程中,不仅要求获得较高的图像质量,算法的运算速度也非常关键。表3列出了分别使用所评估的4种插值方法在三组不同数据集上插值单幅图像所需的时间。可以看出,虽然本发明方法的性能略好于或相当于基于形状和基于光流的算法，但其在计算效率上要好很多,对单个图像进行插值所需的处理时间减少了一个数量级。与线性插值方法时间效率接近,但其插值效果在定性和定量上均不如本发明方法。本发明方法不仅在视觉评价和客观数值上取得了不错的插值效果,且时间复杂度较低。
表3不同插值算法运行时间对比
本发明提出了一种基于三维空间体素拟合及配准的医学图像层间插值方法,并使用本发明方法和三种常用的插值技术:线性插值、基于形状插值和基于TV-L1光流配准的插值进行了定量评估和视觉比较。每种方法的相对性能都是在同一组CT扫描序列上进行评估的。本发明方法在处理层间距较小图像方面,其与模拟真实图像的体素计算显示出比其他三种插值技术更低的误差,同时时间复杂度较低。
本发明还提供一种实现基于三维空间体拟合的医学图像层间插值方法的设备，包括：
存储器，用于存储计算机程序及基于三维空间体拟合的医学图像层间插值方法；
处理器，用于执行所述计算机程序及基于三维空间体拟合的医学图像层间插值方法，以实现基于三维空间体拟合的医学图像层间插值方法的步骤。
本发明还提供一种具有基于三维空间体拟合的医学图像层间插值方法的可读存储介质，可读存储介质上存储有计算机程序，所述计算机程序被处理器执行以实现基于三维空间体拟合的医学图像层间插值方法的步骤。
实现基于三维空间体拟合的医学图像层间插值方法的设备是结合本发明中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。
所属技术领域的技术人员能够理解，实现基于三维空间体拟合的医学图像层间插值方法的设备的各个方面可以实现为系统、方法或程序产品。因此，本公开的各个方面可以具体实现为以下形式，即：完全的硬件实施方式、完全的软件实施方式(包括固件、微代码等)，或硬件和软件方面结合的实施方式，这里可以统称为“电路”、“模块”或“系统”。
在实现基于三维空间体拟合的医学图像层间插值方法的设备的存储介质中，可读信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了可读程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。可读信号介质还可以是可读存储介质以外的任何可读介质，该可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本发明中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本发明所示的这些实施例，而是要符合与本发明所公开的原理和新颖特点相一致的最宽的范围。
Claims (3)
1.一种基于三维空间体拟合的医学图像层间插值方法，其特征在于，方法包括：
步骤一，配置三维空间体素拟合及配准的医学图像层间插值算法的两个约束条件；
设相邻两个断层切片, 之间的距离为1, 待插值层间切片由个体素 , 组成, 待插值体素配置为原空间体上的采样值；
设每个体素是从单位体积上采样得到, 即
； (1)
其中, 是权函数；
定义在
上, 且对于每一个都满足式(1)；
同时, 将权函数设置为1；
步骤二，构建三维空间体；
在三维空间中,空间体,
满足以下条件：
如果式(1)中的是三元二次多项式, 则重建, 即；
令, , , 则在空间的三维空间内写成:
(2)
若由式(2)定义, 由式(1)定义, 则：
(3)
其中,
为使得反映出医学图像的轮廓特性, 用带约束的最小二乘法确定；目标函数为：
其中, 为权函数, 通过联立式(2)和式(4)构建方程组求得:
(6)
当中心点像素沿方向接近线性变化时, 则, 与成反比；
权函数的定义如下：
, (7)
其中, 和为修正参数，取和为1; 最小化目标函数即可求得:
通过式(4)可知, ；式(2)推出如下形式：
用带约束的最小二乘法确定，构建目标函数,
其中, 为权函数, 定义规则同前;
;
设置式(2)中的沿z轴方向呈线性变化, 那么, , 在同一条直线上，用式(7)中的代替,求出其他权重：
步骤三，修正空间体；
当符合以下两种条件时, 配置为，使
采样得到的体素满足条件,
在三维空间,上,由加权平均生成的空间体也满足；
①在点处达到最小值；
②在点处达到最大值；
步骤四，构造逼近空间体，对空间体进行重采样生成待插值医学图像；
在每个区域
k＝1，2，...，t-1
上, 构造一个单位空间体, 将所有的拼接在一起构成;对于构造单位空间体,
当时, 区域上的空间体
由
加权平均生成, 即
其中
是权值, , , ；
对于在边界区域构造空间体
由定义,
当时, 由两个空间体定义, 即
当时, 由四个空间体定义, 即
。
2.一种实现基于三维空间体拟合的医学图像层间插值方法的设备，其特征在于，包括：
存储器，用于存储计算机程序及基于三维空间体拟合的医学图像层间插值方法；
处理器，用于执行所述计算机程序及基于三维空间体拟合的医学图像层间插值方法，以实现如权利要求1所述基于三维空间体拟合的医学图像层间插值方法的步骤。
3.一种具有基于三维空间体拟合的医学图像层间插值方法的可读存储介质，其特征在于，所述可读存储介质上存储有计算机程序，所述计算机程序被处理器执行以实现如权利要求1所述基于三维空间体拟合的医学图像层间插值方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010272244.0A CN111476888B (zh) | 2020-04-09 | 2020-04-09 | 一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010272244.0A CN111476888B (zh) | 2020-04-09 | 2020-04-09 | 一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111476888A CN111476888A (zh) | 2020-07-31 |
| CN111476888B true CN111476888B (zh) | 2023-05-12 |
Family
ID=71751329
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010272244.0A Active CN111476888B (zh) | 2020-04-09 | 2020-04-09 | 一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111476888B (zh) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102835974A (zh) * | 2012-08-23 | 2012-12-26 | 华南理工大学 | 基于并行计算机的医学超声三维成像方法 |
| CN103025229A (zh) * | 2010-04-29 | 2013-04-03 | 麻省理工学院 | 适用于光学相干断层扫描技术的移动修正和图像增强的方法和装置 |
| CN110570515A (zh) * | 2019-09-03 | 2019-12-13 | 天津工业大学 | 一种利用ct图像进行人体骨骼三维建模的方法 |
-
2020
- 2020-04-09 CN CN202010272244.0A patent/CN111476888B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103025229A (zh) * | 2010-04-29 | 2013-04-03 | 麻省理工学院 | 适用于光学相干断层扫描技术的移动修正和图像增强的方法和装置 |
| CN102835974A (zh) * | 2012-08-23 | 2012-12-26 | 华南理工大学 | 基于并行计算机的医学超声三维成像方法 |
| CN110570515A (zh) * | 2019-09-03 | 2019-12-13 | 天津工业大学 | 一种利用ct图像进行人体骨骼三维建模的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111476888A (zh) | 2020-07-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhou et al. | Handbook of medical image computing and computer assisted intervention | |
| Zeng et al. | Simultaneous single-and multi-contrast super-resolution for brain MRI images based on a convolutional neural network | |
| Kainz et al. | Fast volume reconstruction from motion corrupted stacks of 2D slices | |
| CN107610194B (zh) | 基于多尺度融合cnn的磁共振图像超分辨率重建方法 | |
| Trinh et al. | Novel example-based method for super-resolution and denoising of medical images | |
| Mair et al. | Estimation of images and nonrigid deformations in gated emission CT | |
| Liu et al. | Motion artifacts reduction in brain MRI by means of a deep residual network with densely connected multi-resolution blocks (DRN-DCMB) | |
| Albu et al. | A morphology-based approach for interslice interpolation of anatomical slices from volumetric images | |
| CN110070612B (zh) | 一种基于生成对抗网络的ct图像层间插值方法 | |
| CN110211193B (zh) | 三维ct层间图像插值修复与超分辨处理方法及装置 | |
| van ˈt Klooster et al. | Automated registration of multispectral MR vessel wall images of the carotid artery | |
| CN112634265B (zh) | 基于dnn的胰腺全自动分割模型的构建、分割方法及系统 | |
| Sander et al. | Autoencoding low-resolution MRI for semantically smooth interpolation of anisotropic MRI | |
| Chan et al. | An attention-based deep convolutional neural network for ultra-sparse-view CT reconstruction | |
| Safari et al. | MRI motion artifact reduction using a conditional diffusion probabilistic model (MAR‐CDPM) | |
| Prince et al. | Image synthesis and superresolution in medical imaging | |
| Terada et al. | Clinical evaluation of super-resolution for brain MRI images based on generative adversarial networks | |
| CN109741439B (zh) | 一种二维mri胎儿图像的三维重建方法 | |
| CN111476888B (zh) | 一种基于三维空间体拟合的医学图像层间插值方法，设备及可读存储介质 | |
| JP2022088320A (ja) | 体動補正装置、方法及びプログラム | |
| Okamoto et al. | Patch-based artifact reduction for three-dimensional volume projection data of sparse-view micro-computed tomography | |
| CN111932443A (zh) | 多尺度表达结合造影剂提升超声与磁共振配准精度的方法 | |
| Wang et al. | 3D image interpolation based on directional coherence | |
| Xu | A Robust and Efficient Framework for Slice-to-Volume Reconstruction: Application to Fetal MRI | |
| Liu et al. | Slice Interpolation for Medical Image based on Spatial Geometry Polynomial Fitting |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |