CN108898578B - 一种医疗图像的处理方法、装置及计算机存储介质 - Google Patents

一种医疗图像的处理方法、装置及计算机存储介质 Download PDF

Info

Publication number
CN108898578B
CN108898578B CN201810534684.1A CN201810534684A CN108898578B CN 108898578 B CN108898578 B CN 108898578B CN 201810534684 A CN201810534684 A CN 201810534684A CN 108898578 B CN108898578 B CN 108898578B
Authority
CN
China
Prior art keywords
flow field
image
medical image
target object
boundary
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
Application number
CN201810534684.1A
Other languages
English (en)
Other versions
CN108898578A (zh
Inventor
魏润杰
王虎峰
高琪
吴鹏
李飞
王巍
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Shengshi Technology Co ltd
Original Assignee
Hangzhou Shengshi Technology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hangzhou Shengshi Technology Co ltd filed Critical Hangzhou Shengshi Technology Co ltd
Priority to CN201810534684.1A priority Critical patent/CN108898578B/zh
Publication of CN108898578A publication Critical patent/CN108898578A/zh
Application granted granted Critical
Publication of CN108898578B publication Critical patent/CN108898578B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/20Image enhancement or restoration using local operators
    • G06T5/30Erosion or dilatation, e.g. thinning
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/70Denoising; Smoothing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/12Edge-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/187Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10088Magnetic resonance imaging [MRI]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20112Image segmentation details
    • G06T2207/20116Active contour; Active surface; Snakes
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20112Image segmentation details
    • G06T2207/20156Automatic seed setting
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30048Heart; Cardiac
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30101Blood vessel; Artery; Vein; Vascular

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Image Processing (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明实施例公开了一种医疗图像的处理方法、装置及计算机存储介质,该方法包括:针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出所述目标对象边界;将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界。

Description

一种医疗图像的处理方法、装置及计算机存储介质
技术领域
本发明涉及医疗图像技术领域,尤其涉及一种医疗图像的处理方法、装置及计算机存储介质。
背景技术
在临床医疗诊断中,通常会采用电子计算机断层扫描(CT,Computed Tomography)成像、核磁成像(MRI,Magnetic Resonance Image)和X光成像等医疗成像技术对人体进行透射,并最终通过记录透视的灰度信息来表征人体内组织器官的形态信息。医生通过对医疗成像中的器官组织形态信息进行分析来针对病人的病例给出诊断意见。因此,在医疗成像中,对人体器官的清晰表征以及精准形态识别是提高医生诊断准确率的重要条件。
在对器官组织完成医疗成像后,需要进行医疗图像的识别和分割两个过程。具体来说,医疗图像的识别过程是针对医疗图像判别和检测某种器官或组织的过程,而医疗图像的分割是将各种器官或组织从医疗图像中分离的过程。精确的医疗图像识别与分割过程,是后续针对器官或组织进行三维模型重建的基础,是辅助医生进行疾病诊断和治疗的重要前提,也是辅助外科手术成功的重要保障。
目前医疗图像识别和分割过程由于受到与器官和组织固有特性相关的灰度偏差场强度和噪声的影响,会导致医疗图像识别和分割的错误率较高;并且受医疗图像识别和分割粒度的影响,导致医疗图像识别和分割的精度较低。无法准确及可靠地将器官或组织从医疗图像中进行判别和分离。
发明内容
为解决上述技术问题,本发明实施例期望提供一种医疗图像的处理方法、装置及计算机存储介质;能够降低医疗图像识别和分割的错误率,并提高医疗图像识别和分割的精确度,从而可以准确及可靠地将器官或组织从医疗图像中判别和分离出来。
本发明的技术方案是这样实现的:
第一方面,本发明实施例提供了一种医疗图像的处理方法,所述方法包括:
针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;其中,所述原始图像组包括原始医疗图像、x方向流场图像、y方向流场图像和z方向流场图像;
基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出所述目标对象边界;
将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界。
第二方面,本发明实施例提供了一种医疗图像处理装置,所述装置包括:预处理部分、分割部分、修正部分和亚像素拟合部分;其中,
所述预处理部分,配置为针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;其中,所述原始图像组包括原始医疗图像、x方向流场图像、y方向流场图像和z方向流场图像;
所述分割部分,配置为基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出所述目标对象边界;
所述修正部分,配置为将所述目标对象边界的边界点按照预设的修正策略进行修正;
所述亚像素拟合部分,配置为所述修正部分将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界。
第三方面,本发明实施例提供了一种医疗图像处理装置,所述装置包括:存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序,所述计算机程序被所述处理器执行时实现如第一方面所述医疗图像处理方法的步骤。
第四方面,本发明实施例提供了一种计算机存储介质,所述计算机存储介质存储有医疗图像处理程序,当所述医疗图像处理程序被至少一个处理器执行时实现如第一方面所述的图片处理方法的步骤。
本发明实施例提供了一种医疗图像的处理方法、装置及计算机存储介质;首先通过预处理消除灰度偏差场及噪声的影响之后,对医疗图像中目标对象的边界进行分割,并且获得边界点的亚像素位置,从而获得目标对象的亚像素边界,不仅降低医疗图像识别和分割的错误率,而且提高医疗图像识别和分割的精确度,从而可以准确及可靠地将器官或组织从医疗图像中判别和分离出来。
附图说明
图1为本发明实施例提供的一种医疗图像的处理方法流程示意图;
图2(a)为本发明实施例提供的一种原始核磁图像示意图;
图2(b)为本发明实施例提供的一种x方向流场图像示意图;
图2(c)为本发明实施例提供的一种y方向流场图像示意图;
图2(d)为本发明实施例提供的一种z方向流场图像示意图;
图3为本发明实施例提供的去除灰度偏差场的流程示意图;
图4为本发明实施例提供的一种中间医疗图像示意图;
图5为本发明实施例提供的去除噪声的流程示意图;
图6为本发明实施例提供的一种预处理后的医疗图像示意图;
图7为本发明实施例提供的分割边界的流程示意图;
图8(a)为本发明实施例提供的y流场腐蚀示意图;
图8(b)为本发明实施例提供的z流场腐蚀示意图;
图8(c)为本发明实施例提供的y、z流场腐蚀图的重叠示意图;
图9(a)为本发明实施例提供的一种弓形区域二值图;
图9(b)为本发明实施例提供的一种弓形区域二值图的骨架图;
图9(c)为本发明实施例提供的一种弓形区域二值图的骨架及种子点图;
图10为本发明实施例提供的一种边界示意图;
图11为本发明实施例提供的修正边界点的流程示意图;
图12(a)为本发明实施例提供的一种基于边界点划分预设范围的示意图;
图12(b)为本发明实施例提供的一种基于边界点构成方框的示意图;
图12(c)为本发明实施例提供的一种基于边界点进行修正拟合的示意图;
图13为本发明实施例提供的获取亚像素边界的流程示意图;
图14为本发明实施例提供的一种像素边界点以及亚像素边界点示意图;
图15为本发明实施例提供的一种三维重构体示意图;
图16为本发明实施例提供的一种医疗图像处理装置组成示意图;
图17为本发明实施例提供的另一种医疗图像处理装置组成示意图;
图18为本发明实施例提供的医疗图像处理装置的具体硬件结构示意图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述。
参见图1,其示出了本发明实施例提供的一种医疗图像的处理方法,所述方法包括:
S101:针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;
其中,所述原始图像组包括原始医疗图像、x方向流场图像、y方向流场图像和z方向流场图像;
S102:基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于目标对象轮廓从所述预处理后的医疗图像中分割出所述目标对象边界;
S103:将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界。
通过图1所示的技术方案,首先通过预处理消除灰度偏差场及噪声的影响之后,对医疗图像中目标对象的边界进行分割,并且获得边界点的亚像素位置,从而获得目标对象的亚像素边界,不仅降低医疗图像识别和分割的错误率,而且提高医疗图像识别和分割的精确度,从而可以准确及可靠地将器官或组织从医疗图像中判别和分离出来。
对于图1所示的技术方案,需要说明的是,在采集原始图像组的过程中,本发明实施例以核磁成像心脏区域的医疗图像为举例,可以针对一个心跳周期,如25个时刻内的30个连续矢状面进行采集,并基于目前核磁成像数据不但包括器官成像信息,还包括血液流动的三维流场信息,可以得出:对任意一个时刻且在该时刻的每个矢状面都包含4幅图像,从而组成一个原始图像组,以第ti时刻第sj个矢状面为例,参见图2,原始图像组中的4幅图像可以包括:图2(a)所示的第ti时刻第sj个矢状面的核磁图像、图2(b)所示的x方向流场图像、图2(c)所示的y方向流场图像和图2(d)所示的z方向流场图像,且这4幅图像每幅的分辨率均可设置为r0×c0(在本实施例中优选为r0=256,c0=256),总计能够得到25×30×4=3000张图像。通过图1所示的技术方案,可以将从25个时刻中任意给定一个时刻(即第ti时刻),精确地分割出该时刻的每个矢状面的核磁图像中的动脉血管。基于此,就可以将该时刻(即第ti时刻)的所有30个矢状面上的动脉血管重构成为一个高精度的三维动脉血管。可以理解地,图1所示的技术方案除了可以应用于核磁成像心脏区域的医疗图像以外,还可以应用于其他成像技术针对其他器官或组织进行采集的医疗图像,本实施例对此不作具体限定和赘述。
针对图1中所述的S101,在一种可能的实现方式中,所述针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像,包括:
将所述原始医疗图像按照预设的去除策略去除灰度偏差场,获得中间医疗图像;
将所述中间医疗图像按照预设的去噪策略去除噪声后,获得所述预处理后的医疗图像。
对于该实现方式,具体来说,参见图3,所述将所述原始医疗图像按照预设的去除策略去除灰度偏差场,获得中间医疗图像,可以包括:
S301:设定原始医疗图像与灰度偏差场之间满足式1所示的对应关系:
I(x)=B(x)J(x)+nσ(x) (1)
其中,I(x)为所述原始医疗图像在像素点x=(x1,x2)的图像灰度值、B(x)为所述原始医疗图像在像素点x=(x1,x2)的灰度偏差场的值、J(x)为所述原始医疗图像在像素点x=(x1,x2)的真实图像灰度值、nσ(x)为所述原始医疗图像在像素点x=(x1,x2)的噪声值;
需要说明的是,上述I(x),B(x),J(x),nσ(x)都分别是针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像来说的,并且均是非负数。
S302:基于灰度偏差场的均匀变化性,将B(x)通过Nb个光滑的基函数的线性组合进行表示,获得式2所示的表达式:
Figure GDA0002651746540000061
其中,gk(x)是第k个光滑基函数,且1≤k≤Nb,wk是gk(x)的权重系数,
Figure GDA0002651746540000062
在本步骤中,设定灰度偏差场是变化均匀的,因此,B(x)是一个关于x的光滑函数,因此就可以通过光滑的基函数的线性组合进行表示,在本实施例中,优选可以设定Nb=10,分别以下式所示:
Figure GDA0002651746540000071
S303:将J(x)划分为N个部分,且每个部分中所有像素点的灰度值都相同,则获得式3所示的表达式:
Figure GDA0002651746540000072
其中,cq是第q部分的灰度值,uq(x)是cq的从属函数,c=(c1,c2,…,cN)T,u=(u1,u2,…,uN)T
需要说明的是,可以设定J(x)由N个部分构成,在本实施例中,N=3,即设定真实图像J(x)包含有血流的部分、背景部分和其他部分;并且每个部分当中的所有像素点的灰度值都相同。而且在本发明实施例中,从属函数uq(x)的具体取值可以优选为:
Figure GDA0002651746540000073
S304:基于式1建立能量函数Fe(B,J)=∫Ω|I(x)-B(x)J(x)|2dx,并基于式2和式3获得式4所示的能量函数的分解形式:
Figure GDA0002651746540000081
其中,|·|是绝对值符号,Ω是所述像素点x=(x1,x2)的坐标取值范围,Ωq是所述像素点x=(x1,x2)在第q部分的坐标取值范围;
需要说明的是,在获得能量函数的分解形式后,可以通过后续S305至S306求出使能量函数达到最小值时对应的u,c,w的值。
S305:基于式4分别对u,c,w取偏导数并赋值为零,获得关于u,c,w的更新函数;
需要说明的是,在本实施例中,关于c的更新函数,具体可以通过下式获得:
Figure GDA0002651746540000082
关于w的更新函数,具体可以通过
Figure GDA0002651746540000083
获得,其中,
Figure GDA0002651746540000084
关于u的更新函数,具体可以通过下式获得:
Figure GDA0002651746540000085
可以理解地,
Figure GDA0002651746540000086
为“推导出”符号,表示由该符号的左边等式可以推导出该符号的右边等式。
S306:将u,c进行初始化后,按照初始化后的u,c对关于u,c,w的更新函数分别进行迭代更新,直至满足预设的收敛条件时,获得满足预设的收敛条件所对应的w*,c*,u*
需要说明的是,在本实施例中,c的每个元素初始值采用0到1的随机数,u的每个元素的初始值设为1;而预设的收敛条件可以设置为|c(ite)-c(ite-1)|<10-3,其中,c(ite)为第ite次更新后c的值。
S307:根据w*,c*,u*以及式2和式3,分别获得B*(x)=(w*)Tg*,J*(x)=(c*)Tu*
S308:基于式1以及B*(x)=(w*)Tg*,J*(x)=(c*)Tu*,按照式5获得去除灰度偏差场后的中间医疗图像I1(x):
Figure GDA0002651746540000091
其中,ε为预设数值,在本实施例中,优选取ε为10-2
通过图3所示的技术方案,从而能够去除原始医疗图像中的灰度偏差场,以前述图2(a)所示的核磁图像,通过图3所示的技术方案,可以获得图4所示的中间医疗图像。那么在获得中间医疗图像后,就需要针对噪声进行消除,因此,具体来说,参见图5,所述将所述中间医疗图像按照预设的去噪策略去除噪声后,获得所述预处理后的医疗图像,可以包括:
S501:将第sj个矢状面在一个预设周期内的每一时刻的原始医疗图像对应的中间医疗图像转化为行向量,并将所有行向量按照所述预设周期的时间顺序组成数据矩阵X;
具体来说,以前述针对一个心跳周期,如25个时刻内的30个连续矢状面进行采集核磁图像为例,第sj个矢状面上,所有25个时刻的每一幅去除灰度偏差场后的中间核磁图像I1,转化成为一个行向量,因此,每个行向量的列数是256×256=65365,然后将25个行向量按采集的时间顺序罗列起来,组成一个数据矩阵X,大小为25×65365。
S502:对XT做奇异值分解,得到XT=UEVT
其中,符号T代表矩阵转置;
Figure GDA0002651746540000101
σγ是XXT的特征值的平方根,且在对角线上按从大到小排列;V的列是XXT的特征向量;U=XTVE-1;V,U的各自的列向量标准正交;
S503:对奇异值分解后的XT按照式6进行矩阵分解:
Figure GDA0002651746540000102
其中,E1是E的前n1行与n1列组成的矩阵,大小是n1×n1,且1≤n1≤M;M表示所述预设周期内的采样时刻点的数目,E2是E的后(M-n1)行与(M-n1)列,大小是(M-n1)×(M-n1);Ol是(M-n1)×n1的零矩阵;Ou是n1×(M-n1)的零矩阵;U1是由U的前n1列组成的矩阵,U2是由U的其余列组成的矩阵,V1是由V的前n1列组成的矩阵,V2是由V的其余列组成的矩阵;
在本实施例中,所述预设周期为一个心跳周期,该周期内有25个采样时刻点,因此,E2是E的后(25-n1)行与(25-n1)列,大小是(25-n1)×(25-n1);Ol是(25-n1)×n1的零矩阵;Ou是n1×(25-n1)的零矩阵;U1是由U的前n1列组成的矩阵,U2是由U的其余列组成的矩阵,V1是由V的前n1列组成的矩阵,V2是由V的其余列组成的矩阵;
S504:在设定Y*=U1E1V1 T后,根据X*=(Y*)T获取n1模态的数据矩阵X*
S505:将X*的每一行的行向量复原分别为对第sj个矢状面在一个预设周期内每一时刻预处理后的医疗图像;
S506:将所有预处理后的医疗图像的第ti个图像确定为第ti时刻第sj个矢状面预处理后的医疗图像I2
对于S504至S506来说,具体来说,X*=(Y*)T是n1模态的数据矩阵。将X*的25行的每一行的行向量复原为256×256的矩阵,所得矩阵就是第sj个矢状面上,每一个时刻的去除噪音后的核磁图像,共计25张。
可以理解地,参见图6,将图4所示的去除完灰度偏差场后的中间核磁图像按照图5所示的技术方案去除噪音后,能够更加清楚的在图6中呈现目标对象,即心脏血管的区域。
在对图像预处理完毕后,就可以针对目标对象的边界进行区域分割,区域分割的过程是将图像中的特定组织或者器官判别、检测和分离的过程。在本实施例中,进行分割的过程利用了流场图像的信息,基于此,对于图1所示的技术方案,参见图7,S102所述的基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于目标对象轮廓从所述预处理后的医疗图像中分割出所述目标对象边界,可以包括:
S701:基于预设个数nt阈值的大津法分别获取第一流场图像对应的第一流场(nt+1)值图以及第二流场图像对应的第二流场(nt+1)值图;
其中,所述第一流场图像为所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中的任一方向流场图像;所述第二流场图像为除第一流场图像以外的其他任意方向流场图像;
需要说明的是,在本实施例中,第一流场图像为y方向流场图像,第二流场图像为z方向流场图像;可以设定nt=3,根据大津法(Otsu)作用在心脏收缩时刻tS,第sj矢状面的y、z流场图像上,分别得到y、z流场(nt+1)值图,在本实施例中优选为4值图。
S702:分别对第一流场(nt+1)值图以及第二流场(nt+1)值图进行图像腐蚀,得到第一流场腐蚀图以及第二流场腐蚀图;
S703:分别将所述第一流场腐蚀图以及所述第二流场腐蚀图中灰度值取最大值的区域确定为所述第一流场腐蚀图的强信号区域以及所述第二流场腐蚀图的强信号区域;
对于S702和S703,具体来说,本实施例针对y、z流场(nt+1)值图采用了半径为2的圆盘为腐蚀的结构元素对象进行腐蚀,得到y、z流场腐蚀图,分别如图8(a)、8(b)所示,在所得到的y、z流场腐蚀图中,取灰度值最大值的区域为y、z流场图像的强信号区域。
S704:将第一流场腐蚀图的强信号区域以及第二流场腐蚀图的强信号区域之间的连通区域进行合并后,对合并后的连通区域进行光滑运算,获得目标对象区域二值图;
具体来说,在本实施例中,以目标对象为心脏血管为例,目标对象区域也就是心脏血管的弓形区域。可以寻找y、z流场图像强信号区域的连通区域,即判断y、z流场图像信号较强区域的重叠性,如图8(c),并将连通区域合并,非连通区域删除;随后,可以在合并后的连通区域上做图像光滑运算,本实施例优选采用图像膨胀,得到一个如图9(a)所示的弓形区域二值图。
S705:基于所述目标对象区域二值图的骨架后,以所述目标对象区域图像中心点为顶点,确定所述目标对象区域二值图骨架上的种子点;
需要说明的是,仍然以目标对象区域是心脏血管的弓形区域为例,弓形区域二值图的骨架如图9(b)所示,以弓形区域图像的中心点为顶点,以正向水平方向为起始方向,以一个固定角度,本实施例为45度,逆时针旋转360度,并沿着每一个经过方向做射线。射线与图9(b)所示的弓形区域二值图骨架相交的点,记为种子点。如图9(c)。
可以理解地,如果目标对象区域图像大小为100×80,那么其图像中心点位置为50.5和40.5。
S706:基于所述种子点对所述预处理后的医疗图像进行轮廓分割,获得每个种子点对应的最终轮廓;
需要说明的是,本实施例可以采用主动轮廓的方法在第ti时刻,第sj个矢状面的预处理后的核磁图像I2上做动脉血管分割,具体来说,以图9(c)中每一个种子点为圆心,以特定的半径长度(本实施例选取半径为10)的圆环作为初始轮廓,做特定次数的迭代(本实施例选取迭代次数为150),从而获得通过各个种子点生成的最终轮廓。
S707:将所有种子点对应的最终轮廓进行合并,获得预处理后的医疗图像中目标对象的边界。
需要说明的是,通过图7所示的方案,可以获得图10所示的边界示意图。而该边界示意图由于受到噪音、灰度差别和器官复杂性等原因的影响,部分边界区域会产生偏差,因此,需要对偏差进行修正。所以,参见图11,S103中所述的将所述目标对象边界的边界点按照预设的修正策略进行修正,可以包括:
S1101:按照预设的排列顺序对所有边界点进行排列及编号,并获取各边界点的曲率;
具体来说,针对图10所示的边界示意图,在本实施例中,将边界点沿着逆时针排列编号,并计算出所有边界点的曲率,找出所有曲率当中取得最大值的边界点,并在图中用星号标出。
S1102:将所有边界点中曲率最大的边界点进行标记后执行以下步骤:
步骤a,以任一曲率最大的边界点Q0为中心,在预设的范围内搜索排列在Q0之后且与Q0的编号差别小于预设常数的所有最大曲率值点Q=(Q1,Q2,...,Qs);
具体来说,如图12(a)所示,以任一取得最大曲率值的一点Q0为中心,做一个特定长度半径的圆(本实施例采用半径为10),在该圆内部寻找编号在Q0之后且与Q0的编号差别小于预设常数(本实施例该常数为40)的取得最大曲率值的所有点Q=(Q1,Q2,...,Qs)。
步骤b,以Q0为第一个点,以Q中每一个点Qκ∈Q为第二个点,Q0Qκ为一边按照预设的方向组成方框;
具体来说,如图12(b)所示,以Q0为第一个点,Qκ∈Q为第二个点,Q0Qκ为一边,沿着顺时针的方向,依次找出第三个点与第四个点,从而构成一个正方形方框。
步骤c,获取所述方框的分形维数;
具体来说,本实施例中对于步骤c的分形维数可以采用计盒维数表示,也就是说,获取所述方框的计盒维数,并将所述方框的计盒维数作为所述方框的分形维数。具体获取计盒维数的过程如下:在上述正方形方框内,将属于边界内的点设置为1,其他点为0;随后用不同边长li(本实施例li={1,2,3,4})的小正方形去覆盖该正方形边框,记录覆盖该正方形边框内的所有值为1的点所需的最少小正方形的个数;最后以点(log(li),log(Nli)),li={1,2,3,4}拟合出直线,并求出斜率,那么该斜率的相反数就是计盒维数。
步骤d,当基于所述分形维数确定Q0Qκ两点之间的边界异常时,在Q0Qκ两点之间生成拟合曲线作为修正后的边界;
需要说明的是,当分形维数通过计盒维数表示时,可以通过将计盒维数与预设的阈值进行比较,如果大于则确定Q0Qκ两点之间的边界异常。本实施例中,所述预设的阈值优选为1.6586;
此外,当Q0Qκ两点之间的边界异常时,可以拟合一条曲线,连接Q0Qκ两点,作为修正后的边界,在本实施例优选采用的是将Q0Qκ两点用直线连接进行拟合,如图12(c)所示。
S1103:针对每个曲率最大的边界点执行步骤a至步骤d,获得修正后的边界点。
在通过图11所示方案对异常的边界进行修正后,考虑到这些边界为像素级别的边界(即像素位置只可以取整数),粒度较大,在描述一些现实中的图形时,往往会出现失真。因此,本实施例会通过获取每个边界点的亚像素位置(即像素位置可以取实数)来描述目标对象的边界,从而会更加准确和可靠地描述边界。基于此,参见图13,S103中所述的获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界,包括:
S1300:在对所述目标对象边界的边界点修正完毕后,针对每个边界点执行以下步骤:
可以理解地,S1301至S1310针对单个边界点获取对应的亚像素位置,所以就需要针对目标对象的每一个边界点按照S1301至S1310获取对应的亚像素为止,从而获得亚像素边界。
S1301:设定在以边界点(x0,y0)为原点的坐标轴下,直线y=a+bx为真实边界线;
S1302:以所述边界点(x0,y0)为中心确定大小为5×3个方格的周边范围;
S1303:通过所述直线y=a+bx将所述周边范围划分为两个区域;其中,所述直线y=a+bx上方的第一区域的平均灰度值为A1,所述直线y=a+bx下方的第二区域的平均灰度值为A2
需要说明的是,在具体实施过程中,设定点(x0,y0)是预处理后的医疗图像I2中的一个边界点,在点(x0,y0)附近确定如图14所示的5×3方格,将I2图像所有点做坐标平移,使得点(x0,y0)新坐标为(0,0),又设定直线y=a+bx是真实边界,将图像分为两个区域,直线左上方是区域1,其平均灰度假设为A1,直线右下方是区域2,其平均灰度假设为A2
S1304:设定所述周边范围内每个点(x',y')满足
Figure GDA0002651746540000151
其中,h是像素边长,Px',y'是(x',y')位置的像素在所述直线y=a+bx下方所占的面积;
S1305:按照式7分别获取所述周边范围的左、中、右列在所述直线y=a+bx下方的面积L,M,R:
Figure GDA0002651746540000152
S1306:按照式8分别获取所述周边范围的左、中、右列的灰度值之和SL,SM,SR
Figure GDA0002651746540000161
S1307:按照式9对所述第一区域的平均灰度值A1和所述第二区域的平均灰度值A2进行拟合:
Figure GDA0002651746540000162
需要说明的是,S1307中,可以设定区域1和区域2的灰度值分别可以用左上边角三个和右下边角三个的像素平均灰度值分别去近似,从而得到式9。
S1308:基于式7、式8和式9获得所述直线y=a+bx中的参数
Figure GDA0002651746540000163
以及
Figure GDA0002651746540000164
S1309:在以边界点(x0,y0)为原点的坐标轴下,将所述原点的亚像素位置确认为(0,a);
S1310:根据所述原点的亚像素位置(0,a)确定所述边界点(x0,y0)的亚像素坐标在预处理后的医疗图像中为(x0,y0+a);
可以理解地,以图14所示的5×3方格图中的亚像素位置认为是(0,a)。将平移后的坐标进行恢复,从而可以得到I2中该边界点(x0,y0)的亚像素坐标为(x0,y0+a)。
上述方案,是针对第ti时刻第sj个矢状面的原始图像组进行处理,而在具体实现过程中,需要对第ti时刻的所有30张矢状面上的图像组,分别按照上述方案进行处理,从而得到第ti时刻所有矢状面的原始图像组对应目标对象的亚像素边界,因此,就能够进行三维重构。基于此,在一种可能的实现方式中,图1所示的技术方案还可以包括:
在获取到第ti时刻每个矢状面的原始图像组对应目标对象的亚像素边界后,基于所述目标对象的所有亚像素边界进行三维组装,获得所述目标对象的原始三维体;
将所述原始三维体进行平滑后,获得所述目标对象的三维重构体。
以前述举例来说,在得到第ti时刻30个矢状面的亚像素边界后,可以将每个矢状面上具有若干亚像素边界点,以该矢状面的顺序数为x坐标(即第sj个矢状面的x坐标是sj),以亚像素边界点的坐标为z和y坐标,可以生成三维空间中的点云;接着可以采用泊松表面重构算法,获得光滑、闭合的三维曲面,该三维曲面是对点云的最佳拟合曲面,同时该三维曲面所包含的实体就是三维重构体,如图15所示。由于通过亚像素边界点进行目标对象的重构,从而比像素级别的边界点具有更加精准的描述性,提高了识别和分割的精确度。
为实现前述图1所示的医疗图像的处理方法,参见图16,其示出了本发明实施例提供的一种医疗图像处理装置160,该装置可以包括:预处理部分1601、分割部分1602、修正部分1603和亚像素拟合部分1604;其中,
所述预处理部分1601,配置为针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;其中,所述原始图像组包括原始医疗图像、x方向流场图像、y方向流场图像和z方向流场图像;
所述分割部分1602,配置为基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出所述目标对象边界;
所述修正部分1603,配置为将所述目标对象边界的边界点按照预设的修正策略进行修正;
所述亚像素拟合部分1604,配置为所述修正部分1603将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界。
在上述方案中,所述预处理部分1601,配置为:
将所述原始医疗图像按照预设的去除策略去除灰度偏差场,获得中间医疗图像;
将所述中间医疗图像按照预设的去噪策略去除噪声后,获得所述预处理后的医疗图像。
在上述方案中,所述预处理部分1601,配置为:
设定原始医疗图像与灰度偏差场之间满足式1所示的对应关系;
基于灰度偏差场的均匀变化性,将B(x)通过Nb个光滑的基函数的线性组合进行表示,获得式2所示的表达式;
将J(x)划分为N个部分,且每个部分中所有像素点的灰度值都相同,则获得式3所示的表达式;
基于式1建立能量函数Fe(B,J)=∫Ω|I(x)-B(x)J(x)|2dx,并基于式2和式3获得式4所示的能量函数的分解形式;
基于式4分别对u,c,w取偏导数并赋值为零,获得关于u,c,w的更新函数;
将u,c进行初始化后,按照初始化后的u,c对关于u,c,w的更新函数分别进行迭代更新,直至满足预设的收敛条件时,获得满足预设的收敛条件所对应的w*,c*,u*
根据w*,c*,u*以及式2和式3,分别获得B*(x)=(w*)Tg*,J*(x)=(c*)Tu*
基于式1以及B*(x)=(w*)Tg*,J*(x)=(c*)Tu*,按照式5获得去除灰度偏差场后的中间医疗图像I1(x)。
在上述方案中,所述预处理部分1601,配置为:
将第sj个矢状面在一个预设周期内的每一时刻的原始医疗图像对应的中间医疗图像转化为行向量,并将所有行向量按照所述预设周期的时间顺序组成数据矩阵X;
对XT做奇异值分解,得到XT=UEVT
对奇异值分解后的XT按照式6进行矩阵分解;
在设定Y*=U1E1V1 T后,根据X*=(Y*)T获取n1模态的数据矩阵X*
将X*的每一行的行向量复原分别为对第sj个矢状面在一个预设周期内每一时刻预处理后的医疗图像;
将所有预处理后的医疗图像的第ti个图像确定为第ti时刻第sj个矢状面预处理后的医疗图像I2
在上述方案中,所述分割部分1602,配置为包括:
基于预设个数nt阈值的大津法分别获取第一流场图像对应的第一流场(nt+1)值图以及第二流场图像对应的第二流场(nt+1)值图;其中,所述第一流场图像为所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中的任一方向流场图像;所述第二流场图像为除第一流场图像以外的其他任意方向流场图像;
分别对第一流场(nt+1)值图以及第二流场(nt+1)值图进行图像腐蚀,得到第一流场腐蚀图以及第二流场腐蚀图;
分别将所述第一流场腐蚀图以及所述第二流场腐蚀图中灰度值取最大值的区域确定为所述第一流场腐蚀图的强信号区域以及所述第二流场腐蚀图的强信号区域;
将所述第一流场腐蚀图的强信号区域以及所述第二流场腐蚀图的强信号区域之间的连通区域进行合并后,对合并后的连通区域进行光滑运算,获得目标对象区域二值图;
基于所述目标对象区域二值图的骨架后,以所述目标对象区域中心点为顶点,确定所述目标对象区域二值图骨架上的种子点;
基于所述种子点对所述预处理后的医疗图像进行轮廓分割,获得每个种子点对应的最终轮廓;
将所有种子点对应的最终轮廓进行合并,获得预处理后的医疗图像中目标对象的边界。
在上述方案中,所述修正部分1603,配置为:
按照预设的排列顺序对所有边界点进行排列及编号,并获取各边界点的曲率;
将所有边界点中曲率最大的边界点进行标记后执行以下步骤:
步骤a,以任一曲率最大的边界点Q0为中心,在预设的范围内搜索排列在Q0之后且与Q0的编号差别小于预设常数的所有最大曲率值点Q=(Q1,Q2,...,Qs);
步骤b,以Q0为第一个点,以Q中每一个点Qκ∈Q为第二个点,Q0Qκ为一边按照预设的方向各组成方框;
步骤c,获取所述方框的分形维数;
步骤d,当基于所述分形维数确定Q0Qκ两点之间的边界异常时,在Q0Qκ两点之间生成拟合曲线作为修正后的边界;
针对每个曲率最大的边界点执行步骤a至步骤d,获得修正后的边界点。
在上述方案中,所述修正部分1603,配置为:
获取所述方框的计盒维数,并将所述方框的计盒维数作为所述方框的分形维数。
在上述方案中,所述亚像素拟合部分1604,配置为:
在对所述目标对象边界的边界点修正完毕后,针对每个边界点执行以下步骤:
设定在以边界点(x0,y0)为原点的坐标轴下,直线y=a+bx为真实边界线;
以所述边界点(x0,y0)为中心确定大小为5×3个方格的周边范围;
通过所述直线y=a+bx将所述周边范围划分为两个区域;其中,所述直线y=a+bx上方的第一区域的平均灰度值为A1,所述直线y=a+bx下方的第二区域的平均灰度值为A2
设定所述周边范围内每个点(x',y')满足
Figure GDA0002651746540000201
其中,h是像素边长,Px',y'是(x',y')位置的像素在所述直线y=a+bx下方所占的面积;
按照式7分别获取所述周边范围的左、中、右列在所述直线y=a+bx下方的面积L,M,R;
按照式8分别获取所述周边范围的左、中、右列的灰度值之和SL,SM,SR
按照式9对所述第一区域的平均灰度值A1和所述第二区域的平均灰度值A2进行拟合;
基于式7、式8和式9获得所述直线y=a+bx中的参数
Figure GDA0002651746540000211
以及
Figure GDA0002651746540000212
在以边界点(x0,y0)为原点的坐标轴下,将所述原点的亚像素位置确认为(0,a);
根据所述原点的亚像素位置(0,a)确定所述边界点(x0,y0)的亚像素坐标在预处理后的医疗图像中为(x0,y0+a)。
在上述方案中,参见图17,所述装置160还包括:三维重构部分1605,配置为:
在获取到第ti时刻每个矢状面的原始图像组对应目标对象的亚像素边界后,基于所述目标对象的所有亚像素边界进行三维组装,获得所述目标对象的原始三维体;以及,
将所述原始三维体进行平滑后,获得所述目标对象的三维重构体。
可以理解地,在本实施例中,“部分”可以是部分电路、部分处理器、部分程序或软件等等,当然也可以是单元,还可以是模块也可以是非模块化的。
另外,在本实施例中的各组成部分可以集成在一个处理单元中,也可以是各个单元单独物理存在,也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现,也可以采用软件功能模块的形式实现。
所述集成的单元如果以软件功能模块的形式实现并非作为独立的产品进行销售或使用时,可以存储在一个计算机可读取存储介质中,基于这样的理解,本实施例的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来,该计算机软件产品存储在一个存储介质中,包括若干指令用以使得一台计算机设备(可以是个人计算机,服务器,或者网络设备等)或processor(处理器)执行本实施例所述方法的全部或部分步骤。而前述的存储介质包括:U盘、移动硬盘、只读存储器(ROM,Read Only Memory)、随机存取存储器(RAM,Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
因此,本实施例提供了一种计算机可读存储介质,该计算机可读存储介质存储有医疗图像的处理程序,所述医疗图像的处理程序被至少一个处理器执行时实现上述图1所示技术方案所述的医疗图像的处理方法的步骤。
基于上述医疗图像处理装置160以及计算机可读存储介质,参见图18,其示出了本发明实施例提供的一种医疗图像处理装置160的具体硬件结构,可以包括:存储器1801、处理器1802及存储在所述存储器1801上并可在所述处理器1802上运行的计算机程序1803。各个组件通过总线系统1804耦合在一起。可理解,总线系统1804用于实现这些组件之间的连接通信。总线系统1804除包括数据总线之外,还包括电源总线、控制总线和状态信号总线。但是为了清楚说明起见,在图18中将各种总线都标为总线系统1804。
所述处理器1802,配置为执行所述计算机程序1803以实现以下步骤:
针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;其中,所述原始图像组包括原始医疗图像、x方向流场图像、y方向流场图像和z方向流场图像;
基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出所述目标对象边界;
将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界。
可以理解,本发明实施例中的存储器1801可以是易失性存储器或非易失性存储器,或可包括易失性和非易失性存储器两者。其中,非易失性存储器可以是只读存储器(Read-Only Memory,ROM)、可编程只读存储器(Programmable ROM,PROM)、可擦除可编程只读存储器(Erasable PROM,EPROM)、电可擦除可编程只读存储器(Electrically EPROM,EEPROM)或闪存。易失性存储器可以是随机存取存储器(Random Access Memory,RAM),其用作外部高速缓存。通过示例性但不是限制性说明,许多形式的RAM可用,例如静态随机存取存储器(Static RAM,SRAM)、动态随机存取存储器(Dynamic RAM,DRAM)、同步动态随机存取存储器(Synchronous DRAM,SDRAM)、双倍数据速率同步动态随机存取存储器(Double DataRate SDRAM,DDRSDRAM)、增强型同步动态随机存取存储器(Enhanced SDRAM,ESDRAM)、同步连接动态随机存取存储器(Synchlink DRAM,SLDRAM)和直接内存总线随机存取存储器(Direct Rambus RAM,DRRAM)。本文描述的系统和方法的存储器1801旨在包括但不限于这些和任意其它适合类型的存储器。
而处理器1802可能是一种集成电路芯片,具有信号的处理能力。在实现过程中,上述方法的各步骤可以通过处理器1802中的硬件的集成逻辑电路或者软件形式的指令完成。上述的处理器1802可以是通用处理器、数字信号处理器(Digital Signal Processor,DSP)、专用集成电路(Application Specific Integrated Circuit,ASIC)、现成可编程门阵列(Field Programmable Gate Array,FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。结合本发明实施例所公开的方法的步骤可以直接体现为硬件译码处理器执行完成,或者用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器,闪存、只读存储器,可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器1801,处理器1802读取存储器1801中的信息,结合其硬件完成上述方法的步骤。
可以理解的是,本文描述的这些实施例可以用硬件、软件、固件、中间件、微码或其组合来实现。对于硬件实现,处理单元可以实现在一个或多个专用集成电路(ApplicationSpecific Integrated Circuits,ASIC)、数字信号处理器(Digital Signal Processing,DSP)、数字信号处理设备(DSP Device,DSPD)、可编程逻辑设备(Programmable LogicDevice,PLD)、现场可编程门阵列(Field-Programmable Gate Array,FPGA)、通用处理器、控制器、微控制器、微处理器、用于执行本申请所述功能的其它电子单元或其组合中。
对于软件实现,可通过执行本文所述功能的模块(例如过程、函数等)来实现本文所述的技术。软件代码可存储在存储器中并通过处理器执行。存储器可以在处理器中或在处理器外部实现。
具体来说,所述处理器1802,还配置为运行所述计算机程序时,执行前述图1所示的方法步骤,这里不再进行赘述
需要说明的是:本发明实施例所记载的技术方案之间,在不冲突的情况下,可以任意组合。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应以所述权利要求的保护范围为准。

Claims (13)

1.一种医疗图像的处理方法,其特征在于,所述方法包括:
针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;其中,所述原始图像组包括原始医疗图像、x方向流场图像、y方向流场图像和z方向流场图像;
基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出目标对象边界;
将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到目标对象的亚像素边界;
其中,所述基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出目标对象边界,包括:对所述任意两个方向流场图进行图像腐蚀,得到第一流场腐蚀图以及第二流场腐蚀图,且将所述第一流场腐蚀图以及所述第二流场腐蚀图中灰度值取最大值的区域确定为第一流场腐蚀图的强信号区域以及第二流场腐蚀图的强信号区域,将所述第一流场腐蚀图的强信号区域以及所述第二流场腐蚀图的强信号区域之间的连通区域进行合并后,对合并后的连通区域进行光滑运算,获得目标对象区域二值图;将所述第一流场腐蚀图的强信号区域以及所述第二流场腐蚀图的强信号区域之间的连通区域进行合并后,对合并后的连通区域进行光滑运算,获得目标对象区域二值图;基于所述目标对象区域二值图的骨架后,以目标对象区域中心点为顶点,确定所述目标对象区域二值图骨架上的种子点;基于所述种子点对所述预处理后的医疗图像进行轮廓分割,获得每个种子点对应的最终轮廓。
2.根据权利要求1所述的方法,其特征在于,所述针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像,包括:
将所述原始医疗图像按照预设的去除策略去除灰度偏差场,获得中间医疗图像;
将所述中间医疗图像按照预设的去噪策略去除噪声后,获得所述预处理后的医疗图像。
3.根据权利要求2所述的方法,其特征在于,所述将所述原始医疗图像按照预设的去除策略去除灰度偏差场,获得中间医疗图像,包括:
设定原始医疗图像与灰度偏差场之间满足式1所示的对应关系:
I(x)=B(x)J(x)+nσ(x) (1)
其中,I(x)为所述原始医疗图像在像素点x=(x1,x2)的图像灰度值、B(x)为所述原始医疗图像在像素点x=(x1,x2)的灰度偏差场的值、J(x)为所述原始医疗图像在像素点x=(x1,x2)的真实图像灰度值、nσ(x)为所述原始医疗图像在像素点x=(x1,x2)的噪声值;
基于灰度偏差场的均匀变化性,将B(x)通过Nb个光滑的基函数的线性组合进行表示,获得式2所示的表达式:
Figure FDA0002651746530000021
其中,gk(x)是第k个光滑基函数,且1≤k≤Nb,wk是gk(x)的权重系数,
Figure FDA0002651746530000022
T为转置运算符;
将J(x)划分为N个部分,且每个部分中所有像素点的灰度值都相同,则获得式3所示的表达式:
Figure FDA0002651746530000023
其中,cq是第q部分的灰度值,uq(x)是cq的从属函数,c=(c1,c2,…,cN)T,u=(u1,u2,…,uN)T
基于式1建立能量函数Fe(B,J)=∫Ω|I(x)-B(x)J(x)|2dx,并基于式2和式3获得式4所示的能量函数的分解形式:
Figure FDA0002651746530000031
其中,|·|是绝对值符号,Ω是所述像素点x=(x1,x2)的坐标取值范围,Ωq是所述像素点x=(x1,x2)在第q部分的坐标取值范围;
基于式4分别对u,c,w取偏导数并赋值为零,获得关于u,c,w的更新函数;
将u,c进行初始化后,按照初始化后的u,c对关于u,c,w的更新函数分别进行迭代更新,直至满足预设的收敛条件时,获得满足预设的收敛条件所对应的w*,c*,u*
根据w*,c*,u*以及式2和式3,分别获得B*(x)=(w*)Tg*,J*(x)=(c*)Tu*
基于式1以及B*(x)=(w*)Tg*,J*(x)=(c*)Tu*,按照式5获得去除灰度偏差场后的中间医疗图像I1(x):
Figure FDA0002651746530000032
其中,ε为预设数值。
4.根据权利要求2所述的方法,其特征在于,所述将所述中间医疗图像按照预设的去噪策略去除噪声后,获得所述预处理后的医疗图像,包括:
将第sj个矢状面在一个预设周期内的每一时刻的原始医疗图像对应的中间医疗图像转化为行向量,并将所有行向量按照所述预设周期的时间顺序组成数据矩阵X;
对XT做奇异值分解,得到XT=UEVT;其中,符号T代表矩阵转置;
Figure FDA0002651746530000041
σγ是XXT的特征值的平方根,且在对角线上按从大到小排列;V的列是XXT的特征向量;U=XTVE-1;V,U的各自的列向量标准正交;
对奇异值分解后的XT按照式6进行矩阵分解:
Figure FDA0002651746530000042
其中,E1是E的前n1行与n1列组成的矩阵,大小是n1×n1,且1≤n1≤M;M表示所述预设周期内的采样时刻点的数目,E2是E的后(M-n1)行与(M-n1)列,大小是(M-n1)×(M-n1);Ol是(M-n1)×n1的零矩阵;Ou是n1×(M-n1)的零矩阵;U1是由U的前n1列组成的矩阵,U2是由U的其余列组成的矩阵,V1是由V的前n1列组成的矩阵,V2是由V的其余列组成的矩阵;
在设定Y*=U1E1V1 T后,根据X*=(Y*)T获取n1模态的数据矩阵X*
将X*的每一行的行向量复原分别为对第sj个矢状面在一个预设周期内每一时刻预处理后的医疗图像;
将所有预处理后的医疗图像的第ti个图像确定为第ti时刻第sj个矢状面预处理后的医疗图像I2
5.根据权利要求1所述的方法,其特征在于,所述对所述任意两个方向流场图进行图像腐蚀,得到第一流场腐蚀图以及第二流场腐蚀图,包括:
基于预设个数nt阈值的大津法分别获取第一流场图像对应的第一流场(nt+1)值图以及第二流场图像对应的第二流场(nt+1)值图;其中,所述第一流场图像为所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中的任一方向流场图像;所述第二流场图像为除第一流场图像以外的其他任意方向流场图像;
分别对第一流场(nt+1)值图以及第二流场(nt+1)值图进行图像腐蚀,得到第一流场腐蚀图以及第二流场腐蚀图。
6.根据权利要求1所述的方法,其特征在于,所述将所述目标对象边界的边界点按照预设的修正策略进行修正,包括:
按照预设的排列顺序对所有边界点进行排列及编号,并获取各边界点的曲率;
将所有边界点中曲率最大的边界点进行标记后执行以下步骤:
步骤a,以任一曲率最大的边界点Q0为中心,在预设的范围内搜索排列在Q0之后且与Q0的编号差别小于预设常数的所有最大曲率值点Q=(Q1,Q2,...,Qs);
步骤b,以Q0为第一个点,以Q中每一个点Qκ∈Q为第二个点,Q0Qκ为一边按照预设的方向各组成方框;
步骤c,获取所述方框的分形维数;
步骤d,当基于所述分形维数确定Q0Qκ两点之间的边界异常时,在Q0Qκ两点之间生成拟合曲线作为修正后的边界;
针对每个曲率最大的边界点执行步骤a至步骤d,获得修正后的边界点。
7.根据权利要求6所述的方法,其特征在于,所述获取所述方框的分形维数,包括:
获取所述方框的计盒维数,并将所述方框的计盒维数作为所述方框的分形维数。
8.根据权利要求1所述的方法,其特征在于,所述获取每个边界点的亚像素位置,得到所述目标对象的亚像素边界,包括:
在对所述目标对象边界的边界点修正完毕后,针对每个边界点执行以下步骤:
设定在以边界点(x0,y0)为原点的坐标轴下,直线y=a+bx为真实边界线;
以所述边界点(x0,y0)为中心确定大小为5×3个方格的周边范围;
通过所述直线y=a+bx将所述周边范围划分为两个区域;其中,所述直线y=a+bx上方的第一区域的平均灰度值为A1,所述直线y=a+bx下方的第二区域的平均灰度值为A2
设定所述周边范围内每个点(x',y')满足
Figure FDA0002651746530000061
其中,h是像素边长,Px',y'是(x',y')位置的像素在所述直线y=a+bx下方所占的面积;
按照式7分别获取所述周边范围的左、中、右列在所述直线y=a+bx下方的面积L,M,R:
Figure FDA0002651746530000062
按照式8分别获取所述周边范围的左、中、右列的灰度值之和SL,SM,SR
Figure FDA0002651746530000063
按照式9对所述第一区域的平均灰度值A1和所述第二区域的平均灰度值A2进行拟合:
Figure FDA0002651746530000064
基于式7、式8和式9获得所述直线y=a+bx中的参数
Figure FDA0002651746530000065
以及
Figure FDA0002651746530000066
在以边界点(x0,y0)为原点的坐标轴下,将所述原点的亚像素位置确认为(0,a);
根据所述原点的亚像素位置(0,a)确定所述边界点(x0,y0)的亚像素坐标在预处理后的医疗图像中为(x0,y0+a)。
9.根据权利要求1所述的方法,其特征在于,所述方法还包括:
在获取到第ti时刻每个矢状面的原始图像组对应目标对象的亚像素边界后,基于所述目标对象的所有亚像素边界进行三维组装,获得所述目标对象的原始三维体;
将所述原始三维体进行平滑后,获得所述目标对象的三维重构体。
10.一种医疗图像处理装置,其特征在于,所述装置包括:预处理部分、分割部分、修正部分和亚像素拟合部分;其中,
所述预处理部分,配置为针对第ti时刻第sj个矢状面的原始图像组中的原始医疗图像按照预设的图像预处理策略去除灰度偏差场及噪声,获得预处理后的医疗图像;其中,所述原始图像组包括原始医疗图像、x方向流场图像、y方向流场图像和z方向流场图像;
所述分割部分,配置为基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出目标对象边界;
其中,所述基于所述x方向流场图像、所述y方向流场图像和所述z方向流场图像中任意两个方向流场图像从所述预处理后的医疗图像中获取目标对象轮廓,并基于所述目标对象轮廓从所述预处理后的医疗图像中分割出目标对象边界,包括:对所述任意两个方向流场图进行图像腐蚀,得到第一流场腐蚀图以及第二流场腐蚀图,且将所述第一流场腐蚀图以及所述第二流场腐蚀图中灰度值取最大值的区域确定为第一流场腐蚀图的强信号区域以及第二流场腐蚀图的强信号区域,将所述第一流场腐蚀图的强信号区域以及所述第二流场腐蚀图的强信号区域之间的连通区域进行合并后,对合并后的连通区域进行光滑运算,获得目标对象区域二值图;基于所述目标对象区域二值图的骨架后,以目标对象区域中心点为顶点,确定所述目标对象区域二值图骨架上的种子点;基于所述种子点对所述预处理后的医疗图像进行轮廓分割,获得每个种子点对应的最终轮廓;
所述修正部分,配置为将所述目标对象边界的边界点按照预设的修正策略进行修正;
所述亚像素拟合部分,配置为所述修正部分将所述目标对象边界的边界点按照预设的修正策略进行修正后,获取每个边界点的亚像素位置,得到目标对象的亚像素边界。
11.根据权利要求10所述的装置,其特征在于,所述装置还包括:三维重构部分,配置为:
在获取到第ti时刻每个矢状面的原始图像组对应目标对象的亚像素边界后,基于所述目标对象的所有亚像素边界进行三维组装,获得所述目标对象的原始三维体;以及,
将所述原始三维体进行平滑后,获得所述目标对象的三维重构体。
12.一种医疗图像处理装置,其特征在于,所述装置包括:存储器、处理器及存储在所述存储器上并能够在所述处理器上运行的计算机程序,所述计算机程序被所述处理器执行时实现如权利要求1至9任一项所述医疗图像处理方法的步骤。
13.一种计算机存储介质,其特征在于,所述计算机存储介质存储有医疗图像处理程序,当所述医疗图像处理程序被至少一个处理器执行时实现如权利要求1至9任一项所述的医疗图像处理方法的步骤。
CN201810534684.1A 2018-05-29 2018-05-29 一种医疗图像的处理方法、装置及计算机存储介质 Active CN108898578B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810534684.1A CN108898578B (zh) 2018-05-29 2018-05-29 一种医疗图像的处理方法、装置及计算机存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810534684.1A CN108898578B (zh) 2018-05-29 2018-05-29 一种医疗图像的处理方法、装置及计算机存储介质

Publications (2)

Publication Number Publication Date
CN108898578A CN108898578A (zh) 2018-11-27
CN108898578B true CN108898578B (zh) 2020-12-08

Family

ID=64343540

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810534684.1A Active CN108898578B (zh) 2018-05-29 2018-05-29 一种医疗图像的处理方法、装置及计算机存储介质

Country Status (1)

Country Link
CN (1) CN108898578B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110046646B (zh) * 2019-03-07 2023-06-30 深圳先进技术研究院 图像处理方法、系统、计算设备及存储介质
CN110008895B (zh) * 2019-04-01 2023-01-17 中南林业科技大学 一种赛道特征识别方法及智能赛车
CN110047078B (zh) * 2019-04-18 2021-11-09 北京市商汤科技开发有限公司 图像处理方法及装置、电子设备和存储介质
CN112734773B (zh) * 2021-01-28 2023-03-21 依未科技(北京)有限公司 一种亚像素级眼底血管分割方法、装置、介质和设备

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101706843A (zh) * 2009-11-16 2010-05-12 杭州电子科技大学 一种乳腺cr图像交互式读片方法
CN107154036A (zh) * 2017-03-24 2017-09-12 中南大学 基于亚像素的血管分割方法及其分割系统
CN107680107A (zh) * 2017-10-30 2018-02-09 西北工业大学 一种基于多图谱的扩散张量磁共振图像的自动分割方法
CN108038862A (zh) * 2017-12-11 2018-05-15 深圳市图智能科技有限公司 一种交互式医学图像智能分割建模方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7782056B2 (en) * 2007-12-13 2010-08-24 Isis Innovation Ltd. Systems and methods for correction of inhomogeneities in magnetic resonance images
US9558561B2 (en) * 2015-01-06 2017-01-31 Varian Medical Systems International Ag Semiautomatic drawing tool for image segmentation

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101706843A (zh) * 2009-11-16 2010-05-12 杭州电子科技大学 一种乳腺cr图像交互式读片方法
CN107154036A (zh) * 2017-03-24 2017-09-12 中南大学 基于亚像素的血管分割方法及其分割系统
CN107680107A (zh) * 2017-10-30 2018-02-09 西北工业大学 一种基于多图谱的扩散张量磁共振图像的自动分割方法
CN108038862A (zh) * 2017-12-11 2018-05-15 深圳市图智能科技有限公司 一种交互式医学图像智能分割建模方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
A Watershed Approach for Improving Medical Image Segmentation;E A Zanaty等;《Comput Methods Biomech Biomed Engin》;20120330;第16卷(第12期);第1262-1272页 *
MR图像偏差场校正与分割的研究;yjdoc669上传;《道客巴巴》;20140810;第9-120页 *
SegNet: A Deep Convolutional Encoder-Decoder Architecture for Image Segmentation;Vijay Badrinarayanan等;《IEEE Transactions on Pattern Analysis & Machine Intelligence》;20171231;第39卷(第12期);第2484-2486页 *
图像分形维数计算及其边缘提取;李鹏飞等;《吉林大学学报(信息科学版)》;20110331;第29卷(第2期);第152-157页 *
基于局部偏差的图像分割与偏置场校正模型;黄臣程等;《四川理工学院学报(自然科学版)》;20171031;第30卷(第5期);摘要、第45页 *

Also Published As

Publication number Publication date
CN108898578A (zh) 2018-11-27

Similar Documents

Publication Publication Date Title
CN108898578B (zh) 一种医疗图像的处理方法、装置及计算机存储介质
US20190172205A1 (en) Method and system for extracting vasculature
Papież et al. An implicit sliding-motion preserving regularisation via bilateral filtering for deformable image registration
CN112508965B (zh) 医学影像中正常器官的轮廓线自动勾画系统
CN107133946B (zh) 医学图像处理方法、装置及设备
JP4885138B2 (ja) 一連の画像における動き修正のための方法およびシステム
Martínez et al. Segmentation of pelvic structures for planning CT using a geometrical shape model tuned by a multi-scale edge detector
CN110751651B (zh) 基于多尺度迁移学习的mri胰腺图像分割方法
Park et al. Deformable registration of CT and cone-beam CT with local intensity matching
Saadatmand-Tarzjan Self-affine snake for medical image segmentation
CN108038840B (zh) 一种图像处理方法、装置、图像处理设备及存储介质
CN115830016B (zh) 医学图像配准模型训练方法及设备
Yu et al. Accelerated gradient-based free form deformable registration for online adaptive radiotherapy
Matsopoulos et al. Thoracic non-rigid registration combining self-organizing maps and radial basis functions
CN113421226B (zh) 基于互信息的ct-dr多模态食管图像配准方法及系统
Perez–Gonzalez et al. Probabilistic learning coherent point drift for 3D ultrasound fetal head registration
CN116703994B (zh) 医学图像配准的方法、计算设备和计算机可读存储介质
WO2023232067A1 (en) Systems and methods for lesion region identification
Shi et al. Fast shading correction for cone-beam CT via partitioned tissue classification
CN112102327B (zh) 一种图像处理方法、装置及计算机可读存储介质
CN113822796A (zh) 基于改进的图像金字塔的多模态脑部影像配准方法
CN109671131B (zh) 影像校正方法、装置、医疗影像设备及存储介质
Al Abboodi et al. Supervised Transfer Learning for Multi Organs 3D Segmentation With Registration Tools for Metal Artifact Reduction in CT Images
Ninomiya et al. Feasibility of anatomical feature points for the estimation of prostate locations in the Bayesian delineation frameworks for prostate cancer radiotherapy
CN111612867A (zh) 运动伪影校正方法、装置、计算机设备和可读存储介质

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
CB02 Change of applicant information
CB02 Change of applicant information

Address after: Room 707, building 2, Xizi international jinzuo, Nanyuan street, Yuhang District, Hangzhou City, Zhejiang Province

Applicant after: HANGZHOU SHENGSHI TECHNOLOGY Co.,Ltd.

Address before: 311100 Hangzhou, Yuhang, Zhejiang Linping new town Nanyuan street, new far CBC1 block 1502 room.

Applicant before: HANGZHOU SHENGSHI TECHNOLOGY Co.,Ltd.

GR01 Patent grant
GR01 Patent grant
PE01 Entry into force of the registration of the contract for pledge of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: The invention relates to a medical image processing method, a device and a computer storage medium

Effective date of registration: 20210602

Granted publication date: 20201208

Pledgee: Hangzhou United Rural Commercial Bank Limited by Share Ltd. Baoshan Branch

Pledgor: HANGZHOU SHENGSHI TECHNOLOGY Co.,Ltd.

Registration number: Y2021980004299

PC01 Cancellation of the registration of the contract for pledge of patent right
PC01 Cancellation of the registration of the contract for pledge of patent right

Date of cancellation: 20231113

Granted publication date: 20201208

Pledgee: Hangzhou United Rural Commercial Bank Limited by Share Ltd. Baoshan Branch

Pledgor: HANGZHOU SHENGSHI TECHNOLOGY Co.,Ltd.

Registration number: Y2021980004299