基于多层螺旋计算机层析数据的冠脉三维运动估计

基于多层螺旋计算机层析数据的冠脉三维运动估计

1血管模型的建立和动力学冠状动脉粥样硬化的估计是心脏疾病研究领域的一个热点问题之一。x线成像图像用于跟踪血管，并分析运动的变化。由于X射线造影成像的本质是透视投影,不同深度的空间血管段在投影图像上具有不同的放大率,因此仅从二维图像估计运动是不准确的。孙正等人提出将三维重建和运动估计结合起来进行,由两幅不同角度的单面二维造影图像重建出三维血管,之后在三维空间内对血管运动进行分析,但算法的结果在很大程度上取决于三维重建的精度,而在二维图像的获取过程中,平台的移动不可避免,需要通过优化变换矩阵减小重建误差。近年,多层螺旋计算机层析(MSCT)成像技术发展迅速,特别在冠心病诊治中,MSCT冠脉造影(CTA)作为一种无创、三维成像的检查方式,在临床上应用日益广泛。不断提高的时空分辨率,使得利用CTA数据分析运动的血管成为可能。目前,国内外研究多集中在利用CTA数据构建三维静态血管模型进行狭窄段分析。本文着眼于动态模型的建立及运动场的估计。对4D-CTA数据进行分析,利用局部血管增强及自适应区域生长方法提取出各个期相的冠脉树并进行细化,对骨架点进行基于连贯点漂移(CPD)算法的配准,从而估计运动场。通过对模拟数据和真实全期相CTA数据的处理,验证本文提出方法的精度,并定量描述主要分支的位移幅度与速率。2血管预分割和溶栓设置本文在前期工作的基础上,应用基于Hessian矩阵的局部血管增强与自适应阈值区域生长相结合的方法进行血管分割。由于造影剂的作用,CTA数据中的血管的灰度通常大于周围组织,表现为较为明亮的细长管状体。对体素数据进行差分计算,获得其二阶导数矩阵,即Hessian矩阵。对于血管区域内的点,其Hessian矩阵特征值特点为λ1≈0,λ2<0,λ3<0且|λ2|≈|λ3|≫|λ1|,其中λ1对应的特征方向为血管方向,通过分析矩阵的特征值,区分血管和非血管区域内的体素,对数据进行滤波,从而达到增强血管、抑制其他组织的目的。为更好地适应血管灰度的不均匀性,将数据空间划分为一系列的感兴趣区域(ROI),在每一ROI中先构造血管增强滤波函数进行预处理,再进行血管分割。首先指定初始种子点,选择一个特定大小的ROI,计算该区域内Hessian矩阵的特征值,利用特征值构造血管增强滤波函数,抑制非血管结构。对增强后的数据进行预分割,从初始值开始,逐次降低阈值,统计分割出的区域体积,当阈值低到一定程度时,分割出的结果会泄漏到血管周围区域,引起一个比较大的体积突变,则上一次的阈值为最佳阈值。将最佳阈值作为最终阈值进行区域生长,利用质心法计算分割出的血管在ROI立方体各个截面上的质心,作为新的种子点,生成新的ROI,重复分割过程。每一次分割的局部结果及时进行全局融合,最终得到完整的冠状动脉树。ROI的大小由多次实验确定,经验值为边长为31voxel的立方体。分割后的冠脉采用基于数学形态学的三维细化算法,从三维空间6个方向上进行迭代判断,由外向内删除“简化点”,保留端点、分支点等特征点,直至得到26连通的单体素曲线,经过剪枝等后处理,消除细小毛刺,获得保持原血管树拓扑结构及连通特性的骨架线。对各个期相的数据均进行同样的冠脉分割及细化处理步骤,得到同一冠脉分支不同时刻的骨架点序列,记为Xt={xt,i}(i=1,2,……nt),t为表征不同时刻的序号,nt为该时刻分支上点的个数。3高斯模型权重的确定对于运动场的估计,一般转化为点集之间的配准问题进行求解。假设点集Xt与Xt+1分别为相邻两个时刻某分支上的骨架点,为方便,Xt记为Y,Xt+1记为X。对两点集配准,就是确定两个点集之间的对应点对{xi,yj}及空间位置变换关系T,使得xi=T(yj)。图像配准算法很多,具体到冠脉运动,因为轨迹复杂,且在运动过程中发生形变,所以T为非刚体变换。较为常见的非刚体配准方法有自由形变模型(FFD)、鲁棒点配准(RPM)和CPD算法等,其中CPD算法将运动连贯理论(MCT)应用于配准过程,保证了相邻点之间变换的平滑性,在非刚体点配准中取得了较为满意的效果。假设X与Y分别为数据点集与模型点集,记为XN×D=(x1,x2,…xN),YM×D=(y1,y2,…yM),维度D=3。记Y变换后的点集为T(Y,θ),θ为变换参数。若将X与T(Y,θ)视为高斯混合模型的点分布与相应的模型中心,则X的分布服从p(xn)=Μ+1∑i=1Ρ(m)p(xn|m)=w1Ν+(1-w)Μ∑i=1Ρ(m)1(2πσ2)D/2exp(-∥xn-Τ(ym‚θ)∥22σ2)(1)p(xn)=∑i=1M+1P(m)p(xn|m)=w1N+(1−w)∑i=1MP(m)1(2πσ2)D/2exp(−∥xn−T(ym‚θ)∥22σ2)(1)其中:P(m)为各个高斯模型的权重,简化起见,规定P(m)=1/M;w为无匹配点在点集中所占的比例;σ2为高斯模型的协方差。xn与ym间的对应关系由后验概率计算,即P(m|xn)=P(m)p(xn|m)/p(xn)(2)式(1)中,未知参数为θ、σ2。构造目标函数,当函数具有最小值时的参数取值为最优解,即E(θ‚σ2)=-Ν∑n=1logΜ+1∑m=1Ρ(m)p(xn|m)(3)E(θ‚σ2)=−∑n=1Nlog∑m=1M+1P(m)p(xn|m)(3)改写式(3),利用期望最大化方法(EM)算法求解Q(θ‚σ2)=12σ2Μ‚Ν∑m‚n=1Ρ(m|xn)∥xn-(ym+Τ(ym))∥2+ΝΡD2logσ2(4)Q(θ‚σ2)=12σ2∑m‚n=1M‚NP(m|xn)∥xn−(ym+T(ym))∥2+NPD2logσ2(4)其中,Νp=Ν∑n=1Μ∑m=1Ρ(m|xn)Np=∑n=1N∑m=1MP(m|xn)对于刚体变换,式(4)中T可进一步分解为平移、旋转和缩放等矩阵;而对于非刚体变换则需要选择更为复杂的函数形式v,对形变进行满足平滑性和准确性的逼近。对式(4)加入正则化项,进行平滑性约束,有Q(v‚σ2)=12σ2Μ‚Ν∑m‚n=1Ρold(m|xn)∥xn-(ym+v(ym))∥2+ΝΡD2logσ2+λ2∥Ρv∥(5)Q(v‚σ2)=12σ2∑m‚n=1M‚NPold(m|xn)∥xn−(ym+v(ym))∥2+NPD2logσ2+λ2∥Pv∥(5)根据MCT理论,相邻点的运动是连贯的,即运动场v是平滑的,若将v用参数为β的高斯核的线性组合表示,则正好满足这种平滑性要求,即v(yi)=Μ∑m=1wmexp(-12∥yi-ymβ∥)=1σ2λΝ∑n=1Ρ(mxn)(xn-(ym+v(ym))exp(-12)∥yi-ymβ∥)(6)v(yi)=∑m=1Mwmexp(−12∥yi−ymβ∥)=1σ2λ∑n=1NP(mxn)(xn−(ym+v(ym))exp(−12)∥yi−ymβ∥)(6)求解最小二乘问题及变分方程,得到v及σ2的最优解;再根据式(2),计算点集间的对应关系矩阵,找到使得后验概率最大的对应点对,从而得到点集Y到X的形变场。4结果与分析4.1oxel1333实验数据为64层螺旋CT机扫描图像,分为11个期相,每个期相的数据均为512voxel×512voxel×215voxel,层间分辨率为0.625mm,层内分辨率为0.3899mm×0.3899mm。图1是单期相分割出的左冠脉及细化后的骨架。将多期相的左、右冠脉都分割出来,并进行细化,骨架线如图2所示。4.2位移估计误差为了定量评价配准算法精度,对某一时刻的回旋支LCX骨架构造已知参数的运动场,得到形变后的曲线,然后对这两条曲线进行配准进而估计运动场,比较估计出的运动场与已知运动场的参数,评价算法的误差。心脏的运动主要为非刚体的膨胀和收缩,可近似用“旋转+平移+放大/收缩+扭转”模拟。假设Rz、Ry和Rx分别表示绕三维空间z、y、x3个轴的旋转矩阵,E表示缩放矩阵,S表示扭转变形矩阵,T为平移矩阵,则某点y(y1,y2,y3)运动后的坐标x(x1,x2,x3)可表示为[x1x2x3]=RzRyRxES[y1y2y3]+Τ=[cosγsinγ0-sinγcosγ0001][cosβ0-sinβ010sinβ0cosβ][1000cosαsinα0-sinαcosα][Ex000Ey000Ez][1Κxx-Κyy-Κxx1ΚzzΚyy-Κzz][y1y2y3]+[ΤxΤyΤz](7)式(7)中,各项参数的取值见表1。膨胀、收缩时构造出的运动场及估计出的运动场如图3、图4所示。直观上,图3、图4估计出的运动场都是比较平滑和均匀的,且和构造的运动场比较相似。对两种情况下位移情况进行计算,其中Drmax、Drmin、Dravex、Demax、Demin和Deavex分别表示构造的运动场与估计出的运动场的最大位移、最小位移和平均位移,见表2。以最大位移为例,定义最大位移估计误差为Emax=Demax-DrmaxDrmax×100％(8)其他两种位移误差定义与式(8)类似,配准误差以均方根误差(RMSE)衡量。两种情况下的配准RMSE分别为0.4744mm和0.3512mm,小于1voxel,位移估计误差小于1%。由计算结果可知,估计出的运动场与构造的运动场差别很小,采用CPD进行运动估计,具有较高的精度。4.3降压降位计算时各各自为基的平均移动距离、平均位移幅度及速率对全期相的左右冠脉主要分支前降支LAD、回旋支LCX、对角支D及右冠脉RCA分别进行配准,并计算其在各个时刻的3D位移幅度和速率。位移幅度定义为其他时刻相对于R-R间期0%时各分支上点的平均移动距离;速度定义为Vi=Di/t(i=1,2,……10),其中Di是该期相到下一期相的平均位移幅度,t为期相间的时间间隔,以心率60bpm计算,则t=0.1s。计算结果见表3、表4。对表3、4中的数据进行三次样条插值,得到整个R-R间期内任意时刻的各分支的3D平均位移幅度及速率,如图5、图6所示。从图可以看出,右冠脉相对于左冠脉运动更为明显,具有更大的位移幅度及更高的速率,同时左、右冠脉均有3个速率较高的时刻,对应心脏运动中的收缩、舒张期充盈及心房收