一种亚像素精度的快速三维超声心动图散斑跟踪方法

1.本发明涉及一种三维(3d)超声心动图散斑跟踪方法，尤其涉及一种亚像素精度的快速3d超声心动图散斑跟踪方法，属于医学超声图像处理技术领域，旨在实现心脏功能精确评价和定量诊疗的3d ste成像技术(3d speckle trackingechocardiography)。


背景技术：

2.心力衰竭是冠心病、高血压、先天性心脏病、瓣膜性心脏病、心肌病等心血管疾病发展的最后阶段，是造成心血管疾病患者死亡的主要原因和临床治疗的最后战场，被认为是21世纪心血管界面临的最大挑战。对易患心力衰竭的高危人群，早期进行心脏功能的定量评价，及早识别心脏收缩与舒张功能异常并进行有效的干预，有利于预防和延缓心力衰竭的发生和发展，从而有助于延长心血管疾病患者的寿命和改善预后。散斑跟踪三维超声心动图(3d speckle trackingechocardiography，3d ste)正是这样一种新兴的、旨在实现心脏功能精确评价和定量诊疗的超声成像技术，它通过对3d超声心动图散斑(speckle)或散斑模式(speckle pattern)的跟踪和运动分析，能精确地实现心肌组织的速度向量成像、应变和应变率成像，以及左心室扭转角位移、扭转力矩、扭转角速度等运动参数的测量与显像，从而为临床医生客观评估心肌的整体和局部功能、心室的收缩与舒张功能提供定量分析依据，也为临床医生定量研究心肌组织的生物弹性力学特性提供技术手段。
3.超声心动图散斑跟踪方法中要跟踪分析的散斑或散斑模式，在图像处理实践中则是位于心肌组织内的若干感兴趣区域(roi)，亦即传统块匹配处理的块或核(block/kernel)。而见诸文献报道的超声心动图散斑跟踪方法主要有两类：块匹配的散斑跟踪方法和光流场的散斑跟踪方法。光流场的散斑跟踪方法(optical flow based speckle tracking)能实现较高精度的心肌组织运动跟踪，但对图像质量要求苛刻，临床实用性差。块匹配的散斑跟踪方法(block matching based speckle tracking)临床实用性较好，是3d ste成像技术领域主流的方法，但传统的块匹配方法难以实现高精度的心肌组织运动跟踪。另外，以往方法的计算处理效率以及自适应特性和鲁棒性等临床应用性能普遍欠佳。


技术实现要素：

4.针对以往方法或者跟踪精度低，或者临床实用性差，普遍难以满足心脏功能精确评价和定量诊疗临床应用需求的弊端，本着性能良好、临床实用的技术原则，本发明提出了一种具有良好自适应特性和鲁棒性、1/16亚像素精度的快速3d超声心动图散斑跟踪方法。
5.术语解释：
6.1、三次样条插值，cubic spline interpolation，其中，样条指穿过一系列形值点的光滑曲线，样条插值则是计算获得此类曲线的数学方法。
7.2、全搜索算法，full search，图像处理中实现块匹配运算的一种策略，通过评估搜索区域内的全部候选点获得匹配块。
8.3、sad，sum of absolute differences，绝对误差和。
9.4、成像系统3d图像空间，3d医学成像领域的技术术语，在本发明中指3d超声心动图序列图像所在的笛卡尔空间，如图3所示，三维是指笛卡尔空间的x、y、z坐标三维。
10.本发明的技术方案如下：
11.一种3d超声心动图散斑跟踪的方法，由一台配置了心脏功能精确评价和定量诊疗功能的计算机工作站来实现。设3d超声心动图序列图像由n帧图像组成，分别记作i
k
，k＝1,2,...,n，其中i
k
为3d超声心动图序列图像的第k帧图像，i1为初始帧，i
n
为末帧，i1,i3,...,i
n
为奇数帧及末帧构成的序列图像；将第t帧图像进行块匹配所使用的参考块记作r
t
，t＝1,3,5,...,n，第t帧图像进行块匹配所使用的、经像素值校准后的参考块记作将在第t帧图像上获得的最佳匹配块记作m
t
，将在第t帧图像上获得的跟踪点在成像系统3d图像空间中的位置记作r
t
(x
t
,y
t
,z
t
)，其中(x
t
,y
t
,z
t
)表示跟踪点在3d图像空间中的坐标，该方法的步骤如下：
12.s1)由临床医师在初始帧i1上手动选取跟踪点；
13.在3d超声心动图序列图像的初始帧i1上，由临床医师手动选取跟踪点，该跟踪点在成像系统3d图像空间中的位置为r1(x1,y1,z1)；
14.s2)设定初始参数；
15.设定以各帧的跟踪点为中心的心肌组织感兴趣区域的大小n
x
×
n
y
×
n
z
，其中n
x
、n
y
、n
z
分别在成像系统3d图像空间中为沿x、y、z轴方向的大小；设定搜索区域大小；设定参考块像素值校准的因子m；令当前帧序号j＝1；设定初始参考块大小，以r
j
(x
j
,y
j
,z
j
)点为中心形成初始参考块r
j
；令r
j
是指在第j帧图像进行块匹配所使用的参考块，是指第j帧图像进行块匹配所使用的、经像素值校准后的参考块，m
j
是指在第j帧图像上获得的最佳匹配块；
16.s3)参考块像素值校准；
17.j增加2；令r
j
＝m
j
‑2；
18.对参考块r
j
的像素值进行校准，得到第j帧图像进行块匹配所使用的、经像素值校准后的参考块
19.s4)进行1/16亚像素精度的块匹配计算，获得当前帧跟踪点在成像系统3d图像空间中的位置的坐标r
j
(x
j
,y
j
,z
j
)；
20.s5)若i
j
为末帧，则执行s6，否则执行s3；
21.s6)计算偶数帧跟踪点在成像系统3d图像空间中的位置
22.对奇数帧和末帧跟踪点在成像系统3d图像空间中的坐标，分别取其x轴、y轴、z轴坐标，通过三次样条插值计算偶数帧跟踪点在成像系统3d图像空间中的x轴、y轴、z轴坐标，获得偶数帧跟踪点在成像系统3d图像空间中的位置，进而得到所有帧跟踪点在成像系统3d图像空间中的位置；
23.s7)形成并记录所有帧的心肌组织感兴趣区域
24.以所有帧的跟踪点为中心，形成并记录各帧大小为n
x
×
n
y
×
n
z
的感兴趣区域，实现3d超声心动图的散斑跟踪。
25.根据本发明优选的，对参考块r
j
的像素值进行校准，如式(ⅰ)所示：
[0026][0027]
式(ⅰ)中，0<m,p,q<1，m+p+q＝1，p＝0.9(1
‑
m)/(n
‑
3)2·
(j
‑
3)2，q＝1
‑
m
‑
p。
[0028]
根据本发明优选的，以为参考块，执行1/16亚像素精度的块匹配计算，实现过程如下：
[0029]
s41)在i
j
上对以r
j
‑2(x
j
‑2,y
j
‑2,z
j
‑2)点为中心的搜索区域进行第一次三次样条插值，使得2个相邻节点之间3个等间隔的内点被插值，获得1/4亚像素精度的搜索区域s11；
[0030]
s42)采用全搜索算法及sad(sum of absolute differences)匹配准则，并使用sad迭代计算提前终止策略，在搜索区域s11中进行一次匹配，记最小sad值对应的点为一次1/4亚像素精度跟踪点r11；
[0031]
s43)以一次1/4亚像素精度跟踪点r11为中心点，对周围8个1/4体素块进行第二次三次样条插值，继而2个相邻节点之间3个等间隔的内点被插值，获得1/16亚像素精度的搜索区域s12。
[0032]
s44)通过计算1/8亚像素精度跟踪点对搜索区域s12进行压缩处理：
[0033]
计算8个1/4体素块中心点b1、b2、b3、b4、b5、b6、b7、b8的sad值，取其中最小的sad值对应的点(以b5为例)为中心划定1/4体素块为二次搜索区域s13；
[0034]
s44)采用全搜索算法及sad匹配准则，并使用sad迭代计算提前终止策略，进行二次匹配，最小sad值对应的点即为当前帧(第j帧)1/16亚像素精度的跟踪点r
j
(x
j
,y
j
,z
j
)，以该点为中心的块为当前帧(第j帧)上的最佳匹配块m
j
。
[0035]
根据本发明优选的，步骤s42)及步骤s44)中，具体实现过程如下：
[0036]
设ca
max
为全搜索候选块总数，l
max
为全搜索候选块z轴方向层数，sad
min
为当前最小sad值，r代表当前计算第r个候选块，sad
rs
代表第r个候选块的前s层误差和，sad
s
代表第s层误差和，具体计算流程如下：
[0037]
s421)计算初始候选块sad作为sad
min
，令r＝1；
[0038]
s422)判断是否为最后一个候选块，即判断r>ca
max
是否成立，若是，则当前sad
min
为最小sad，计算结束，否则，r增加1，s＝1，sad
rs
＝0，并执行步骤s423)；
[0039]
s423)判断是否为当前候选块的最后一层，即判断s>l
max
是否成立，若是，则记当前sad
rs
为sad
min
，转入步骤s422)，否则，计算第s层sad误差sad
s
，并令sad
rs
＝sad
rs
+sad
s
，并执行步骤s424)；
[0040]
s424)判断前s层误差和sad
rs
是否大于sad
min
，若是，则执行步骤s422)，否则，s增加1，执行步骤s423)。
[0041]
本发明的有益效果为：
[0042]
本发明采用隔帧块匹配处理框架；结合心脏运动的医学知识，自适应地校准参考块的像素值以避免跟踪误差的传播与积累；对搜索区域通过三次样条插值获得1/16亚像素跟踪精度；采用搜索区域压缩、sad迭代计算提前终止等策略提升块匹配处理的计算效率。实现了一种高精度的快速3d超声心动图散斑跟踪方法，且该方法具有良好的自适应特性和鲁棒性。本发明方法适用于旨在实现心脏功能精确评价和定量诊疗的3d ste成像技术(3d speckle trackingechocardiography)。将益于提高我国冠心病、高血压等常见心血管疾病心功能定量评价的临床诊疗水平。
附图说明
[0043]
图1(a)为得到的一次1/4亚像素精度跟踪点r11及其周围8个1/4体素块示意图；
[0044]
图1(b)为1/16亚像素精度的搜索区域s12的示意图；
[0045]
图1(c)为8个1/4体素块中心点b1、b2、b3、b4、b5、b6、b7、b8的示意图；
[0046]
图1(d)为二次搜索区域s13的示意图；
[0047]
图2为3d超声心动图散斑跟踪的方法的流程示意图；
[0048]
图3为成像系统3d图像空间坐标系示意图。
具体实施方式
[0049]
下面结合说明书附图和实施例对本发明作进一步限定，但不限于此。
[0050]
实施例1
[0051]
一种3d超声心动图散斑跟踪的方法，本实施例中，某家医院的心内科超声检查室配置有国产心脏超声诊断仪和一台计算机工作站。该工作站的主要功能是基于本发明方法实现3d ste成像技术。即，利用本发明方法对3d超声心动图影像数据进行散斑跟踪，实现心肌组织的速度向量成像、应变和应变率成像，以及左心室扭转角位移、扭转力矩、扭转角速度等运动参数的测量与显像等。通过u盘等移动存储介质，临床医师将超声诊断仪采集获得的3d超声心动图影像数据转存到工作站，设3d超声心动图序列图像由n帧图像组成，分别记作i
k
，k＝1,2,...,n，其中i
k
为3d超声心动图序列图像的第k帧图像，i1为初始帧，i
n
为末帧，i1,i3,...,i
n
为奇数帧及末帧构成的序列图像；将第t帧图像进行块匹配所使用的参考块记作r
t
，t＝1,3,5,...,n，第t帧图像进行块匹配所使用的、经像素值校准后的参考块记作将在第t帧图像上获得的最佳匹配块记作m
t
，将在第t帧图像上获得的跟踪点在成像系统3d图像空间中的位置记作r
t
(x
t
,y
t
,z
t
)，其中(x
t
,y
t
,z
t
)表示跟踪点在3d图像空间中的坐标(如图3所示)，如图2所示，该方法的步骤如下：
[0052]
s1)由临床医师在初始帧i1上手动选取跟踪点；
[0053]
在3d超声心动图序列图像的初始帧i1上，由临床医师手动选取跟踪点，该跟踪点在成像系统3d图像空间中的位置为r1(x1,y1,z1)；
[0054]
s2)设定初始参数；
[0055]
设定以各帧的跟踪点为中心的心肌组织感兴趣区域的大小n
x
×
n
y
×
n
z
，其中n
x
、n
y
、n
z
分别在成像系统3d图像空间中为沿x、y、z轴方向的大小；设定搜索区域大小；设定参考块像素值校准的因子m；令当前帧序号j＝1；设定初始参考块大小，以r
j
(x
j
,y
j
,z
j
)点为中心形成初始参考块r
j
；令r
j
是指在第j帧图像进行块匹配所使用的参考块，是指第j帧图像进行块匹配所使用的、经像素值校准后的参考块，m
j
是指在第j帧图像上获得的最佳匹配块；
[0056]
s3)参考块像素值校准；
[0057]
j增加2；令r
j
＝m
j
‑2；
[0058]
对参考块r
j
的像素值进行校准，得到第j帧图像进行块匹配所使用的、经像素值校准后的参考块
[0059]
s4)进行1/16亚像素精度的块匹配计算，获得当前帧跟踪点在成像系统3d图像空
间中的位置的坐标r
j
(x
j
,y
j
,z
j
)；
[0060]
s5)若i
j
为末帧，则执行s6，否则执行s3；
[0061]
s6)计算偶数帧跟踪点在成像系统3d图像空间中的位置
[0062]
对奇数帧和末帧跟踪点在成像系统3d图像空间中的坐标，分别取其x轴、y轴、z轴坐标，通过三次样条插值计算偶数帧跟踪点在成像系统3d图像空间中的x轴、y轴、z轴坐标，获得偶数帧跟踪点在成像系统3d图像空间中的位置，进而得到所有帧跟踪点在成像系统3d图像空间中的位置；
[0063]
s7)形成并记录所有帧的心肌组织感兴趣区域
[0064]
以所有帧的跟踪点为中心，形成并记录各帧大小为n
x
×
n
y
×
n
z
的感兴趣区域，实现3d超声心动图的散斑跟踪。
[0065]
利用已记录的感兴趣区域，工作站软件通过心肌组织的速度向量成像、应变和应变率成像，以及左心室扭转运动参数测量与显像等方式或手段，可靠、便捷、高效地辅助临床医师对患者的心脏功能进行精确评价和定量诊疗。
[0066]
实施例2
[0067]
根据实施例1所述的一种3d超声心动图散斑跟踪的方法，其区别在于：
[0068]
对参考块r
j
的像素值进行校准，如式(ⅰ)所示：
[0069][0070]
式(ⅰ)中，0<m,p,q<1，m+p+q＝1，p＝0.9(1
‑
m)/(n
‑
3)2·
(j
‑
3)2，q＝1
‑
m
‑
p。
[0071]
实施例3
[0072]
根据实施例1所述的一种3d超声心动图散斑跟踪的方法，其区别在于：
[0073]
以为参考块，执行1/16亚像素精度的块匹配计算，实现过程如下：
[0074]
s41)在i
j
上对以r
j
‑2(x
j
‑2,y
j
‑2,z
j
‑2)点为中心的搜索区域进行第一次三次样条插值，使得2个相邻节点之间3个等间隔的内点被插值，获得1/4亚像素精度的搜索区域s11；
[0075]
s42)采用全搜索算法及sad(sum of absolute differences)匹配准则，并使用sad迭代计算提前终止策略，在搜索区域s11中进行一次匹配，记最小sad值对应的点为一次1/4亚像素精度跟踪点r11；
[0076]
s43)以一次1/4亚像素精度跟踪点r11为中心点，对周围8个1/4体素块(如图1(a)中a1
‑
a8为顶点的区域)进行第二次三次样条插值，继而2个相邻节点之间3个等间隔的内点被插值，获得1/16亚像素精度的搜索区域s12(如图1(b))。
[0077]
s44)通过计算1/8亚像素精度跟踪点对搜索区域s12进行压缩处理：
[0078]
计算8个1/4体素块中心点b1、b2、b3、b4、b5、b6、b7、b8(如图1(c))的sad值，取其中最小的sad值对应的点(以b5为例)为中心划定1/4体素块为二次搜索区域s13(如图1(d)以c1
‑
c7及r11为顶点的区域)；
[0079]
s44)采用全搜索算法及sad匹配准则，并使用sad迭代计算提前终止策略，进行二次匹配，最小sad值对应的点即为当前帧(第j帧)1/16亚像素精度的跟踪点r
j
(x
j
,y
j
,z
j
)，以该点为中心的块为当前帧(第j帧)上的最佳匹配块m
j
。
[0080]
本发明方法的2d变体是一种亚像素精度的快速2d超声心动图散斑跟踪方法，即本发明方法也适用于对2d超声心动图进行散斑跟踪。
[0081]
实施例4
[0082]
根据实施例1所述的一种3d超声心动图散斑跟踪的方法，其区别在于：
[0083]
步骤s42)及步骤s44)中，具体实现过程如下：
[0084]
设ca
max
为全搜索候选块总数，l
max
为全搜索候选块z轴方向层数，sad
min
为当前最小sad值，r代表当前计算第r个候选块，sad
rs
代表第r个候选块的前s层误差和，sad
s
代表第s层误差和，具体计算流程如下：
[0085]
s421)计算初始候选块sad作为sad
min
，令r＝1；
[0086]
s422)判断是否为最后一个候选块，即判断r>ca
max
是否成立，若是，则当前sad
min
为最小sad，计算结束，否则，r增加1，s＝1，sad
rs
＝0，并执行步骤s423)；
[0087]
s423)判断是否为当前候选块的最后一层，即判断s>l
max
是否成立，若是，则记当前sad
rs
为sad
min
，转入步骤s422)，否则，计算第s层sad误差sad
s
，并令sad
rs
＝sad
rs
+sad
s
，并执行步骤s424)；
[0088]
s424)判断前s层误差和sad
rs
是否大于sad
min
，若是，则执行步骤s422)，否则，s增加1，执行步骤s423)。