本发明涉及医学图像处理技术领域，具体涉及一种基于cta图像的动态心脏各腔室分割方法。
背景技术：
心脏是人体心血管系统的关键器官，内部组织结构较复杂，活体处于不断的舒缩中，而目前心脏的影像资料多为正交的二维平面影像像，每个平面都是在某个时刻的心脏运动过程中的某个采样，临床医师较难通过二维图像去想象病变部位的三维形态和空间毗邻关系，这就决定了心脏疾病在影像学上定位和定性诊断较困难。中国可视化人(cvh)项目，于2002年10月，由第三军医大学发起，该项目是从尸体获得包含丰富解剖学信息的高分辨率、连续断层切片图像的数据集。此数据集的优点是，它不仅可以清楚地观察到的心脏纤维骨架(该数据集中连续断层图像中各个组织呈自然的纹理颜色，容易区分不同组织的边界)，而且保持了解剖结构的原始空间关系。cvh数据集层间距小于临床连续断层影像。水平分辨率越大，层间距越小，能够呈现的解剖细节就越丰富，特别是二尖瓣、三尖瓣等结构需要在心脏轴面来观察的结构。通过数字化人体数据集建立高精度的心脏可视化解剖模型可以反应心脏毗邻解剖结构原位的空间关系[1]，但是如何将这些信息用于活体的断层解剖图像，以指导临床的诊断和资料，这是心脏可视化面临的一个挑战。
借助医学影像技术，定性和定量分析心脏的解剖及运动规律，进而研究其与心脏疾病之间的关系是现代医学研究的热点。利用计算机模型对形态进行仿真，模拟真实心脏的运动过程及病理状态，进一步认识心血管系统的运动规律和本质，深入了解心血管系统疾病和两次表现之间的关系。关键问题是形变模型要适应心脏组织的周期性变化和不同的心脏个体，即以模型为引导的分割，用模型去适应图像特征的变化。心脏cta(cardiaccomputedtomographyarteriography，计算机断层造影术)以成像速度快、分辨率高、信息量大，作为心脏检查的重要手段。基于cta图像的量化评估对心功能的评估上具有重要的作用，特别是在动态定量评估方面能通过较高的时间分辨率捕捉心动周期内心脏内部形态变化特征，进而评估心脏生理功能。因此，基于cta的心脏量化评估引起国内外学者的高度关注。定量评估的核心是对心脏区域的分割。心脏分割的目的主要是从复杂的心脏图像中提取出心脏感兴趣区域的轮廓曲线，然而由于医学图像的复杂性与多样性，成像质量常受到噪声、运动伪影、局部体效应的干扰，在很多场合下还是会选择人工手动分割的方式进行处理[2]。目前国内外学者提出了许多针对心脏中某些解剖结构的分割算法，其主要关注点在于心室的分割[3]。
由于心脏形态结构的复杂性，心脏分割的准确性和实时性一直是研究的热点。在cta图像中,单一心腔的分割相对比较容易特别是左心室，但是对于四个腔室的同时分割比较困难。其原因是心腔在解剖结构上具有特殊性，使医学图像呈现出一些不容易分割的状态：1)单纯依靠图像特征难以完全实现解剖结构的分割。心腔的cta图像中的灰度变化与腔体的分隔线并非完全对应，某些在灰度上连通的区域在心脏解剖结构上是分离的，而心肌与周围组织的灰度近似，仅依靠灰度难以区分；2)心脏一直处于运动状态，成像伪影及噪声会影响成像效果；3)由于心脏结构的特殊性，心脏图像中会出现曲率很大的凹陷区域，这部分图像难以分割[4]。
因此，如何借助解剖知识实现心脏的分割是进行心脏各腔室量化评估的关键问题。心脏纤维骨架在心脏运动中起支点和稳定的作用，心脏纤维骨架既是几何形态的中心，又是位于各方向运动的中心，其“双中心”的地位对理解心脏运动规律有极其重要的意义。超声能实时捕捉到瓣膜的运动但是对周围解剖结构的显示却不理想，ct、mri能显示更详细的解剖结构，但是重建出瓣膜的三维形态非常困难，关键在于其形态走行复杂从无法准确判断二维的影像上起止点，特别是在心肌附着处难以区分其边界。因此，建立人体心脏原位的心脏纤维骨架的三维模型是非常困难的。人们尝试着使用动物和使用不同的技术，包括解剖或影像学检查人的研究，已有学者研究了二尖瓣环和三尖瓣环的形态学特征及动态特征[5]，然而已有的研究中少有对同一心脏四个瓣膜的情况进行分析的。对当前的成像后处理方式而言，同时定位和识别4个瓣膜仍然面临困难。数字化人体是高精度的彩色纹理连续断层图像，对心腔内部精细结构的鉴别和区分具有较大优势，在此基础上建立的数字化模型维持了心脏内部毗邻解剖结构的原位空间关系，组织间没有相对的位移和形变，便于和临床影像有机结合。可以将已建立的详细的解剖结构模型结合cta、mri等多种模式下的可视化方式，从而全面了解心脏在生理、病理下的状态。
1.郭燕丽,人体心脏薄层断面解剖和三维可视化与多平面tee的对照研究.2003,第三军医大学.
2.付增良,陈晓军,叶铭,李峰,王成焘,心脏ct图像分割方法.计算机工程,2009.35(12):p.189-191.
3.zhou,y.,w.r.shi,w.chen,y.l.chen,y.li,etal.,activecontoursdrivenbylocalizingregionandedge-basedintensityfittingenergywithapplicationtosegmentationoftheleftventricleincardiacctimages.neurocomputing,2015.156(c):p.199-210.
4.贺娟,基于ct图像的全心脏分割算法研究.2015,浙江大学.
5.maffessanti,f.,p.gripari,g.pontone,d.andreini,e.bertella,etal.,three-dimensionaldynamicassessmentoftricuspidandmitralannuliusingcardiovascularmagneticresonance.europeanheartjournalcardiovascularimaging,2013.14(10):p.986-995.
技术实现要素：
鉴于以上技术问题，本发明提出了一种基于cta图像的动态心脏各腔室分割方法，该发明通过结合解剖学与生理学知识，用cvh高分辨率切片重建的心脏纤维骨架三维模型确定其在心脏三维空间中的位置，实现连续影像数据心脏平面的快速构建，以解决对心脏四个腔室同时进行分割及不同个体cta图像在可视化量化分析的一致性问题。
本发明是通过如下技术方案来实现的：
一种基于cta图像的动态心脏各腔室分割方法，包括以下步骤：
s1，获取心脏cta数据和cvh数据集；
s2，针对获得的心脏cta图像，获得ct值的灰度直方图，针对该直方图设定一个阈值，利用阈值法获取整个心脏血液分布的分割数据，实现心脏血液分布的最大连通；
s3，截断主动脉根部与左心室的连通关系，实现冠状动脉的分割；
s4，利用cvh数据集获得的心脏纤维骨架三维模型实现对心脏各个腔室的空间进行划分；
s5，将心脏纤维骨架与心脏cta数据的主动脉根部之间进行位置、大小、角度等的变换，实现两个图像的配准，获得cvh空间到心脏cta图像的几何变换矩阵；
s6，在cvh数据集上定义好四个瓣环的空间坐标，然后将上述步骤s5获得的几何变换矩阵应用到cta数据集上，获得定义好的各个腔室的空间坐标，并通过三维数据集上的重采样重建出具有解剖细节的cta数据集的四个瓣环平面图像，实现心脏各个腔室的体数据分割；
s7，对心脏cta数据至少两个时相断层影像的每个时相重复步骤s1到步骤s6，实现动态心脏的影像分割。
作为优选，所述步骤s3具体为，截断主动脉根部和左心室的连接部分，对主动脉根部使用连通域的方法，分割出主动脉根部图像，获得冠状动脉的分割图像。
作为优选，所述步骤s4具体为，针对cvh数据集，人工进行心脏纤维骨架解剖结构的分割，获得心脏纤维骨架的三维模型，实现对心脏各个腔室空间的划分。
作为优选，所述步骤s5具体包括以下步骤：
s51,由于心脏cta成像的结果与心脏纤维骨架数据重建模型呈三维的镜像关系，所以针对心脏cta数据，首先需要通过镜像变换调整主动脉根部图像的空间位置，使之与心脏纤维骨架的位置相对应；
s52,以主动脉瓣环三皇冠三个顶点为准拟合得到的圆环(环1)与心脏纤维骨架环主动脉瓣环三皇冠三个底点拟合圆环(环2)之间的半径比作为图像尺度变换的参数，通过尺度变换实现对不同病人心脏体积与纤维骨架环大小的匹配；
s53,以环1所在平面与环2所在平面之间的夹角作为旋转变换的参数，以此来消除由于个体差异和心脏运动状态的不同导致的角度变化；
s54,通过平移变换实现对整个冠状动脉和纤维骨架空间位置的对应，完成图像的配准。配准后，通过以下公式获得主动脉根部到四个瓣环平面的几何变换矩阵g，f′(x′,y′,z′)＝f(x,y,z)*g
其中，f(x,y,z)和f(x′,y′,z′)分别是cvh和cta主动脉瓣中点集所组成的向量，x,y,z和x′,y′,z′分别代表变换前和变换后的x，y，z轴的坐标。
g为几何变换矩阵，g(s,t,r)，其中s为尺度变换矩阵，t是平移变换矩阵，r为旋转矩阵，
其中，dx,dy,dz为心脏纤维骨架分别在x，y，z轴的变化量，rx(θ),ry(θ),rz(θ)分别代表在x，y，z轴上旋转的角度，
由以下公式得到几何变换之后的坐标：
f′(x′,y′,z′)＝f(x,y,z)*g，f(x,y,z)∈f(x,y,z)
其中，f(x,y,z)是由cvh模型的三维坐标和四个瓣环平面组成，f′(x′,y′,z′)是将cvh定义的坐标关系f(x,y,z)转换到cta图像空间坐标的集合。
作为优选，所述步骤s6具体包括以下步骤：
s61,计算四个瓣环平面与原始cta数据图像空间的交点；
s62,获得四个瓣环平面x轴和y轴的取值位置；
s63,计算四个瓣环平面的法向量n，其中，n＝(a,b,c)，a、
s64,m和m′为平面上任意两点，mm′＝(x-x0,y-y0,z-z0)；根据空间平面的点法式方程及以下公式获得平面上每个像素点z取值：
s65,构建四个瓣环平面图像矩阵，读取原始图像空间中的x，y，z坐标对应的像素值来填充图像矩阵对应位置的值，其余部分用非数值nan格式来填充。
作为优选，所述步骤s7中，将心脏cta数据划分成20个时相断层影像，对这20个时相断层影像的每个时相重复步骤s1到步骤s6，实现动态心脏影像的分割。
作为优选，所述四个瓣环包括二尖瓣环、三尖瓣环、主动脉瓣环和肺动脉瓣环。
本发明的有益效果是：
1、本发明能够同时对心脏的四个腔室的血液(即左心室、右心室、左心房、右心房)进行分割，算法简单、速度快；
2、对心脏分割的难点即右心室、右心房具有较好的分割效果；
3、量化计算结果符合临床研究，程序鲁棒性较好；
4、解决了不同个体cta图像在可视化量化分析的一致性问题。
具体实施方式
下面将结合具体实施例详细地描述：
本发明首先以瓣环的空间关系来设计瓣环平面，获得主动脉根部到四个瓣环平面的几何变换关系；然后构造出心脏cta图像与cvh心脏纤维骨架模型之间的空间关系，对心脏cta图像与心脏纤维骨架模型进行配准；再对变换后的模型应用相对的几何变换关系以构建个体的心脏瓣环平面，对全心脏cta图像进行分割；对分割结果通过像素体积量化累加得到各腔室体积。
一种基于cta图像的动态心脏各腔室分割方法，包括以下步骤：
s1，获取心脏cta数据和cvh数据集；
cta数据的获取：ct血管造影数据由西南医院放射科提供，采集参数是探测器准直0.6mm，旋转时间330ms，管电压120kv，管电流400ma，采用回顾性心电图门控方式在整个ecg周期获得同步化连续扫描图像形成一个对比增强的体数据集，可以在心动周期的任何时相进行图像重建。图像重建层厚0.75mm，重建间距0.4mm，重建卷积核b26f。在心动周期的重建窗口的位置可以根据个体进行选择，以尽量减少误差。单幅图像的分辨率为512*512像素，像素分辨率为0.12mm*0.12mm。原始序列断层图像用于在特定的解剖标志产生心脏平面。
cvh的获取方法：将人体标本,经外形测量、血管灌注后,用5％明胶包埋,置入-30℃盐水池中冰冻1周,然后在-25℃低温实验室中用tk-6350型数控铣床(铣切精度为0.001mm)从头至足逐层铣切.逐层用高清晰度数码相机获取连续断层图像数据，得到cvh数据集。
在我们的研究中，心脏建模所需的原始图像来自于cvh5(第5例中国数字化人体)数据集。该数据集为真实人体组织的彩色图像。从中选取胸部共计780张连续横断面图像作为待分割图像，每个层面均能连续完整地显示心内细微结构。每张断层图像的分辨率为3072*2048。其像素分辨率为0.12mm*0.12mm，层间距为0.2mm。
s2，针对获得的心脏cta图像，获得ct值的灰度直方图，针对该直方图设定一个阈值，利用阈值法获取整个心脏血液分布的分割数据，实现心脏血液分布的最大连通；
灰度直方图是关于灰度级分布的函数，是对图像中灰度级分布的统计。灰度直方图是将数字图像中的所有像素，按照灰度值的大小，统计其容积。
根据cta图像的灰度直方图，设定一个阈值(根据不同个体,一般在160-220之间)，在本研究中取的阈值为180，就可以得到心脏血液分布的分割数据。
s3，截断主动脉根部与左心室的连通关系，实现冠状动脉的分割；
所述步骤s3具体为，截断主动脉根部和左心室的连接部分，对主动脉根部使用连通域的方法，分割出主动脉根部图像，获得冠状动脉的分割图像。
s4，利用cvh数据集获得的心脏纤维骨架三维模型对心脏各个腔室的空间进行划分；
所述步骤s4具体为，针对cvh数据集，人工进行心脏纤维骨架解剖结构的分割，获得心脏纤维骨架的三维模型，实现对心脏各个腔室空间的划分。
cvh数据集为彩色图像，各个组织结构均以其自然的纹理颜色呈现，更容易区分组织结构之间的边界。彩色图像一般采取r、g、b三个通道组合来表示一幅图像，为了减少对内存的消耗，只选择一个通道的图像导入amira软件，将连续断层堆积起来形成体数据，分割前对每个待分割的结构建立一个标签，其中包括颜色和名称的对应关系，然后选择对应的标签人工进行解剖结构的分割。最后将分割好的形态添加到该标签的数据结构中，作为待重建数据。在amira软件中选择compute->surfacegen，重建平滑参数设置为不平滑，其余保持默认参数设置。即生成已分割结构的三维面绘制模型。
所述人工分割方法为在amira软件中，手动勾画出心脏纤维骨架的解剖结构。具体步骤：点击opendata，将780张连续断层图像全部导入，通道转换选择channel1红色通道，转换为软件可以处理的灰度图像。voxelsize设置为0.12×0.12×0.2，以保证分割重建模型长宽高比例的正确。右键单击数据集，选择labelling-labelfield，得到am格式分割文件，选中后单击segmentationeditor，即可进入分割模式。在materials选项中新建分割区，并设置其颜色。
s5，将心脏纤维骨架与心脏cta数据的主动脉根部之间进行平移、缩放、旋转角度等的变换，实现两个图像的配准，获得cvh空间到心脏cta图像的几何变换矩阵；
所述步骤s5具体包括以下步骤，
s51,由于心脏cta成像的结果与心脏纤维骨架数据重建模型呈三维的镜像关系，所以针对心脏cta数据，首先需要通过镜像变换调整主动脉根部图像的空间位置，使之与心脏纤维骨架的位置相对应；
s52,以主动脉瓣环三皇冠三个顶点为准拟合得到的圆环(环1)与心脏纤维骨架环主动脉瓣环三皇冠三个底点拟合圆环(环2)之间的半径比作为图像尺度变换的参数，通过尺度变换实现对不同病人心脏体积与纤维骨架环大小的匹配；
s53,以环1所在平面与环2所在平面之间的夹角作为旋转变换的参数，以此来消除由于个体差异和心脏运动状态的不同导致的角度变化；
s54,通过平移变换实现对整个冠状动脉和纤维骨架空间位置的对应，完成图像的配准。配准后，通过以下公式获得主动脉根部到四个瓣环平面的几何变换矩阵g，
f′(x′,y′,z′)＝f(x,y,z)*g
其中，f(x,y,z)和f(x′,y′,z′)分别是cvh和cta主动脉瓣中点集所组成的向量，x,y,z和x′,y′,z′分别代表变换前和变换后的x，y，z轴的坐标。
g为几何变换矩阵，g(s,t,r)，其中s为尺度变换矩阵，t是平移变换矩阵，r为旋转矩阵。
其中，dx,dy,dz为心脏纤维骨架分别在x，y，z轴的变化量，rx(θ),ry(θ),rz(θ)分别代表在x，y，z轴上旋转的角度。
由以下公式得到几何变换之后的坐标：
f′(x′,y′,z′)＝f(x,y,z)*g，f(x,y,z)∈f(x,y,z)
其中，f(x,y,z)是由cvh模型的三维坐标和四个瓣环平面组成，f′(x′,y′,z′)是将cvh定义的坐标关系f(x,y,z)转换到cta图像空间坐标的集合。
s6，在cvh数据集上定义好四个瓣环的空间坐标，然后将上述步骤s5获得的几何变换矩阵应用到cta数据集上，获得定义好的各个腔室的空间坐标，并通过三维数据集上的重采样重建出具有解剖细节的cta数据集的四个瓣环平面图像，实现心脏各个腔室的体数据分割。所述步骤s6具体包括以下步骤，
s61,计算四个瓣环平面与原始cta数据图像空间的交点；
s62,获得四个瓣环平面x轴和y轴的取值位置；
s63,计算四个瓣环平面的法向量n，其中，n＝(a,b,c)，a、
s64,m和m′为平面上任意两点，mm′＝(x-x0,y-y0,z-z0)；根据空间平面的点法式方程及以下公式获得平面上每个像素点z取值：
s65,构建四个瓣环平面图像矩阵，读取原始图像空间中的x，y，z坐标对应的像素值来填充图像矩阵对应位置的值，其余部分用非数值nan格式来填充。
s7，对心脏cta数据至少两个时相断层影像的每个时相重复步骤s1-s6，实现动态心脏的影像分割。
所述步骤s7中：将心脏cta数据划分成20个时相断层影像，对这20个时相断层影像的每个时相重复步骤s1到步骤s6，实现动态心脏影像的分割。
所述四个瓣环包括二尖瓣环、三尖瓣环、主动脉瓣环和肺动脉瓣环。
采用上述方法实现了以下有益效果：
1、本发明能够同时对心脏的四个腔室的血液(即左心室、右心室、左心房、右心房)进行分割，算法简单、速度快；
cta图像数据像素数均为512×512×293，通过本方式处理过程耗时最长为182s，最短为167s。
2、对心脏分割的难点即右心室、右心房具有较好的分割效果；
心房的划分相对于心室而言，其难点在于流入道的分隔。左心房流入道包括四条肺静脉，右心房流入道有上下腔静脉，周边的血管网错综复杂。因此，在分割过程中，我们依据心脏各腔室连通通路与心脏纤维骨架相应瓣环的关系，通过拟合平面将各腔室连接的通路截断，有效实现了心脏各个腔室的分割；
3、量化计算结果符合临床研究，程序鲁棒性较好；
本研究分析与处理不同时期不同状态的心脏cta图像共21例，其中19例图像心脏四腔室完全分割，剩余两例由于cta图像心房灌注部分连通，无法进行有效的截断，只能得到其心室造影的分割图像。排除成像质量的因素以外，所处理的数据中尚未出现过大的分割误差。
4、解决了不同个体cta图像在可视化量化分析的一致性问题。由于心脏内部形态不规则，各个腔室的起止位置的确定是分割中的难点，且不同患者的心脏形态有差异，靠人为界定临界位置的起止点容易引入量化误差。本发明利用人体解剖结构共性特征去匹配患者的个性化影像数据，按统一准则进行腔室的分割，有效解决量化分割的一致性问题。
应当理解的是，本方法在具体应用中并不局限于cta影像，凡涉及三维或四维心脏各个腔室及腔内血液分割的应用，本发明均适用。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作任何各种改动和修改，这些等价形式同样落于本申请所附权利要求书所限制。