专利名称：一种旋转x射线造影图像迭代重建方法
技术领域：
本发明属于计算机技术领域，涉及ー种旋转X射线造影图像迭代重建方法。
背景技术：
旋转X射线冠状动脉造影成像技术是继双平面血管造影成像以来又ー项得到广泛关注的冠状动脉成像技木。通过该技木，医师能够从多个角度提供完整，准确地观察血管的形状以及运动方式，对脑血管肿瘤，冠心病的诊疗具有重要意义。但仅凭旋转图像序列，难以对冠状动脉的具体结构有一个直观准确的认识，因此血管造影的三维重建成为目前学术领域的研究热点，也是医疗器械厂商亟待解决的重要问题。旋转X射线造影技术的成像几何同三维重建同锥束CT类似。由于旋转造影成像过程中伴随着心脏运动，为了获得准确的血管三维结构，需要采用迭代重建算法。但是由于该成像系统的空间分辨率高于一般的 锥束CT，该重建问题的计算量远远大于锥束CT图象重建。迭代算法的运算主要集中在投影、反投影两个环节，投影、反投影的快速计算问题成为目前旋转X射线血管造影成像三维重建问题的主要瓶颈，投影矩阵的使用可以大大加速投影，反投影操作的运算速度，但受到硬件条件限制，目前的计算机内存存储能力不足以存储旋转造影系统投影矩阵。
发明内容
技术问题本发明提供了一种对硬件系统要求较低，能够大幅度提高投影、反投影速度，便于并性计算优化，可快速准确获得血管三维结构的旋转X射线造影图像迭代重建方法。技术方案本发明的旋转X射线造影图像迭代重建方法，包括以下步骤I)从旋转X射线造影设备读取扫描数据文件，保存投影序列图像并记录如下參数射线源到探測板距离SDD，射线源到C臂旋转中心距离S0D，各投影采样的旋转角度，投影图像像素边长h，投影图像长度U个像素単位和宽度V个像素単位，所述投影图像尺寸与ニ维投影空间尺寸一致；确定ニ维投影空间坐标轴U、V，所述坐标轴U、V分别平行于ニ维投影空间的长、宽方向，根据用户的精度需求，分别设定三维图像空间长度为X个体素単位、宽度为Y个体素単位、高度为Z个体素単位，三维图像空间体素边长1，确定三维图像空间坐标轴x、y、z，所述坐标轴x、y、z分别平行于三维图像空间的长、宽、高方向，且均通过三维图像空间的中心位置；2)根据步骤I)中设定的三维图像空间长宽高对三维图像空间进行降采样操作，具体方法为选定降采样倍数dif，所述降采样倍数dif为能被三维图像空间长度X、宽度Y、高度Z分别整除的整数，将三维图像空间中由dif*dif*diff体素组成的、边长为dif*l的立方体内的体素，归并为降采样体素，对整个三维图像空间完成归并处理后，按照归并顺序排列所述降采样体素，得到降采样三维图像空间，降采样体素边长为三维图像空间体素边长的dif倍;根据步骤I)中记录的投影图像长宽，对ニ维投影空间进行降采样操作，具体方法为选定降采样倍数dpf，所述降采样倍数为能被ニ维投影空间长度U和宽度V分别整除的整数，将ニ维投影空间中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影空间降采样像素，对整个ニ维投影空间完成归并处理后，按照归并顺序排列所述投影空间降采样像素，得到降采样ニ维投影空间，投影空间降采样像素边长为ニ维投影空间像素边长的Clpf倍；根据步骤I)中记录的投影图像长宽，对步骤I)获得的投影序列图像进行降采样操作，得到降采样序列投影图像，具体方法为对投影序列中每个投影图像，将其中由cU*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影图像降采样像素，并且对所述正方形内的像素值进行累加求和，作为投影图像降采样像素值，对整个投影图像完成归并处理后，按照归并顺序排列所述投影图像降采样像素，得到降采样投影图像，投影图像降采样像素边长为投影图像像素边长的dpf倍；3)根据步骤I)中所记录的探測板距离SDD，射线源到C臂旋转中心距离S0D，初始采样旋转角度，针对步骤2)得到的降采样三维图像空间和降采样ニ维投影空间，利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并将得到的初始扫描方向低分 辨率体素索引投影矩阵以三元组存放方式保存于内存设备中；4)针对步骤2)提供的降采样三维图像空间，构造各旋转角度下降采样三维图像空间旋转矩阵R，具体方法为首先根据图像旋转角度位置关系，利用线性插值算法构造并存储降采样图像空间切面层旋转矩阵R°，随后根据降采样三维图像空间旋转矩阵R的分块对
称性，利用公式R =得到降采样图像空间旋转矩阵R ；5)对步骤2)得到的降采样投影序列图像采用顶帽滤波方法进行滤波处理，对滤波处理结果ニ值化得到第一分割結果，对步骤2)得到的降采样投影序列图像采用Frangi血管滤波方法进行滤波处理，对滤波处理结果ニ值化得到第二分割結果，对第一分割结果和第二分割结果求并集，得到低分辨率投影序列图像血管分割結果；6)利用步骤3)所得到的三元组格式存储的初始扫描方向低分辨率体素索引投影矩阵和步骤4)所得到的各旋转角度降采样图像空间旋转矩阵，对步骤5)所得到的低分辨率投影序列图像血管分割结果进行反投影操作具体方法为对低分辨率投影图像血管分割结果，根据像素编号排列拉伸低分辨率投影图像血管分割结果，得到低分辨率投影图像血管分割结果向量，将所得向量同初始扫描方向低分辨率体素索引投影矩阵的转置矩阵相乘得到中间低分辨率反投影結果；将该中间低分辨率反投影结果同对应旋转角度降采样图像空间旋转矩阵的转置矩阵相乘，得到该旋转角度下的低分辨率反投影结果；7)利用步骤6)所得到的低分辨率反投影结果确定低分辨率三维图像空间血管掩摸，具体方法为当投影数目小于等于5吋，对各旋转角度下低分辨率反投影结果求交集，作为低分辨率三维图像空间血管掩模；当投影数目大于5时，利用对各旋转角度下低分辨率反投影结果求和并进行阈值划分，作为低分辨率三维图像空间血管掩模；8)对步骤7)得到的低分辨率三维图像空间血管掩模，根据步骤2)所设定的降采样倍数dif进行升采样操作，得到三维图像空间血管掩模，具体方法为将低分辨率三维图像空间血管掩模的每ー个低分辨率体素在x,y,z方向等分为dif份，拆分成为dif*dif*dif个高分辨率体素，将该低分辨率体素的值作为其拆分成的高分辨率体素的值，并按照拆分顺序排列高分辨率体素，当低分辨率三维图像空间血管掩模的所有低分辨率体素完成上述操作后，即得到三维图像空间血管掩模；9)根据步骤I)中所记录的射线源到探測板距离SDD，射线源到C臂旋转中心距离S0D，初始采样旋转角度，投影图像长度U、宽度V，投影图像像素边长h，设定的三维图像空间长度X、宽度Y、高度Z，三维图像空间体素边长1，以及步骤8)所得的到三维图像空间血管掩模，利用距离驱动算法计算各旋转角度下的完整投影矩阵；10)根据步骤I)所记录的投影序列图像，根据步骤9)所得到的各旋转角度下的完整投影矩阵，完成三维血管结构重建，具体方法为101)设定ー个重建结果向量，所述重建结果向量长度为图像空间体素总个数XXYXZ，重建结果向量元素值全部为O;设定ー个步长向量，步长向量长度为图像空间体素总个数XXYXZ，步长向量元素值全部为O ;设定ー个单元向量，単元向量长度为图像空间体素总个数XXYXZ，単元向量元素值全部为I ; 102)对于每个投影角度，进行如下操作将单元向量与当前投影角度下的完整投影矩阵相乘，将相乘结果同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前投影角度的单方向步长向量，将该单方向步长向量累加至步长向量；对所有投影角度完成上述操作后，将累加完成的步长向量记录保存为迭代步长向量；103)执行重建结果向量迭代更新步骤100至300次，完成指定次数的重建结果向量迭代更新步骤后，将重建结果向量保存为重建结果三维血管体数据，所述重建结果向量迭代更新的方法如下对于每个投影角度，进行如下操作将当前的重建结果向量同当前投影角度下的完整投影矩阵相乘，得到当前方向临时投影向量，用当前方向投影图像向量减去所述当前方向临时投影向量，得到当前方向临时投影差向量，将所述当前方向临时投影差向量同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前方向差异反投影向量，用所述当前方向差异反投影向量点除以迭代步长向量，得到当前方向迭代步进向量，将当前方向迭代步进向量累加至重建结果向量；对所有角度完成上述操作后，完成一次重建结果向量更新。本发明中，步骤3)中利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并将得到的初始扫描方向低分辨率体素索引投影矩阵以三元组存放方式保存于内存设备中的具体方法为31)设定发射源的空间坐标为(_sSQD，0，0)，设定探測板中心点坐标为(sSDD-sSQD，0，O), Ssod表示射线源到C臂旋转中心距离，Ssdd表示射线源到探測板距离；32)对所有降采样体素作如下操作从步骤31)设定的发射源向序号为jd的体素各角点引出射线，各射线将与降采样ニ维投影空间相交，其八个交点的连线所包围区域面积记<_，所述包围区域同序号为id的像素相交面积为Sd' μ，则au=sd' JiZsdj-即为低分辨率投影矩阵第id行，第jd列的元素值；33)将所得到的所有具有非零值的元素的行号、列号以及元素值保存，得到三元组格式存储的初始扫描方向低分辨率体素索引投影矩阵。本发明的步骤4)中构造并存储降采样图像空间切面层旋转矩阵R°的具体方法为记投影采样的旋转角度为α，将三维图像空间按照ζ方向体素単位分组，共计Z个切面层，选取任一切面层，将其中的每个体素围绕Z轴旋转-a角度，记录旋转后的位置的邻域体素编号Itl, I1, 12，I3，按照线性插值方法确定矩阵第k行，第Itl列，第I1列，第I2列，第I3列的四个元素的元素值，其中k为体素的序号，将所得到的所有具有非零值的元素的行号，列号以及元素值保存，得到三元组格式存储的降采样图像空间切面层旋转矩阵R°。本发明中，步骤9)中利用距离驱动算法计算各旋转角度下的完整投影矩阵的具体过程为91)设定发射源的空间坐标为(-SsqdXcoS a，-sS0DXsina , O)，设定探测板中心点坐丰不为((Ssdd-Ssod) X cos Ct，(sSDD—sS0D) X sin ct，O);92)对所有三维图像空间体素作如下操作如果在步骤8所得到的三维图像空间血管掩模中，序号为j的体素没有非零值，则不做任何操作直接进入下一个体素的操作，如果具有非零值，则根据步骤91)设定的发射源向该体素各角点引出射线，各射线将与降采样 ニ维投影空间相交，其八个交点的连线所包围区域面积记为ち，同序号为i的像素相交面积为t, ノも即为投影矩阵第i行，第j列的元素值；93)将所得到的所有具有非零值的元素的行号，列号以及元素值保存，得到三元组格式存储的α角度方向投影矩阵。本发明是一种旋转X射线造影图像迭代重建方法，该方法分为3个阶段第一阶段构造低分辨率投影矩阵，并将完整矩阵拆解为单一角度矩阵和旋转矩阵2个分量进行简化存储。在该阶段首先在単一投影方向上，利用体素索引投影算法构造低分辨率投影矩阵；之后构造3D体数据旋转矩阵，利用该旋转矩阵自身対称性进行简化存储。第二阶段在第一阶段得到低分辨率投影矩阵的基础上，进行进一步基于投影内容的简化。在该阶段，首先对采集投影数据进行降采样及快速分割；之后利用低分辨率投影矩阵进行反投影，确定3D热区，其后对该区域进行升采样。最后利用升采样区域构造投影矩阵掩摸，结合系统參数计算精确投影矩阵。第三阶段在第二阶段得到的精确投影矩阵的基础上进行旋转X射线造影系统的三维血管重建。在该阶段，首先利用単元向量计算步长向量；之后通过一定次数的迭代步骤，反复计算修正当前结果向量同投影数据之间的差异，井根据步长向量设定每次的具体修正量，在迭代步骤完成后获得三维血管重建数据。有益效果本发明同现有技术方法相比，具有如下优点本方法确定的投影矩阵构造方法大幅降低了矩阵存储空间，一方面，利用稀疏矩阵存储技术，可以降低所需存储空间，另一方面，同一般稀疏矩阵存储方式相比，考虑了血管掩模后，所需存储空间可进一歩大幅降低。以如下数据为例ニ维投影大小为512*512，三维图像大小为256*256*256 :单ー投影角度完整投影矩阵大小为512*512*256*256*256，以每个元素值占用4字节存储空间为例，共需要存储空间16384GB ;如果使用稀疏矩阵存储，由于平均每个体素对应相关像素为12个左右，则矩阵非零元素个数为256*256*256*9，每个非零元素的存储需要耗费4字节元素值，以及4字节行列号，共需要存储空间约为I. 7GB ;利用本发明提出的算法，在考虑图像掩模的情况下，通常血管掩模的非零元素个数在2*105左右，每个非零元素对应12个相关像素，共需要存储空间约为20Mb字节。同一般稀疏矩阵相比，本方法所需存储空间降至1%左右，利用获得的投影矩阵同投影向量相乘，则完成投影操作，利用获得的投影矩阵的转置，同图像向量相乘，则完成反投影操作。同一般稀疏矩阵相比本方法确定的投影矩阵构造方法大幅提高了投影、反投影计算速度，记录通过矩阵相乘运算完成投影、反投影操作所需时间，理论上，矩阵尺寸减小至1%，则速度提高应为100倍，实际应用当中，由于数据传递，数据交互等底层原因，实验证明速度提升约为20-40倍。获得简化后的投影矩阵后，可以利用其完成旋转X射线造影系统的三维血管重建，得到准确的三维血管结构。
图I是本发明方法的旋转造影系统投影生成示意图。图2a是本发明方法的距离驱动像素索引投影矩阵几何示意图。图2b是本发明方法的距离驱动像素索引投影矩阵元素权重计算示意图。 图3是本发明方法的ニ值化投影分割反投影叠加示意图。图4是本发明方法的整体流程图。图5是本发明方法步骤3)的子流程图。图6是本发明方法步骤4)中构造并存储降采样图像空间切面层旋转矩阵R°的子流程图。图7是本发明方法步骤9)中利用距离驱动算法计算各旋转角度下的完整投影矩阵的子流程图。
具体实施例方式本发明的旋转X射线造影图像迭代重建方法，包括下列步骤I)从旋转X射线造影设备读取扫描数据文件，保存投影序列图像并记录如下參数射线源到探測板距离SDD，射线源到C臂旋转中心距离S0D，各投影采样的旋转角度，投影图像像素边长h，投影图像像长度U个像素単位和宽度V个像素単位，投影图像尺寸与ニ维投影空间像素尺寸一致，每个象素单位为长度为h的正方形，确定ニ维投影空间坐标轴U，V，坐标轴U，V分别平行于ニ维投影空间的长，宽方向，根据用户的精度需求，分别设定三维图像空间长度为X个体素単位、宽度为Y个体素単位、高度为Z个体素単位，每个体素单位为边长I的立方体，确定三维图像空间坐标轴X，y, z，坐标轴X，y, z分别平行于三维图像空间的长、宽、高方向，且均通过三维图像空间的中心位置，也即坐标原点为三维图像空间的中心位置。系统模拟图如图I所示，图中S为发射源位置，O为旋转中心，D为探测板中心位置。射线源到探測板距离=SD，射线源到C臂旋转中心距离=S0，记图像空间体素尺寸IX 1X1，图像空间包含体素个数XXYXZ ;记投影图像像素尺寸为hXh，投影空间包含像素个数UX V;记向量X/)和X轴夹角为a，每次以O为中心旋转角度为Λ a，总投影角度数为K。2)根据步骤I)中设定的三维图像空间长宽高对三维图像空间进行降采样操作，具体方法为选定降采样倍数dif，所述降采样倍数dif为能被三维图像空间长度X、宽度Y、高度Z分别整除的整数，通常可以选择2，3，4，将三维图像空间中由dif*dif*dif个体素组成的、边长为dif*l的立方体内的体素，归并为降采样体素，对整个三维图像空间完成归并处理后，按照归并顺序排列所述降采样体素，得到降采样三维图像空间，降采样体素边长为三维图像空间体素边长的dif倍，记Xif = X/dif, Yif = Y/dif, Zif = Z/dif则降采样图像空间体素尺寸IdifX IdifXldif,降采样图像空间包含体素个数为Xif*Yif*Zif ；根据步骤I)中记录的投影图像长宽，对ニ维投影空间进行降采样操作，具体方法为选定降采样倍数dpf，所述降采样倍数为能被ニ维投影空间长度U和宽度V分别整除的整数，通常可以选择2，3 ，4，将ニ维投影空间中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影空间降采样像素，对整个ニ维投影空间完成归并处理后，按照归并顺序排列所述投影空间降采样像素，得到降采样ニ维投影空间，投影空间降采样像素边长为ニ维投影空间像素边长的Clpf倍，记Upf = U/dpf, Vpf = V/dpf则降采样投影空间像素尺寸hdpf Xhdpf,降采样投影空间包含像素个数UpfXVpf ；根据步骤I)中记录的投影图像长宽，对步骤I)获得的投影序列图像进行降采样操作，得到降采样序列投影图像，具体方法为对投影序列中每个投影图像，将其中由cU*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影图像降采样像素，并且对所述正方形内的像素值进行累加求和，作为投影图像降采样像素值，对整个投影图像完成归并处理后，按照归并顺序排列所述投影图像降采样像素，得到降采样投影图像，投影图像降采样像素边长为投影图像像素边长的Clpf倍，降采样投影图像像素尺寸hdpf Xhdpf,降采样投影图像包含像素个数Upf X Vpf。3)根据步骤I)中所记录的探測板距离SDD，射线源到C臂旋转中心距离S0D，初始采样旋转角度，针对步骤2)得到的降采样三维图像空间和降采样ニ维投影空间，利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，具体方法为取初始方向a=0，设定发射源的空间坐标为(_sS(I)，0，0)，探測板中心坐标为(sSDD-sS(I)，0，0)，按照X，y, z的顺序遍历降采样图像空间体素。如图2 (a)所示，对于序号为(VX，Vy，VZ)的体素，记其编号jd = vzXXifXYif+vyXXif+vx,作点源S和编号jd体素中心点连线,其延长线与探测板相交点在探測平面的坐标为(pud，pvd)，作点源S和编号jd体素各角点连线，其延长线与探測板相交点所包围矩形面积记为cC. ,.，其面积为<_，若序号为(m，η)的投影像素点，其编号记为id = pvdXU+pud,与,相交，且相交面积为sd/ u，，如图2 (b)所示,则记低分辨率体素索引投影矩阵矩阵第id行，jd列的元素值~=パ^./Sdj ；将所有具有非零值的元素的行号，列号以及元素值以三元组存放方式保存于内存设备中，得到的初始扫描方向低分辨率体素索引投影矩阵；4)针对步骤2)提供的降采样三维图像空间，构造各旋转角度下降采样三维图像空间旋转矩阵R具体方法为首先对图像空间沿z轴(旋转轴)进行划分，根据图像旋转角度位置关系，利用线性插值算法构造降采样图像空间切面层旋转矩阵R°，该矩阵大小为，XifYifXXifYif为稀疏矩阵，对于所选切面层，可看作ニ维图像处理，以图像中心为原点，对于
像素(row，col)，其2D坐标分别为(a°，b°)，其中％ ==し"/-记其序
Z ο
号为j = col XXif+row。对应稀疏旋转矩阵R°的第j列。根据2D旋转坐标变换公式，经过旋转_a角度后该像素点的坐标位置为
权利要求
1.一种旋转X射线造影图像迭代重建方法，其特征在于，包括以下步骤 1)从旋转X射线造影设备读取扫描数据文件，保存投影序列图像并记录如下参数射线源到探测板距离SDD，射线源到C臂旋转中心距离SOD，各投影采样的旋转角度，投影图像像素边长h，投影图像长度U个像素单位和宽度V个像素单位，所述投影图像尺寸与二维投影空间尺寸一致；确定二维投影空间坐标轴U、V，所述坐标轴U、V分别平行于二维投影空间的长、宽方向，根据用户的精度需求，分别设定三维图像空间长度为X个体素单位、宽度为Y个体素单位、高度为Z个体素单位，三维图像空间体素边长1，确定三维图像空间坐标轴χ、y、z，所述坐标轴x、y、z分别平行于三维图像空间的长、宽、高方向，且均通过三维图像空间的中心位置； 2)根据步骤I)中设定的三维图像空间长宽高对三维图像空间进行降采样操作，具体方法为选定降采样倍数dif，所述降采样倍数dif为能被三维图像空间长度X、宽度Y、高度Z分 别整除的整数，将三维图像空间中由dif*dif*dif个体素组成的、边长为dif*l的立方体内的体素，归并为降采样体素，对整个三维图像空间完成归并处理后,按照归并顺序排列所述降采样体素，得到降采样三维图像空间，降采样体素边长为三维图像空间体素边长的dif倍； 根据步骤I)中记录的投影图像长宽，对二维投影空间进行降采样操作，具体方法为选定降采样倍数Clpf，所述降采样倍数为能被二维投影空间长度U和宽度V分别整除的整数，将二维投影空间中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影空间降采样像素，对整个二维投影空间完成归并处理后，按照归并顺序排列所述投影空间降采样像素，得到降采样二维投影空间，投影空间降采样像素边长为二维投影空间像素边长的Clpf倍； 根据步骤I)中记录的投影图像长宽，对步骤I)获得的投影序列图像进行降采样操作，得到降采样序列投影图像，具体方法为对投影序列中每个投影图像，将其中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影图像降采样像素，并且对所述正方形内的像素值进行累加求和，作为投影图像降采样像素值，对整个投影图像完成归并处理后，按照归并顺序排列所述投影图像降采样像素，得到降采样投影图像，投影图像降采样像素边长为投影图像像素边长的Clpf倍； 3)根据步骤I)中所记录的探测板距离SDD，射线源到C臂旋转中心距离SOD，初始采样旋转角度，针对步骤2)得到的降采样三维图像空间和降采样二维投影空间，利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并将得到的初始扫描方向低分辨率体素索引投影矩阵以三元组存放方式保存于内存设备中； 4)针对步骤2)提供的降采样三维图像空间，构造各旋转角度下降采样三维图像空间旋转矩阵R，具体方法为首先根据图像旋转角度位置关系，利用线性插值算法构造并存储降采样图像空间切面层旋转矩阵R°，随后根据降采样三维图像空间旋转矩阵R的分块对称性，利用公式Λ =到降采样图像空间旋转矩阵R ； Z 5)对步骤2)得到的降采样投影序列图像采用顶帽滤波方法进行滤波处理，对滤波处理结果二值化得到第一分割结果，对步骤2)得到的降采样投影序列图像采用Frangi血管滤波方法进行滤波处理，对滤波处理结果二值化得到第二分割结果，对第一分割结果和第二分割结果求并集，得到低分辨率投影序列图像血管分割结果；6)利用步骤3)所得到的三元组格式存储的初始扫描方向低分辨率体素索引投影矩阵和步骤4)所得到的各旋转角度降采样图像空间旋转矩阵，对步骤5)所得到的低分辨率投影序列图像血管分割结果进行反投影操作具体方法为对低分辨率投影图像血管分割结果，根据像素编号排列拉伸低分辨率投影图像血管分割结果，得到低分辨率投影图像血管分割结果向量，将所得向量同初始扫描方向低分辨率体素索引投影矩阵的转置矩阵相乘得到中间低分辨率反投影结果；将该中间低分辨率反投影结果同对应旋转角度降采样图像空间旋转矩阵的转置矩阵相乘，得到该旋转角度下的低分辨率反投影结果； 7)利用步骤6)所得到的低分辨率反投影结果确定低分辨率三维图像空间血管掩模，具体方法为当投影数目小于等于5时，对各旋转角度下低分辨率反投影结果求交集，作为低分辨率三维图像空间血管掩模；当投影数目大于5时，利用对各旋转角度下低分辨率反投影结果求和并进行阈值划分，作为低分辨率三维图像空间血管掩模； 8)对步骤7)得到的低分辨率三维图像空间血管掩模，根据步骤2)所设定的降采样倍数dif进行升采样操作，得到三维图像空间血管掩模，具体方法为将低分辨率三维图像空 间血管掩模的每一个低分辨率体素在X，1，Z方向等分为dif份，拆分成为dif*dif*dif个高分辨率体素，将该低分辨率体素的值作为其拆分成的高分辨率体素的值，并按照拆分顺序排列高分辨率体素，当低分辨率三维图像空间血管掩模的所有低分辨率体素完成上述操作后，即得到三维图像空间血管掩模； 9)根据步骤I)中所记录的射线源到探测板距离SDD，射线源到C臂旋转中心距离SOD，初始采样旋转角度，投影图像长度U、宽度V，投影图像像素边长h，设定的三维图像空间长度X、宽度Y、高度Z，三维图像空间体素边长1，以及步骤8)所得的到三维图像空间血管掩模，利用距离驱动算法计算各旋转角度下的完整投影矩阵； 10)根据步骤I)所记录的投影序列图像，根据步骤9)所得到的各旋转角度下的完整投影矩阵，完成三维血管结构重建，具体方法为 101)设定一个重建结果向量，所述重建结果向量长度为图像空间体素总个数XXYXZ，重建结果向量元素值全部为O ;设定一个步长向量，所述步长向量长度为图像空间体素总个数XXYXZ，步长向量元素值全部为O ;设定一个单元向量，所述单元向量长度为图像空间体素总个数XXYXZ，单元向量元素值全部为I ; 102)对于每个投影角度，进行如下操作将单元向量与当前投影角度下的完整投影矩阵相乘，将相乘结果同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前投影角度的单方向步长向量，将该单方向步长向量累加至步长向量；对所有投影角度完成上述操作后，将累加完成的步长向量记录保存为迭代步长向量； 103)执行重建结果向量迭代更新步骤100至300次，完成指定次数的重建结果向量迭代更新步骤后，将重建结果向量保存为重建结果三维血管体数据，所述重建结果向量迭代更新的方法如下 对于每个投影角度，进行如下操作将当前的重建结果向量同当前投影角度下的完整投影矩阵相乘，得到当前方向临时投影向量，用当前方向投影图像向量减去所述当前方向临时投影向量，得到当前方向临时投影差向量，将所述当前方向临时投影差向量同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前方向差异反投影向量，用所述当前方向差异反投影向量点除以迭代步长向量，得到当前方向迭代步进向量，将当前方向迭代步进向量累加至重建结果向量；对所有角度完成上述操作后，完成一次重建结果向量更新。
2.根据权利要求I所述的旋转X射线造影图像迭代重建方法，其特征在于，所述步骤3)中利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并将得到的初始扫描方向低分辨率体素索引投影矩阵以三元组存放方式保存于内存设备中的具体方法为 31)设定发射源的空间坐标为(_sSQD，0，0)，设定探测板中心点坐标为(sSDD-sSQD，0，0)，Ssod表示射线源到C臂旋转中心距离，Ssdd表示射线源到探测板距离； 32)对所有降采样体素作如下操作从步骤31)设定的发射源向序号为jd的体素各角点引出射线，各射线将与降采样二维投影空间相交，其八个交点的连线所包围区域面积记<，所述包围区域同序号为id的像素相交面积为Sd' μ则au=sd' JiZsdj-即为低分辨率投影矩阵第0行，第jd列的元素值；· 33)将所得到的所有具有非零值的元素的行号、列号以及元素值保存，得到三元组格式存储的初始扫描方向低分辨率体素索引投影矩阵。
3.根据权利要求I所述的旋转X射线造影图像迭代重建方法，其特征在于，所述的步骤4)中构造并存储降采样图像空间切面层旋转矩阵R°的具体方法为记投影采样的旋转角度为α，将三维图像空间按照ζ方向体素单位分组，共计Z个切面层，选取任一切面层，将其中的每个体素围绕ζ轴旋转-α角度，记录旋转后的位置的邻域体素编号Itl, I1, I2, I3，按照线性插值方法确定矩阵第k行，第Itl列，第I1列，第I2列，第I3列的四个元素的元素值，其中k为体素的序号，将所得到的所有具有非零值的元素的行号，列号以及元素值保存，得到三元组格式存储的降采样图像空间切面层旋转矩阵R°。
4.根据权利要求I所述的旋转X射线造影图像迭代重建方法，其特征在于，所述步骤9)中利用距离驱动算法计算各旋转角度下的完整投影矩阵的具体过程为 91)设定发射源的空间坐标为(-SsodXc0Sa，-sS0DXsina, O)，设定探测板中心点坐标为((sSDD-sS0D) Xcos a , (sSDD-sS0D) Xsina , O)； 92)对所有三维图像空间体素作如下操作如果在步骤8所得到的三维图像空间血管掩模中，序号为j的体素没有非零值，则不做任何操作直接进入下一个体素的操作，如果具有非零值，则根据步骤91)设定的发射源向该体素各角点引出射线，各射线将与降采样二维投影空间相交，其八个交点的连线所包围区域面积记为\_，同序号为i的像素相交面积为t/ ji，则au=t' /%即为投影矩阵第i行，第j列的元素值； 93)将所得到的所有具有非零值的元素的行号，列号以及元素值保存，得到三元组格式存储的α角度方向投影矩阵。
全文摘要
本发明公开了一种旋转X射线造影图像迭代重建方法，首先在第一阶段构造低分辨率投影矩阵，并将完整矩阵拆解为单一角度矩阵和旋转矩阵2个分量进行简化存储，然后第二阶段在第一阶段得到低分辨率投影矩阵的基础上，进行进一步基于投影内容的简化，最后在第三阶段进行三维血管重建本方法采用的经过掩模简化的投影矩阵解决了旋转X射线造影系统投影，反投影计算量过大，计算时间过长的问题，能够有效获得三维血管结构，帮助临床医师进行诊断。
文档编号G06T11/00GK102842141SQ20121022816
公开日2012年12月26日 申请日期2012年7月3日 优先权日2012年7月3日
发明者胡轶宁, 谢理哲, 沈傲东, 罗立民 申请人:东南大学