CN104240287A - 一种利用ct图像生成冠脉全景图的方法及系统 - Google Patents
一种利用ct图像生成冠脉全景图的方法及系统 Download PDFInfo
- Publication number
- CN104240287A CN104240287A CN201310228555.7A CN201310228555A CN104240287A CN 104240287 A CN104240287 A CN 104240287A CN 201310228555 A CN201310228555 A CN 201310228555A CN 104240287 A CN104240287 A CN 104240287A
- Authority
- CN
- China
- Prior art keywords
- mtd
- image
- model
- mrow
- standard
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000013507 mapping Methods 0.000 claims abstract description 63
- 230000000747 cardiac effect Effects 0.000 claims description 100
- 230000009466 transformation Effects 0.000 claims description 47
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 43
- 238000005070 sampling Methods 0.000 claims description 31
- 210000003484 anatomy Anatomy 0.000 claims description 21
- 210000004351 coronary vessel Anatomy 0.000 claims description 16
- 238000013519 translation Methods 0.000 claims description 7
- 238000005259 measurement Methods 0.000 claims description 4
- 210000000709 aorta Anatomy 0.000 claims description 3
- 210000005246 left atrium Anatomy 0.000 claims description 3
- 210000005240 left ventricle Anatomy 0.000 claims description 3
- 210000004165 myocardium Anatomy 0.000 claims description 3
- 210000001147 pulmonary artery Anatomy 0.000 claims description 3
- 210000005245 right atrium Anatomy 0.000 claims description 3
- 210000005241 right ventricle Anatomy 0.000 claims description 3
- 230000004069 differentiation Effects 0.000 abstract 1
- 238000002591 computed tomography Methods 0.000 description 55
- 238000010586 diagram Methods 0.000 description 18
- 238000005457 optimization Methods 0.000 description 6
- 238000005192 partition Methods 0.000 description 4
- 230000011218 segmentation Effects 0.000 description 3
- 238000000638 solvent extraction Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 210000003516 pericardium Anatomy 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000011524 similarity measure Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 208000029078 coronary artery disease Diseases 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000004092 self-diagnosis Methods 0.000 description 1
Landscapes
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开了一种利用CT图像生成冠脉全景图的方法及系统,该方法步骤:获取标准心脏模型,标准心脏模型区分为多个解剖区域;获取全景模型;建立标准心脏模型表面各点与全景模型的表面各点的映射关系,将标准心脏模型表面各点映射到全景模型的表面的对应位置,建立标准全景模型;将待分析心脏CT图像与标准心脏模型进行配准,使得配准后的待分析心脏CT图像具有与标准心脏模型相同的空间位置关系;根据配准后的待分析心脏CT图像与标准心脏模型的空间位置映射关系,及,该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将配准后的待分析心脏CT图像映射到该标准全景模型上,从而获得冠脉全景图和/或待分析心脏的解剖区域区分图。
Description
技术领域
本发明涉及医学图像的处理,特别是涉及一种利用CT图像生成冠脉全景图的方法及系统。
背景技术
当前,为了获取医学图像,大量利用了X射线断层扫描技术(computedtomography,CT)系统,利用该系统可获得多切片的CT图像数据。
在现有技术中,医生通常利用由病人心脏的多切片CT图像数据形成的如图7所示的病人的3D冠脉图进行冠脉病的诊断,以获取病人的冠脉情况。
然而,3D冠脉图的缺点在于,无法将整个冠脉树显示在以心脏为背景的基面上,无法非常直观地了解冠状动脉的走行、分布和形态。
发明内容
本发明解决的技术问题在于,利用待分析心脏的CT图像生成冠脉全景图。
更进一步的,建立待分析心脏的CT图像与标准全景模型之间的映射。进而在该标准全景模型上显示整个冠脉树,方便医生更清晰直观的查看冠脉树。
本发明公开了一种利用CT图像生成冠脉全景图的方法,包括如下步骤:
步骤A、获取标准心脏模型,该标准心脏模型区分为多个解剖区域;
步骤B、获取全景模型;
步骤C、建立该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将标准心脏模型表面各点映射到该全景模型的表面的对应位置,从而建立标准全景模型;
步骤D、将待分析心脏CT图像与该标准心脏模型进行配准,使得配准后的待分析心脏CT图像具有与该标准心脏模型相同的空间位置关系;
步骤E、根据配准后的待分析心脏CT图像与该标准心脏模型的空间位置映射关系,以及,该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将配准后的待分析心脏CT图像映射到该标准全景模型上,从而获得冠脉全景图和/或待分析心脏的解剖区域区分图。
该解剖区域包括左心房、右心房、左心室、右心室、心肌、升主动脉和肺动脉。
该步骤C中建立该标准心脏模型表面各点与该全景模型的表面各点的映射关系的步骤进一步包括:
步骤C1,将该标准心脏模型的表面进行网格划分;
步骤C2,根据划分得到的网格分布位置以及全景模型的表面的像素点的分布位置,建立该标准心脏模型的每个网格与该全景模型的表面的像素点的映射关系;
步骤C3,根据该映射关系,将标准心脏模型的每个网格的图像,对应映射到该全景模型的表面,从而建立标准全景模型。
步骤C1进一步包括:
步骤C11,对该标准心脏模型的表面进行初始化网格划分;
步骤C12,利用自适应网格生成法,对初始化网格进行细分。
步骤C12进一步包括:
步骤C121,计算每个初始化网格的一特征曲率值;
步骤C122,依次判断每个该初始化网格的特征曲率值是否符合一阈值区间,对于不符合该阈值区间的初始化网格执行步骤C123,直到所有网格的特征曲率值均符合该阈值区间;
步骤C123,在该网格中设置一新节点,利用该新节点,分割该网格。
该步骤C进一步包括:将该标准全景模型进行二维平面映射,获得该标准全景模型的二维图。
步骤D进一步包括:
将该待分析心脏CT图像与该标准心脏模型进行仿射配准;
基于仿射配准的结果,将待分析心脏CT图像与该标准心脏模型进行弹性配准。
该仿射配准进一步包括:
步骤D1,对待分析心脏CT图像与标准心脏模型的图像进行等分辨率采样;
步骤D2,根据仿射变换公式将待分析心脏CT图像的采样图像中的每个像素点映射到由标准心脏模型的采样图像所建立的物理坐标系中;
步骤D3,根据映射结果获取MSD测度,对MSD测度进行寻优,来确定最优的仿射变换公式;
步骤D4,根据该最优的仿射变换公式,对待分析心脏CT图像的采样图像进行仿射变换。
该弹性配准进一步包括:
步骤D5,设置B样条的采样率,将步骤D4中经仿射变换后的待分析心脏CT图像的采样图像作为B样条弹性配准的初始图像,获得控制点的位置;
步骤D6,根据每个控制点所控制的局部图像,建立局部图像互信息测度;
步骤D7,根据局部图像互信息测度,采用LBFGS算法对每个控制点的位置进行优化;;
步骤D8,根据优化的控制点的位置,对步骤D4中经仿射变换后的待分析心脏CT图像的采样图像进行线性插值,得到弹性配准的结果图像。
该仿射变换公式为:
其中,
θx,θy,θz分别表示围绕X,Y,Z坐标轴的连续旋转角度,T为沿着三个轴的平移量,t1,t2,t3分别表示沿着三个坐标轴的平移距离,s1,s2,s3分别表示沿着三个坐标轴的缩放因子,x表示待分析心脏CT图像的采样图像中每个坐标点,y表示该坐标点经过仿射变换后映射到标准心脏模型的采样图像所建立的物理坐标系中的位置,CenterF为待分析心脏CT图像的采样图像的物理中心坐标,CenterM为标准心脏模型的采样图像的物理中心坐标;
MSD(IF,IM)代表该MSD测度,IF表示待分析心脏CT图像的采样图像,IM表示标准心脏模型的采样图像,T是该仿射变换公式,Ω表示IF中像素点的总个数。
该互信息测度为:
IF表示待分析心脏CT图像的采样图像中每个控制点控制范围内的图像,H(IF)表示IF范围内图像信息熵,IM表示标准心脏模型的采样图像中每个控制点控制范围内的图像,H(IM)表示IM范围内图像信息熵,H(IF,IM)表示两幅图像的联合熵。
本发明还公开了一种利用CT图像生成冠脉全景图的系统,包括:
标准心脏模型获取单元,用于获取标准心脏模型,该标准心脏模型区分为多个解剖区域;
全景模型获取单元,用于获取全景模型;
标准全景模型建立单元,用于建立该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将标准心脏模型表面各点映射到该全景模型的表面的对应位置,从而建立标准全景模型;
配准单元,用于将待分析心脏CT图像与该标准心脏模型进行配准,使得配准后的待分析心脏CT图像具有与该标准心脏模型相同的空间位置关系;
冠脉全景图获取单元,用于根据配准后的待分析心脏CT图像与该标准心脏模型的空间位置映射关系,以及,该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将配准后的待分析心脏CT图像映射到该标准全景模型上,从而获得冠脉全景图和/或待分析心脏的解剖区域区分图。
本发明的技术效果在于,获取使冠脉全景图。另外,还可获得自动区分出各个解剖区域的心脏的分区结果示意图。进一步的,本发明可提高冠脉全景图的精确性。
附图说明
图1所示为本发明的利用CT图像生成冠脉全景图的方法的流程图。
图2所示为本发明的心脏模型的MIP图。
图3所示为区分各解剖区域的标准心脏模型的示意图。
图4所示为标准心脏模型经网格化后的结果示意图。
图5A所示为待分析的心脏CT图像。
图5B所示为标准心脏模型的二维示意图。
图5C为仿射变换的配准结果示意图。
图5D为对图5C进行弹性配准的初始B样条的网格位置示意图。
图5E、5F,分别为配准完成后B样条网格位置以及弹性配准结果示意图。
图6A、6B、6C为心脏的分区结果示意图,其中,图6A为轴断图、图6B为矢状图、图6C为冠状图。
图7为病人的3D冠脉图。
图8A为个体病人的冠脉全景图的二维图,图8B为个体病人的冠脉全景图的三维图。
图9A为纹理图像示意图。
图9B为映射区域像示意图。
图10为映射关系示意图。
具体实施方式
如图1所示为本发明的利用CT图像生成冠脉全景图的方法的流程图。其中包括:
步骤A、获取标准心脏模型,该标准心脏模型已区分为多个解剖区域;
步骤B、获取全景模型;
步骤C、建立该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将标准心脏模型表面各点映射到该全景模型的表面的对应位置,从而建立标准全景模型;
步骤D、将待分析心脏CT图像与该标准心脏模型进行配准,使得配准后的待分析心脏CT图像具有与该标准心脏模型相同的空间位置关系;
步骤E、根据配准后的待分析心脏CT图像与该标准心脏模型的空间位置映射关系,以及,该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将配准后的待分析心脏CT图像映射到该标准全景模型上,从而获得冠脉全景图和/或待分析心脏的解剖区域区分图。
以下具体对上述步骤进行详细说明。
为了获取冠脉全景图,首先需要获取标准心脏模型。获取标准心脏模型的步骤进一步包括建立心脏模型的步骤以及对心脏模型进行分区的步骤。
具体而言,为了获取标准心脏模型,需大量的搜集、筛选、整理得到多例心脏CT图像数据。随后,利用手动分割的方式,从该多例心脏CT图像数据中分别提取出心脏。将提取出的多例心脏进行合并,以建立起心脏模型。该心脏模型是建立标准心脏模型的重要基础,如图2所示为心脏模型的MIP图。该心脏模型为三维模型。
建立心脏模型的具体步骤包括:对每一例心脏的解剖切片进行手工分割,提取心包和各腔室等解剖区域的轮廓,并用高斯滤波器对其进行平滑后,使用均匀取样的方法对轮廓进行取样,得到的样点作为心包和腔室的拟合点,通过B样条的曲面拟合获得完整轮廓。由于不同个体的心脏的构造具有差异性,则心脏模型的建立是一个大量且长期的数据计算、合并、校正的过程。在这个过程中需综合空间、自由度、概率等约束条件。心脏模型的建立过程属于本领域的技术人员能够掌握的。
随后,通过手动分割的方式对该心脏模型进行分区,分区后的心脏模型即为标准心脏模型。具体来说,手动将心脏模型分割为7个解剖区域,即左心房,右心房,左心室,右心室,心肌,升主动脉和肺动脉。如图3所示为包含各解剖区域的标准心脏模型的示意图。对于步骤B,医生在诊疗过程中,从给定的全景模型模板中选择有利于自己诊断的全景模型。例如,可选择球体或椭圆体的全景模型。
该步骤C进一步包括:
步骤C1,在获得标准心脏模型后,对该标准心脏模型的表面进行网格划分;
步骤C2,根据网格的分布位置以及该全景模型的表面的像素点的分布位置,建立该标准心脏模型表面各网格与该全景模型的表面各点的映射关系;
步骤C3,根据该映射关系,将标准心脏模型的表面的图像,网格映射到该全景模型的表面,从而建立标准全景模型。
步骤C1具体包括:
步骤C11,依一预设值对该标准心脏模型的表面进行初始化网格划分。
标准心脏模型的表面呈曲面状,依该预设值进行初始化网格划分后得到的初始化网格较大,较为粗糙。此时不适于直接进行与全景模型的映射。为了提高后续映射的精度,需要进一步进行网格细化。依预设值形成的初始化网格可为四边形或三角形或圆形或多边形等,通过预设值可确定不同的网格形状。
每个初始化网格具有代表该网格的坐标值,以及,用于记录与相邻网格之间的相邻关系的单元联系数组。该坐标值例如是该网格的中心点的坐标。
步骤C12,利用自适应网格生成法,对初始化网格进行细分。
细分即为将每一个网格划分成多个小网格,使得曲面网格划分更加细致,更加符合曲面曲率变化。
步骤C12进一步包括:
步骤C121,计算每个初始化网格的一特征曲率值;
由于该标准心脏模型的表面为曲面,而该曲面采用均匀三次B样条[2]表达,本领域的技术人员利用参数表达曲面的曲率公式可求得曲面上任一点的曲率。
在本发明中,可选择每个初始化网格中的至少一点的曲率值作为该特征曲率值,也可以采用每个初始化网格中的多点的平均值作为该特征曲率值,该特征曲率值可根据需要而定义。
步骤C122,依次判断每个该初始化网格的特征曲率值是否符合一阈值区间,对于不符合阈值区间的初始化网格执行步骤C123,直到所有当前网格的特征曲率值均符合该阈值区间时,继续执行步骤C2。
为了满足曲面的精度不同需求,可调整该阈值区间的取值范围,从而实现不同的容许误差要求。当所有网格的特征曲率值均符合该阈值区间时,意味着当前划分得到的所有网格均满足了容许误差要求,具备了进行建立后续映射关系的精度条件。
在步骤C122中,在对所有初始化网格的特征曲率值是否符合阈值区间的判断结束后,只要有一个不符合,针对不符合的那些初始化网格,启动步骤C123的细分过程,针对这些初始化网格的细分过程结束后,对当前存在的全部网格再次进行特征曲率值的计算并判断其是否符合该阈值区间,如果还有网格不符合,继续执行步骤C123,如此循环,直至当前存在的全部网格的特征曲率值均符合该阈值区间。
步骤C123,针对需要细分的网格,在该网格中设置一新节点,利用该新节点,分割该需要细分的网格。
设置新节点可依据预设规则进行,该预设规则可包括:1、将该网格的中心点作为新节点;2、将该网格的对角线交点作为新节点;3、以该网格的一边作为底边,向网格内部形成一等边三角形,若该等边三角形的顶点处于该网格内部,将该顶点作为该网格的新节点。
本领域的技术人员可根据需求以及网格的形状,具体选择该预设规则。例如,对于圆形的网格可选择将该网格的中心点作为新节点,对于方形的网格可选择将该网格的对角线交点作为新节点,对于梯形的网格可选择上述第3种预设规则选择新节点。
在新节点生成后,基于该新节点并利用预设划分规则,进行网格划分。例如,对于方形网格,该预设划分规则可包括,将该方形网格的原节点、与该原节点相连接的两条边的中点以及该新节点,作为划分后的新网格的四个节点。图4所示为标准心脏模型经网格化后的结果示意图。
所有新网格产生后,计算每个该新网格的坐标值,更新标准心脏模型的表面所包括的当前所有网格的单元联系数组,并且,针对整个曲面的此次网格划分,更新对应当前曲面的参数区域,基于更新后的参数区域,可为曲面上当前的每个网格计算其特征曲率值。
该参数区域可体现出曲率计算的精确程度。设初始化网格划分后,点Qi处曲率为Ai(i=1,…,m),则当时所设定的参数范围U∈[0,1],
每网格划分一次,参数区域也进行对应的更新,此时,点Qi所对应的参数范围Ui更新为Ui=Ti/T,其中,Ti=1/Ai,在更新后的该参数区域Ui中,可进一步计算网格特征曲率值,其得到的特征曲率值的精度更高。
步骤C2进一步包括:
根据当前所生成的所有网格的网格总数以及全景模型表面的像素点总数,基于所有网格的分布位置以及全景模型的表面的像素点的分布位置,建立起该标准心脏模型的表面的每个网格与该全景模型的表面的像素点的映射关系。
具体的,例如将球体的全景模型区分为八个区域,例如首先区分上下半球,每个半球再区分为四个区域(左上、左下、右上、右下),也就是说,全景模型的表面被区分为八个区域,同样的,将标准心脏模型的表面依照同样方式也区分为八个区域,建立起该标准心脏模型的表面的每个区域与该全景模型的表面的每个区域的对应关系。例如,两种表面的上半球左上区域之间存在对应关系。依照同样原理,参照同样的坐标系为每个区域进一步细化,从而为对应的分布位置的网格和全景模型的表面的像素点起建立映射关系。
全景模型的表面的多个像素点可能对应标准心脏模型的表面的一个网格。
在步骤C3中,将标准心脏模型的每个网格的图像映射到全景模型的表面与该网格对应的像素上。
映射以后,每个网格原本具有的单元联系数组及节点编号不变,只是每个网格的节点坐标发生改变,由全景模型的表面的一空间坐标值代替。
如此一来,从全景模型的表面的每个坐标点,可以获取对应的标准心脏模型的网格的数据。全景模型经过该映射过程,成为标准全景模型。该映射过程直接生成标准全景模型的三维模型,也可以根据该标准全景模型的三维模型,生成该标准全景模型的二维模型,供医生参考。
生成该标准全景模型的二维模型的具体方法包括:
步骤C41,将该标准全景模型的三维模型的表面展开;
步骤C42,指定一参考平面;
步骤C43,确定标准全景模型的三维模型的表面与该参考平面的相对位置关系;
步骤C44,根据映射关系,将全景模型的三维表面映射到二维平面上。
求解三角网络模型变型设计的参考平面,将网格定点向参考平面投影,建立网格顶点与投影点间的映射关系,通过调整参考平面控制顶点,根据投影点位移量调整网格定点坐标,采用平面理论进行网格模型光顺处理,获得二维平面图。
例如,图9A为纹理图像,其为该标准全景模型的三维模型的表面的示意图。图9B为映射区域,纹理图像向参考平面投影后生成的,也就是该标准全景模型的二维模型。
图9A的纹理图像的宽txt_width映射为图9B中y轴上的值,图9A的纹理图像的长txt_length映射为图9B中x轴上的值。
映射图上的O,X,Y位置已经从纹理图像,通过旋转算子矩阵、平移向量等关系获得(length,width,)-->(x,y,)。
对于步骤D,将待分析心脏CT图像与标准心脏模型进行配准,配准的目的是使得待分析心脏CT图像与标准心脏模型进行空间对准,使得待分析心脏CT图像具备了与标准心脏模型同样的空间位置。
步骤D大体上分为两个部分:第一部分,将待分析心脏CT图像与该标准心脏模型进行仿射配准;第二部分,基于仿射配准的结果,将待分析心脏CT图像与该标准心脏模型进行弹性配准。
也就是说,步骤D的配准过程中,采取的是一个全局到局部的配准过程,首先以仿射变换对待分析心脏CT图像与标准心脏模型进行全局配准,使得待分析心脏CT图像初步具备与标准心脏模型相同的空间位置关系,然后以B样条弹性配准对图像的局部进行更精确配准,提高配准结果的精确程度。
具体而言,本发明的两种配准均为首先预设一空间变换参数,根据空间变换参数进行空间变换,其次建立一配准成本函数对空间变换结果进行描述或者衡量。空间变换即用于空间对准。本发明即通过建立一个配准成本函数来评价空间对准的正确性,然后通过最优化算法寻找配准成本函数的最小值(或最大值),将与该最小值对应的空间变换参数确定为最优的空间变换参数。
公式1为该配准成本函数的通用形式。这里IF表示待分析心脏CT图像,IM表示标准心脏模型的图像,T表示空间变换参数,表示最优空间变换参数,C是测度函数,表示对象图与模板图像之间的相似度,即,用来测量待分析心脏CT图像与标准心脏模型之间的相似度。具体来说,通过建立测度函数,使空间变换参数的变化与两幅图像相似度变化之间建立联系,从而通过对相似性测度寻优来确定最优化空间变换
步骤D具体包括:
步骤D1、图像采样:对待分析心脏CT图像与标准心脏模型的图像进行等分辨率采样,使两组图像具有相同的物理分辨率。
例如,对两组图像在X,Y,Z上分别以1.0mm进行采样。即图像体像素间隔(Spacing)为1.0mm*1.0mm*1.0mm,采样得到的待分析心脏CT图像与标准心脏模型的采样图像,分别记为ImageF和ImageM。
步骤D2、根据ImageF在ImageM坐标系中的仿射变换公式,将ImageF中的每个像素点映射到由ImageM所建立的物理坐标系中。
利用ImageF和ImageM的物理中心建立仿射变换公式,从而使ImageF中的每个像素点映射到由ImageM所建立的物理坐标系中。
物理中心计算公式为:
Cx=ImageColumns/2
Cy=ImageRows/2
Cz=ImageLayers/2 公式2
其中:Cx,Cy,Cz为该图像的物理中心,ImageColumns表示该图像的列数,ImageRows表示该图像的行数,ImageLayers表示该图像的层数。
仿射变换公式为:
公式3
其中,
θx,θy,θz分别表示围绕X,Y,Z坐标轴的连续旋转角度,T为沿着三个轴的平移。t1,t2,t3分别表示沿着三个坐标轴的平移距离,s1,s2,s3分别表示沿着三个坐标轴的缩放因子。x表示ImageF中每个坐标点,y表示该坐标点经过仿射变换后映射到ImageM所建立的物理坐标系中的位置。CenterF为ImageF的物理中心坐标,CenterM为ImageM的物理中心坐标(由公式1计算得出)。对于待分析心脏CT图像上的每个像素点,通过公式3进行仿射变换后,都可以映射到以ImageM物理中心所建立的坐标系中。仿射变换公式即为一种空间变换参数。
步骤D3、根据该步骤D2的映射结果获取MSD测度,对MSD测度进行寻优,来确定最优的仿射变换公式。
建立两组图像的仿射变换关系后,需要对公式3中的9个变量参数θx,θy,θz,t1,t2,t3,s1,s2,s3进行优化,这里使用MSD(Mean of squaredDifferernces)作为相似性测度,即通过计算两组图像中所有对应点像素值的差方合来表示图像的相似度(如公式4所示),公式4是公式1中测度函数C的一个具体实例。
这里IF表示ImageF,IM表示ImageM,T是仿射变换公式(公式3),Ω表示IF中像素点的总个数。当前,MSD测度公式与仿射变换公式为已知,此时利用Powell算法对公式3中的9个变量参数进行寻优,Powell算法对每每一个变量参数进行优化时,均利用梯度下降算法不断进行迭代,在达到终止条件(使MSD测度达到最小值)后,与该最小值所对应的变量参数为最优的变量参数,从而组成最优的仿射变换矩阵。
步骤D4、根据该最优的仿射变换公式,对ImageF进行仿射变换。
如图5A所示ImageF的二维示意图,图5B所示为ImageM的二维示意图,图5C为执行步骤D4后,对图5A进行仿射变换后的配准结果示意图ImageTF。
步骤D5、根据仿射配准的结果执行弹性配准过程:
以ImageTF作为与ImageM进行步骤D5的B样条弹性配准的初始图像。
弹性配准过程采用二层多分辨率的方法(依次设定B样条控制点的采样率分别为32mm*32mm*32mm、16mm*16mm*16mm。即,执行完毕一次配准后,调整采样率,再次执行弹性配准),多分辨率的方法可以在保证配准精度的同时有效的减少配准的时间。
在B样条弹性配准过程中,主要需要优化的对象是控制点的位置。图5D为对图5C进行弹性配准的初始B样条的网格位置示意图。控制点的初始位置如图5D中网格的交点,当控制点位置改变后,控制点所控制范围内的所有像素点均需要根据控制点的变换而进行插值,从而改变了控制点所控制范围内图像的灰度,所以通过对控制点位置进行改变从而控制影响范围内图像的灰度值,从而ImageTF与ImageM达到配准。本发明采用互信息测度(MI)对每个控制点范围内图像的相似度进行评价(公式5)。
上面的公式中IF表示ImageTF中每个控制点控制范围内的图像,H(IF)表示该范围内图像信息熵,IM表示ImageM中每个控制点控制范围内的图像,H(IM)表示该范围内图像信息熵,H(IF,IM)表示两幅图像的联合熵。
步骤D5包括如下具体过程:
步骤D5.1、设置B样条的采样率,并计算每个控制点所控制的图像范围。
根据三次B样条理论,每个控制点的控制范围为控制点周围的4个网格内的像素点。三维情况下,每个控制点的控制范围为64个网格内的像素点。步骤D5.2、根据每个控制点所控制的局部图像,利用公式5建立局部图像互信息测度MI。每个控制点位置的变换影响到局部图像灰度的变化,从而决定互信息MI的大小。
步骤D5.3、设置LBFGS优化的终止条件,采用LBFGS(拟牛顿法)算法对每个控制点的位置进行优化,以MI的值作为优化测度,通过更改控制点的位置从而对局部图像进行配准。
步骤D5.4、当优化算法终止后,缩小B样条率,进行新一轮的优化,重复步骤D5.1-D5.3。
步骤D5.5、当二次采样率都计算完后,根据最终控制点位置,对ImageTF进行线性插值从而确定最后的输出图像,参见图5E、5F,分别为配准完成后B样条网格位置以及弹性配准结果示意图。
步骤E:
由于步骤C已经获取了该标准心脏模型表面各点与该标准全景模型的表面各点的映射关系,步骤D也已经获取了该标准心脏模型表面各点与待分析心脏CT图像之间的映射关系,故而,利用该标准心脏模型作为媒介,在标准全景模型的表面各点,均可找到与该待分析心脏CT图像的对应关系,故而,利用这种对应关系,可将配准后的待分析心脏CT图像映射到该标准全景模型上,获得冠脉全景图。该冠脉全景图中包括标准心脏模型的各个解剖区域。
图8A为个体病人的冠脉全景图的二维图,图8B为个体病人的冠脉全景图的三维图。本发明得出的病人冠脉全景图可以将整个冠脉树完整的显示在一个平面上,方便医生查看冠脉。
另外,由于标准心脏模型已区分为多个解剖区域,故而,在步骤E的将配准后的待分析心脏CT图像映射到该标准全景模型上之后,除了可输出二维和三维的冠脉全景图,还可输出待识别心脏的分区结果示意图(待分析心脏的解剖区域区分图),其中图6A为轴状图、图6B矢状图、图6C冠状图。
本发明的技术效果在于,通过生成冠脉全景图,将整个冠脉树显示在以心脏为背景的基面上,从而非常直观地了解冠状动脉的走行、分布和形态。以及,获得自动区分出各个解剖区域的心脏的分区结果示意图。
Claims (12)
1.一种利用CT图像生成冠脉全景图的方法,其特征在于,包括如下步骤:
步骤A、获取标准心脏模型,该标准心脏模型区分为多个解剖区域;
步骤B、获取全景模型;
步骤C、建立该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将标准心脏模型表面各点映射到该全景模型的表面的对应位置,从而建立标准全景模型;
步骤D、将待分析心脏CT图像与该标准心脏模型进行配准,使得配准后的待分析心脏CT图像具有与该标准心脏模型相同的空间位置关系;
步骤E、根据配准后的待分析心脏CT图像与该标准心脏模型的空间位置映射关系,以及,该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将配准后的待分析心脏CT图像映射到该标准全景模型上,从而获得冠脉全景图和/或待分析心脏的解剖区域区分图。
2.如权利要求1所述的方法,其特征在于,该解剖区域包括左心房、右心房、左心室、右心室、心肌、升主动脉和肺动脉。
3.如权利要求所述1的方法,其特征在于,该步骤C中建立该标准心脏模型表面各点与该全景模型的表面各点的映射关系的步骤进一步包括:
步骤C1,将该标准心脏模型的表面进行网格划分;
步骤C2,根据划分得到的网格分布位置以及全景模型的表面的像素点的分布位置,建立该标准心脏模型的每个网格与该全景模型的表面的像素点的映射关系;
步骤C3,根据该映射关系,将标准心脏模型的每个网格的图像,对应映射到该全景模型的表面,从而建立标准全景模型。
4.如权利要求3所述的方法,其特征在于,步骤C1进一步包括:
步骤C11,对该标准心脏模型的表面进行初始化网格划分;
步骤C12,利用自适应网格生成法,对初始化网格进行细分。
5.如权利要求3所述的方法,其特征在于,步骤C12进一步包括:
步骤C121,计算每个初始化网格的一特征曲率值;
步骤C122,依次判断每个该初始化网格的特征曲率值是否符合一阈值区间,对于不符合该阈值区间的初始化网格执行步骤C123,直到所有网格的特征曲率值均符合该阈值区间;
步骤C123,在该网格中设置一新节点,利用该新节点,分割该网格。
6.如权利要求1所述的方法,其特征在于,该步骤C进一步包括:将该标准全景模型进行二维平面映射,获得该标准全景模型的二维图。
7.如权利要求1所述的方法,其特征在于,步骤D进一步包括:
将该待分析心脏CT图像与该标准心脏模型进行仿射配准;
基于仿射配准的结果,将待分析心脏CT图像与该标准心脏模型进行弹性配准。
8.如权利要求7所述的方法,其特征在于,该仿射配准进一步包括:
步骤D1,对待分析心脏CT图像与标准心脏模型的图像进行等分辨率采样;
步骤D2,根据仿射变换公式将待分析心脏CT图像的采样图像中的每个像素点映射到由标准心脏模型的采样图像所建立的物理坐标系中;
步骤D3,根据映射结果获取MSD测度,对MSD测度进行寻优,来确定最优的仿射变换公式;
步骤D4,根据该最优的仿射变换公式,对待分析心脏CT图像的采样图像进行仿射变换。
9.如权利要求8所述的方法,其特征在于,该弹性配准进一步包括:
步骤D5,设置B样条的采样率,将步骤D4中经仿射变换后的待分析心脏CT图像的采样图像作为B样条弹性配准的初始图像,获得控制点的位置;
步骤D6,根据每个控制点所控制的局部图像,建立局部图像互信息测度;
步骤D7,根据局部图像互信息测度,采用LBFGS算法对每个控制点的位置进行优化;
步骤D8,根据优化的控制点的位置,对步骤D4中经仿射变换后的待分析心脏CT图像的采样图像进行线性插值,得到弹性配准的结果图像。
10.如权利要求8所述的方法,其特征在于,该仿射变换公式为:
其中,
θx,θy,θz分别表示围绕X,Y,Z坐标轴的连续旋转角度,T为沿着三个轴的平移量,t1,t2,t3分别表示沿着三个坐标轴的平移距离,s1,s2,s3分别表示沿着三个坐标轴的缩放因子,x表示待分析心脏CT图像的采样图像中每个坐标点,y表示该坐标点经过仿射变换后映射到标准心脏模型的采样图像所建立的物理坐标系中的位置,CenterF为待分析心脏CT图像的采样图像的物理中心坐标,CenterM为标准心脏模型的采样图像的物理中心坐标;
MSD(IF,IM)代表该MSD测度,IF表示待分析心脏CT图像的采样图像,IM表示标准心脏模型的采样图像,T是该仿射变换公式,Ω表示IF中像素点的总个数。
11.如权利要求9所述的方法,其特征在于,该互信息测度为:
IF表示待分析心脏CT图像的采样图像中每个控制点控制范围内的图像,H(IF)表示IF范围内图像信息熵,IM表示标准心脏模型的采样图像中每个控制点控制范围内的图像,H(IM)表示IM范围内图像信息熵,H(IF,IM)表示两幅图像的联合熵。
12.一种利用CT图像生成冠脉全景图的系统,其特征在于,包括:
标准心脏模型获取单元,用于获取标准心脏模型,该标准心脏模型区分为多个解剖区域;
全景模型获取单元,用于获取全景模型;
标准全景模型建立单元,用于建立该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将标准心脏模型表面各点映射到该全景模型的表面的对应位置,从而建立标准全景模型;
配准单元,用于将待分析心脏CT图像与该标准心脏模型进行配准,使得配准后的待分析心脏CT图像具有与该标准心脏模型相同的空间位置关系;
冠脉全景图获取单元,用于根据配准后的待分析心脏CT图像与该标准心脏模型的空间位置映射关系,以及,该标准心脏模型表面各点与该全景模型的表面各点的映射关系,将配准后的待分析心脏CT图像映射到该标准全景模型上,从而获得冠脉全景图和/或待分析心脏的解剖区域区分图。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310228555.7A CN104240287B (zh) | 2013-06-08 | 2013-06-08 | 一种利用ct图像生成冠脉全景图的方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310228555.7A CN104240287B (zh) | 2013-06-08 | 2013-06-08 | 一种利用ct图像生成冠脉全景图的方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104240287A true CN104240287A (zh) | 2014-12-24 |
CN104240287B CN104240287B (zh) | 2017-10-20 |
Family
ID=52228289
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310228555.7A Expired - Fee Related CN104240287B (zh) | 2013-06-08 | 2013-06-08 | 一种利用ct图像生成冠脉全景图的方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104240287B (zh) |
Cited By (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105662451A (zh) * | 2016-03-31 | 2016-06-15 | 北京思创贯宇科技开发有限公司 | 一种主动脉图像处理方法和装置 |
CN106934821A (zh) * | 2017-03-13 | 2017-07-07 | 中国科学院合肥物质科学研究院 | 一种基于icp算法和b样条的锥形束ct和ct图像配准方法 |
CN107038706A (zh) * | 2017-05-16 | 2017-08-11 | 西安电子科技大学 | 基于自适应网格的红外图像置信度评估装置及方法 |
CN107087114A (zh) * | 2017-06-29 | 2017-08-22 | 深圳天珑无线科技有限公司 | 一种拍摄的方法及装置 |
CN108875780A (zh) * | 2018-05-07 | 2018-11-23 | 广东省电信规划设计院有限公司 | 基于图像数据的图像间差异对象的获取方法以及装置 |
CN108992193A (zh) * | 2017-06-06 | 2018-12-14 | 苏州笛卡测试技术有限公司 | 一种牙齿修复辅助设计方法 |
CN109427059A (zh) * | 2017-08-18 | 2019-03-05 | 西门子医疗有限公司 | 解剖结构的平面可视化 |
CN109584201A (zh) * | 2018-09-14 | 2019-04-05 | 新影智能科技(昆山)有限公司 | 医学图像配准方法、配准系统、存储介质、及电子设备 |
CN110197713A (zh) * | 2019-05-10 | 2019-09-03 | 上海依智医疗技术有限公司 | 一种医疗影像的处理方法、装置、设备和介质 |
CN111544020A (zh) * | 2020-04-17 | 2020-08-18 | 北京东软医疗设备有限公司 | X射线成像设备的几何校正方法及装置 |
CN112037186A (zh) * | 2020-08-24 | 2020-12-04 | 杭州深睿博联科技有限公司 | 一种基于多视图模型融合的冠脉血管提取方法及装置 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2008122056A2 (en) * | 2007-04-02 | 2008-10-09 | Case Western Reserve University | Medical apparatus and method associated therewith |
CN102208109A (zh) * | 2011-06-23 | 2011-10-05 | 南京林业大学 | X射线图像和激光图像的异源图像配准方法 |
CN102411780A (zh) * | 2011-09-07 | 2012-04-11 | 华南理工大学 | 一种基于配准的ct图像全心脏自动分割系统 |
CN103020976A (zh) * | 2012-12-31 | 2013-04-03 | 中国科学院合肥物质科学研究院 | 一种基于带权模糊互信息的三维医学图像配准方法及系统 |
CN103136739A (zh) * | 2011-11-29 | 2013-06-05 | 北京航天长峰科技工业集团有限公司 | 一种复杂场景下可控摄像机监控视频与三维模型配准方法 |
-
2013
- 2013-06-08 CN CN201310228555.7A patent/CN104240287B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2008122056A2 (en) * | 2007-04-02 | 2008-10-09 | Case Western Reserve University | Medical apparatus and method associated therewith |
CN102208109A (zh) * | 2011-06-23 | 2011-10-05 | 南京林业大学 | X射线图像和激光图像的异源图像配准方法 |
CN102411780A (zh) * | 2011-09-07 | 2012-04-11 | 华南理工大学 | 一种基于配准的ct图像全心脏自动分割系统 |
CN103136739A (zh) * | 2011-11-29 | 2013-06-05 | 北京航天长峰科技工业集团有限公司 | 一种复杂场景下可控摄像机监控视频与三维模型配准方法 |
CN103020976A (zh) * | 2012-12-31 | 2013-04-03 | 中国科学院合肥物质科学研究院 | 一种基于带权模糊互信息的三维医学图像配准方法及系统 |
Non-Patent Citations (2)
Title |
---|
XUESONG LU等: "《Mutual Information-Based Multimodal Non-Rigid Image Registration Using Free-Form Deformation with A New Joint Histogram Estimation》", 《2007 IEEE/ICME INTERNATIONAL CONFERENCE ON COMPLEX MEDICAL ENGINEERING》 * |
王娟等: "《一种柱面全景图像自动拼接算法》", 《计算机仿真》 * |
Cited By (15)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105662451A (zh) * | 2016-03-31 | 2016-06-15 | 北京思创贯宇科技开发有限公司 | 一种主动脉图像处理方法和装置 |
CN105662451B (zh) * | 2016-03-31 | 2018-10-09 | 北京思创贯宇科技开发有限公司 | 一种主动脉图像处理方法和装置 |
CN106934821B (zh) * | 2017-03-13 | 2020-06-23 | 中国科学院合肥物质科学研究院 | 一种基于icp算法和b样条的锥形束ct和ct图像配准方法 |
CN106934821A (zh) * | 2017-03-13 | 2017-07-07 | 中国科学院合肥物质科学研究院 | 一种基于icp算法和b样条的锥形束ct和ct图像配准方法 |
CN107038706A (zh) * | 2017-05-16 | 2017-08-11 | 西安电子科技大学 | 基于自适应网格的红外图像置信度评估装置及方法 |
CN108992193B (zh) * | 2017-06-06 | 2020-12-15 | 苏州笛卡测试技术有限公司 | 一种牙齿修复辅助设计方法 |
CN108992193A (zh) * | 2017-06-06 | 2018-12-14 | 苏州笛卡测试技术有限公司 | 一种牙齿修复辅助设计方法 |
CN107087114A (zh) * | 2017-06-29 | 2017-08-22 | 深圳天珑无线科技有限公司 | 一种拍摄的方法及装置 |
CN109427059A (zh) * | 2017-08-18 | 2019-03-05 | 西门子医疗有限公司 | 解剖结构的平面可视化 |
CN109427059B (zh) * | 2017-08-18 | 2022-05-24 | 西门子医疗有限公司 | 解剖结构的平面可视化 |
CN108875780A (zh) * | 2018-05-07 | 2018-11-23 | 广东省电信规划设计院有限公司 | 基于图像数据的图像间差异对象的获取方法以及装置 |
CN109584201A (zh) * | 2018-09-14 | 2019-04-05 | 新影智能科技(昆山)有限公司 | 医学图像配准方法、配准系统、存储介质、及电子设备 |
CN110197713A (zh) * | 2019-05-10 | 2019-09-03 | 上海依智医疗技术有限公司 | 一种医疗影像的处理方法、装置、设备和介质 |
CN111544020A (zh) * | 2020-04-17 | 2020-08-18 | 北京东软医疗设备有限公司 | X射线成像设备的几何校正方法及装置 |
CN112037186A (zh) * | 2020-08-24 | 2020-12-04 | 杭州深睿博联科技有限公司 | 一种基于多视图模型融合的冠脉血管提取方法及装置 |
Also Published As
Publication number | Publication date |
---|---|
CN104240287B (zh) | 2017-10-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104240287B (zh) | 一种利用ct图像生成冠脉全景图的方法及系统 | |
CN102525662B (zh) | 组织器官三维可视化手术导航系统 | |
CN101790748B (zh) | 对3d数字医学图像中的解剖实体进行分割的方法 | |
Kaus et al. | Assessment of a model-based deformable image registration approach for radiation therapy planning | |
EP2369553B1 (en) | Three-dimensional template transformation method and apparatus | |
CN103229210B (zh) | 图像配准装置 | |
US20070109299A1 (en) | Surface-based characteristic path generation | |
CN107330888A (zh) | 基于cta图像的动态心脏各腔室分割方法 | |
US7570802B2 (en) | Automated centerline detection algorithm for colon-like 3D surfaces | |
CN101933045A (zh) | 短轴后期增强心脏mri的自动3d分割 | |
US8929636B2 (en) | Method and system for image segmentation | |
CN108694007B (zh) | 从磁共振图像展开肋骨 | |
US20150110377A1 (en) | Method and system for image segmentation | |
CN106934810A (zh) | 一种脊椎矫正装置 | |
CN106548476A (zh) | 利用医学图像统计肺部三维特征形状方法 | |
CN109003283A (zh) | 一种基于主动形状模型的主动脉轮廓分割算法 | |
CN111369662A (zh) | Ct图像中血管的三维模型重建方法及系统 | |
CN107516314A (zh) | 医学图像超体素分割方法和装置 | |
CN103903255B (zh) | 一种超声图像分割方法和系统 | |
CN111127488A (zh) | 一种基于统计形状模型自动构建患者解剖结构模型的方法 | |
CN114549553A (zh) | 角度测量方法、装置、计算机设备及可读存储介质 | |
Markel et al. | A 4D biomechanical lung phantom for joint segmentation/registration evaluation | |
CN110689080B (zh) | 一种血管结构影像的平面图谱构建方法 | |
CN107492105A (zh) | 一种基于多统计信息的变分分割方法 | |
Boutchko et al. | Practical implementation of tetrahedral mesh reconstruction in emission tomography |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20171020 Termination date: 20210608 |
|
CF01 | Termination of patent right due to non-payment of annual fee |