一种旋转X射线造影图像迭代重建方法.pdf

本发明公开了一种旋转X射线造影图像迭代重建方法，首先在第一阶段构造低分辨率投影矩阵，并将完整矩阵拆解为单一角度矩阵和旋转矩阵2个分量进行简化存储，然后第二阶段在第一阶段得到低分辨率投影矩阵的基础上，进行进一步基于投影内容的简化，最后在第三阶段进行三维血管重建本方法采用的经过掩模简化的投影矩阵解决了旋转X射线造影系统投影，反投影计算量过大，计算时间过长的问题，能够有效获得三维血管结构，帮助临床医师进行诊断。

1.一种旋转X射线造影图像迭代重建方法，其特征在于，包括以下步骤：
1）从旋转X射线造影设备读取扫描数据文件，保存投影序列图像并记录如下参
数：射线源到探测板距离SDD，射线源到C臂旋转中心距离SOD，各投影采样的旋
转角度，投影图像像素边长h，投影图像长度U个像素单位和宽度V个像素单位，所
述投影图像尺寸与二维投影空间尺寸一致；确定二维投影空间坐标轴u、v，所述坐标
轴u、v分别平行于二维投影空间的长、宽方向，根据用户的精度需求，分别设定三
维图像空间长度为X个体素单位、宽度为Y个体素单位、高度为Z个体素单位，三维
图像空间体素边长l，确定三维图像空间坐标轴x、y、z，所述坐标轴x、y、z分别平
行于三维图像空间的长、宽、高方向，且均通过三维图像空间的中心位置；
2）根据步骤1）中设定的三维图像空间长宽高对三维图像空间进行降采样操作，
具体方法为：选定降采样倍数dif，所述降采样倍数dif为能被三维图像空间长度X、
宽度Y、高度Z分别整除的整数，将三维图像空间中由dif*dif*dif个体素组成的、边长
为dif*l的立方体内的体素，归并为降采样体素，对整个三维图像空间完成归并处理后，
按照归并顺序排列所述降采样体素，得到降采样三维图像空间，降采样体素边长为三
维图像空间体素边长的dif倍；
根据步骤1）中记录的投影图像长宽，对二维投影空间进行降采样操作，具体方
法为：选定降采样倍数dpf，所述降采样倍数为能被二维投影空间长度U和宽度V分
别整除的整数，将二维投影空间中由dpf*dpf个像素组成的、边长为dpf*h的正方形内
的像素，归并为投影空间降采样像素，对整个二维投影空间完成归并处理后，按照归
并顺序排列所述投影空间降采样像素，得到降采样二维投影空间，投影空间降采样像
素边长为二维投影空间像素边长的dpf倍；
根据步骤1）中记录的投影图像长宽，对步骤1）获得的投影序列图像进行降采
样操作，得到降采样序列投影图像，具体方法为：对投影序列中每个投影图像，将其
中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影图像降采样
像素，并且对所述正方形内的像素值进行累加求和，作为投影图像降采样像素值，对
整个投影图像完成归并处理后，按照归并顺序排列所述投影图像降采样像素，得到降
采样投影图像，投影图像降采样像素边长为投影图像像素边长的dpf倍；
3）根据步骤1）中所记录的探测板距离SDD，射线源到C臂旋转中心距离SOD，
初始采样旋转角度，针对步骤2）得到的降采样三维图像空间和降采样二维投影空间，
利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并将得到的初始
扫描方向低分辨率体素索引投影矩阵以三元组存放方式保存于内存设备中；
4）针对步骤2）提供的降采样三维图像空间，构造各旋转角度下降采样三维图像
空间旋转矩阵R，具体方法为：首先根据图像旋转角度位置关系，利用线性插值算法
构造并存储降采样图像空间切面层旋转矩阵R0，随后根据降采样三维图像空间旋转矩
阵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）根据步骤1）中所记录的射线源到探测板距离SDD，射线源到C臂旋转中心
距离SOD，初始采样旋转角度，投影图像长度U、宽度V，投影图像像素边长h，设
定的三维图像空间长度X、宽度Y、高度Z，三维图像空间体素边长l，以及步骤8）
所得的到三维图像空间血管掩模，利用距离驱动算法计算各旋转角度下的完整投影矩
阵；
10）根据步骤1）所记录的投影序列图像，根据步骤9）所得到的各旋转角度下
的完整投影矩阵，完成三维血管结构重建，具体方法为：
101）设定一个重建结果向量，所述重建结果向量长度为图像空间体素总个数
X×Y×Z，重建结果向量元素值全部为0；设定一个步长向量，所述步长向量长度为
图像空间体素总个数X×Y×Z，步长向量元素值全部为0；设定一个单元向量，所述
单元向量长度为图像空间体素总个数X×Y×Z，单元向量元素值全部为1；
102）对于每个投影角度，进行如下操作：将单元向量与当前投影角度下的完整
投影矩阵相乘，将相乘结果同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到
当前投影角度的单方向步长向量，将该单方向步长向量累加至步长向量；对所有投影
角度完成上述操作后，将累加完成的步长向量记录保存为迭代步长向量；
103）执行重建结果向量迭代更新步骤100至300次，完成指定次数的重建结果
向量迭代更新步骤后，将重建结果向量保存为重建结果三维血管体数据，所述重建结
果向量迭代更新的方法如下：
对于每个投影角度，进行如下操作：将当前的重建结果向量同当前投影角度下的
完整投影矩阵相乘，得到当前方向临时投影向量，用当前方向投影图像向量减去所述
当前方向临时投影向量，得到当前方向临时投影差向量，将所述当前方向临时投影差
向量同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前方向差异反投影向
量，用所述当前方向差异反投影向量点除以迭代步长向量，得到当前方向迭代步进向
量，将当前方向迭代步进向量累加至重建结果向量；对所有角度完成上述操作后，完
成一次重建结果向量更新。
2.根据权利要求1所述的旋转X射线造影图像迭代重建方法，其特征在于，所
述步骤3）中利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并
将得到的初始扫描方向低分辨率体素索引投影矩阵以三元组存放方式保存于内存设
备中的具体方法为：
31）设定发射源的空间坐标为（-sSOD，0，0），设定探测板中心点坐标为（sSDD
–sSOD，0，0），sSOD表示射线源到C臂旋转中心距离，sSDD表示射线源到探测板距离；
32）对所有降采样体素作如下操作：从步骤31）设定的发射源向序号为jd的体
素各角点引出射线，各射线将与降采样二维投影空间相交，其八个交点的连线所包围
区域面积记sdj，所述包围区域同序号为id的像素相交面积为sd′ji，则aij=sd′ji/sdj即为低
分辨率投影矩阵第id行，第jd列的元素值；
33）将所得到的所有具有非零值的元素的行号、列号以及元素值保存，得到三元
组格式存储的初始扫描方向低分辨率体素索引投影矩阵。
3.根据权利要求1所述的旋转X射线造影图像迭代重建方法，其特征在于，所
述的步骤4）中构造并存储降采样图像空间切面层旋转矩阵R0的具体方法为：记投影
采样的旋转角度为α，将三维图像空间按照z方向体素单位分组，共计Z个切面层，
选取任一切面层，将其中的每个体素围绕z轴旋转-α角度，记录旋转后的位置的邻域
体素编号l0，l1，l2，l3，按照线性插值方法确定矩阵第k行，第l0列，第l1列，第l2
列，第l3列的四个元素的元素值，其中k为体素的序号，将所得到的所有具有非零值
的元素的行号，列号以及元素值保存，得到三元组格式存储的降采样图像空间切面层
旋转矩阵R0。
4.根据权利要求1所述的旋转X射线造影图像迭代重建方法，其特征在于，所
述步骤9）中利用距离驱动算法计算各旋转角度下的完整投影矩阵的具体过程为：
91）设定发射源的空间坐标为(-sSOD×cosα,-sSOD×sinα,0)，设定探测板中心点
坐标为((sSDD-sSOD)×cosα,(sSDD-sSOD)×sinα,0)；
92）对所有三维图像空间体素作如下操作：如果在步骤8所得到的三维图像空间
血管掩模中，序号为j的体素没有非零值，则不做任何操作直接进入下一个体素的操
作，如果具有非零值，则根据步骤91）设定的发射源向该体素各角点引出射线，各射
线将与降采样二维投影空间相交，其八个交点的连线所包围区域面积记为tj，同序号
为i的像素相交面积为t′ji，则aij=t′ji/tj即为投影矩阵第i行，第j列的元素值；
93）将所得到的所有具有非零值的元素的行号，列号以及元素值保存，得到三元
组格式存储的α角度方向投影矩阵。

一种旋转X射线造影图像迭代重建方法

技术领域

本发明属于计算机技术领域，涉及一种旋转X射线造影图像迭代重建方法。 

背景技术

旋转X射线冠状动脉造影成像技术是继双平面血管造影成像以来又一项得到广泛关注的冠状动脉成像技术。通过该技术，医师能够从多个角度提供完整，准确地观察血管的形状以及运动方式，对脑血管肿瘤，冠心病的诊疗具有重要意义。但仅凭旋转图像序列，难以对冠状动脉的具体结构有一个直观准确的认识，因此血管造影的三维重建成为目前学术领域的研究热点，也是医疗器械厂商亟待解决的重要问题。旋转X射线造影技术的成像几何同三维重建同锥束CT类似。由于旋转造影成像过程中伴随着心脏运动，为了获得准确的血管三维结构，需要采用迭代重建算法。但是由于该成像系统的空间分辨率高于一般的锥束CT，该重建问题的计算量远远大于锥束CT图象重建。迭代算法的运算主要集中在投影、反投影两个环节，投影、反投影的快速计算问题成为目前旋转X射线血管造影成像三维重建问题的主要瓶颈，投影矩阵的使用可以大大加速投影，反投影操作的运算速度，但受到硬件条件限制，目前的计算机内存存储能力不足以存储旋转造影系统投影矩阵。 

发明内容

技术问题：本发明提供了一种对硬件系统要求较低，能够大幅度提高投影、反投影速度，便于并性计算优化，可快速准确获得血管三维结构的旋转X射线造影图像迭代重建方法。 

技术方案：本发明的旋转X射线造影图像迭代重建方法，包括以下步骤： 

1）从旋转X射线造影设备读取扫描数据文件，保存投影序列图像并记录如下参数：射线源到探测板距离SDD，射线源到C臂旋转中心距离SOD，各投影采样的旋转角度，投影图像像素边长h，投影图像长度U个像素单位和宽度V个像素单位，所述投影图像尺寸与二维投影空间尺寸一致；确定二维投影空间坐标轴u、v，所述坐标轴u、v分别平行于二维投影空间的长、宽方向，根据用户的精度需求，分别设定三 维图像空间长度为X个体素单位、宽度为Y个体素单位、高度为Z个体素单位，三维图像空间体素边长l，确定三维图像空间坐标轴x、y、z，所述坐标轴x、y、z分别平行于三维图像空间的长、宽、高方向，且均通过三维图像空间的中心位置； 

2）根据步骤1）中设定的三维图像空间长宽高对三维图像空间进行降采样操作，具体方法为：选定降采样倍数dif，所述降采样倍数dif为能被三维图像空间长度X、宽度Y、高度Z分别整除的整数，将三维图像空间中由dif*dif*dif个体素组成的、边长为dif*l的立方体内的体素，归并为降采样体素，对整个三维图像空间完成归并处理后，按照归并顺序排列所述降采样体素，得到降采样三维图像空间，降采样体素边长为三维图像空间体素边长的dif倍； 

根据步骤1）中记录的投影图像长宽，对二维投影空间进行降采样操作，具体方法为：选定降采样倍数dpf，所述降采样倍数为能被二维投影空间长度U和宽度V分别整除的整数，将二维投影空间中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影空间降采样像素，对整个二维投影空间完成归并处理后，按照归并顺序排列所述投影空间降采样像素，得到降采样二维投影空间，投影空间降采样像素边长为二维投影空间像素边长的dpf倍； 

根据步骤1）中记录的投影图像长宽，对步骤1）获得的投影序列图像进行降采样操作，得到降采样序列投影图像，具体方法为：对投影序列中每个投影图像，将其中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影图像降采样像素，并且对所述正方形内的像素值进行累加求和，作为投影图像降采样像素值，对整个投影图像完成归并处理后，按照归并顺序排列所述投影图像降采样像素，得到降采样投影图像，投影图像降采样像素边长为投影图像像素边长的dpf倍； 

3）根据步骤1）中所记录的探测板距离SDD，射线源到C臂旋转中心距离SOD，初始采样旋转角度，针对步骤2）得到的降采样三维图像空间和降采样二维投影空间，利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并将得到的初始扫描方向低分辨率体素索引投影矩阵以三元组存放方式保存于内存设备中； 

4）针对步骤2）提供的降采样三维图像空间，构造各旋转角度下降采样三维图像空间旋转矩阵R，具体方法为：首先根据图像旋转角度位置关系，利用线性插值算法构造并存储降采样图像空间切面层旋转矩阵R0，随后根据降采样三维图像空间旋转矩阵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）根据步骤1）中所记录的射线源到探测板距离SDD，射线源到C臂旋转中心距离SOD，初始采样旋转角度，投影图像长度U、宽度V，投影图像像素边长h，设定的三维图像空间长度X、宽度Y、高度Z，三维图像空间体素边长l，以及步骤8）所得的到三维图像空间血管掩模，利用距离驱动算法计算各旋转角度下的完整投影矩阵； 

10）根据步骤1）所记录的投影序列图像，根据步骤9）所得到的各旋转角度下的完整投影矩阵，完成三维血管结构重建，具体方法为： 

101）设定一个重建结果向量，所述重建结果向量长度为图像空间体素总个数X×Y×Z，重建结果向量元素值全部为0；设定一个步长向量，步长向量长度为图像空间体素总个数X×Y×Z，步长向量元素值全部为0；设定一个单元向量，单元向量长度为图像空间体素总个数X×Y×Z，单元向量元素值全部为1； 

102）对于每个投影角度，进行如下操作：将单元向量与当前投影角度下的完整投影矩阵相乘，将相乘结果同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前投影角度的单方向步长向量，将该单方向步长向量累加至步长向量；对所有投影角度完成上述操作后，将累加完成的步长向量记录保存为迭代步长向量； 

103）执行重建结果向量迭代更新步骤100至300次，完成指定次数的重建结果向量迭代更新步骤后，将重建结果向量保存为重建结果三维血管体数据，所述重建结果向量迭代更新的方法如下： 

对于每个投影角度，进行如下操作：将当前的重建结果向量同当前投影角度下的完整投影矩阵相乘，得到当前方向临时投影向量，用当前方向投影图像向量减去所述当前方向临时投影向量，得到当前方向临时投影差向量，将所述当前方向临时投影差向量同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前方向差异反投影向量，用所述当前方向差异反投影向量点除以迭代步长向量，得到当前方向迭代步进向量，将当前方向迭代步进向量累加至重建结果向量；对所有角度完成上述操作后，完成一次重建结果向量更新。 

本发明中，步骤3）中利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，并将得到的初始扫描方向低分辨率体素索引投影矩阵以三元组存放方式保存于内存设备中的具体方法为： 

31）设定发射源的空间坐标为（-sSOD，0，0），设定探测板中心点坐标为（sSDD-sSOD，0，0），sSOD表示射线源到C臂旋转中心距离，sSDD表示射线源到探测板距离； 

32）对所有降采样体素作如下操作：从步骤31）设定的发射源向序号为jd的体素各角点引出射线，各射线将与降采样二维投影空间相交，其八个交点的连线所包围区域面积记sdj，所述包围区域同序号为id的像素相交面积为sd′ji，则aij=sd′ji/sdj即为低分辨率投影矩阵第id行，第jd列的元素值； 

33）将所得到的所有具有非零值的元素的行号、列号以及元素值保存，得到三元组格式存储的初始扫描方向低分辨率体素索引投影矩阵。 

本发明的步骤4）中构造并存储降采样图像空间切面层旋转矩阵R0的具体方法为：记投影采样的旋转角度为α，将三维图像空间按照z方向体素单位分组，共计Z个切面层，选取任一切面层，将其中的每个体素围绕z轴旋转-α角度，记录旋转后的位置的邻域体素编号l0，l1，l2，l3，按照线性插值方法确定矩阵第k行，第l0列，第l1列，第l2列，第l3列的四个元素的元素值，其中k为体素的序号，将所得到的所有具有非零值的元素的行号，列号以及元素值保存，得到三元组格式存储的降采样图像空间切面层旋转矩阵R0。 

本发明中，步骤9）中利用距离驱动算法计算各旋转角度下的完整投影矩阵的具体过程为： 

91）设定发射源的空间坐标为(-sSOD×cosα,-sSOD×sinα,0)，设定探测板中心点坐标为((sSDD-sSOD)×cosα,(sSDD-sSOD)×sinα,0)； 

92）对所有三维图像空间体素作如下操作：如果在步骤8所得到的三维图像空间血管掩模中，序号为j的体素没有非零值，则不做任何操作直接进入下一个体素的操作，如果具有非零值，则根据步骤91）设定的发射源向该体素各角点引出射线，各射线将与降采样二维投影空间相交，其八个交点的连线所包围区域面积记为tj，同序号为i的像素相交面积为t′ji，则aij=t′ji/tj即为投影矩阵第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字节行列号，共需要存储空间约为1.7GB；利用本发明提出的算法，在考虑图像掩模的情况下，通常血管掩模的非零元素个数在2*105左右，每个非零元素对应12个相关像素，共需要存储空间约为20Mb字节。 

同一般稀疏矩阵相比，本方法所需存储空间降至1%左右，利用获得的投影矩阵同投影向量相乘，则完成投影操作，利用获得的投影矩阵的转置，同图像向量相乘，则完成反投影操作。同一般稀疏矩阵相比本方法确定的投影矩阵构造方法大幅提高了投影、反投影计算速度，记录通过矩阵相乘运算完成投影、反投影操作所需时间，理论上，矩阵尺寸减小至1%，则速度提高应为100倍，实际应用当中，由于数据传递，数据交互等底层原因，实验证明速度提升约为20-40倍。 

获得简化后的投影矩阵后，可以利用其完成旋转X射线造影系统的三维血管重建，得到准确的三维血管结构。 

附图说明

图1是本发明方法的旋转造影系统投影生成示意图。 

图2a是本发明方法的距离驱动像素索引投影矩阵几何示意图。 

图2b是本发明方法的距离驱动像素索引投影矩阵元素权重计算示意图。 

图3是本发明方法的二值化投影分割反投影叠加示意图。 

图4是本发明方法的整体流程图。 

图5是本发明方法步骤3）的子流程图。 

图6是本发明方法步骤4）中构造并存储降采样图像空间切面层旋转矩阵R0的子 流程图。 

图7是本发明方法步骤9）中利用距离驱动算法计算各旋转角度下的完整投影矩阵的子流程图。 

具体实施方式

本发明的旋转X射线造影图像迭代重建方法，包括下列步骤： 

1）从旋转X射线造影设备读取扫描数据文件，保存投影序列图像并记录如下参数：射线源到探测板距离SDD，射线源到C臂旋转中心距离SOD，各投影采样的旋转角度，投影图像像素边长h，投影图像像长度U个像素单位和宽度V个像素单位，投影图像尺寸与二维投影空间像素尺寸一致，每个象素单位为长度为h的正方形，确定二维投影空间坐标轴u，v，坐标轴u，v分别平行于二维投影空间的长，宽方向，根据用户的精度需求，分别设定三维图像空间长度为X个体素单位、宽度为Y个体素单位、高度为Z个体素单位，每个体素单位为边长l的立方体，确定三维图像空间坐标轴x，y，z，坐标轴x，y，z分别平行于三维图像空间的长、宽、高方向，且均通过三维图像空间的中心位置，也即坐标原点为三维图像空间的中心位置。系统模拟图如图1所示，图中S为发射源位置，O为旋转中心，D为探测板中心位置。射线源到探测板距离=SD，射线源到C臂旋转中心距离=SO，记图像空间体素尺寸l×l×l，图像空间包含体素个数X×Y×Z；记投影图像像素尺寸为h×h，投影空间包含像素个数U×V；记向量 和x轴夹角为a，每次以O为中心旋转角度为Δa，总投影角度数为K。 

2）根据步骤1）中设定的三维图像空间长宽高对三维图像空间进行降采样操作，具体方法为：选定降采样倍数dif，所述降采样倍数dif为能被三维图像空间长度X、宽度Y、高度Z分别整除的整数，通常可以选择2，3，4，将三维图像空间中由dif*dif*dif个体素组成的、边长为dif*l的立方体内的体素，归并为降采样体素，对整个三维图像空间完成归并处理后，按照归并顺序排列所述降采样体素，得到降采样三维图像空间，降采样体素边长为三维图像空间体素边长的dif倍,记 Xif＝X/dif,Yif＝Y/dif,Zif＝Z/dif则降采样图像空间体素尺寸ldif×ldif×ldif，降采样图像空间包含体素个数为Xif*Yif*Zif； 

根据步骤1）中记录的投影图像长宽，对二维投影空间进行降采样操作，具体方法为：选定降采样倍数dpf，所述降采样倍数为能被二维投影空间长度U和宽度V分别整除的整数，通常可以选择2，3，4，将二维投影空间中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影空间降采样像素，对整个二维投影空间完成归并处理后，按照归并顺序排列所述投影空间降采样像素，得到降采样二维投影空间，投影空间降采样像素边长为二维投影空间像素边长的dpf倍,记Upf＝U/dpf,Vpf＝V/dpf则降采样投影空间像素尺寸hdpf×hdpf，降采样投影空间包含像素个数Upf×Vpf； 

根据步骤1）中记录的投影图像长宽，对步骤1）获得的投影序列图像进行降采样操作，得到降采样序列投影图像，具体方法为：对投影序列中每个投影图像，将其中由dpf*dpf个像素组成的、边长为dpf*h的正方形内的像素，归并为投影图像降采样像素，并且对所述正方形内的像素值进行累加求和，作为投影图像降采样像素值，对整个投影图像完成归并处理后，按照归并顺序排列所述投影图像降采样像素，得到降采样投影图像，投影图像降采样像素边长为投影图像像素边长的dpf倍,降采样投影图像像素尺寸hdpf×hdpf，降采样投影图像包含像素个数Upf×Vpf。 

3）根据步骤1）中所记录的探测板距离SDD，射线源到C臂旋转中心距离SOD，初始采样旋转角度，针对步骤2）得到的降采样三维图像空间和降采样二维投影空间，利用距离驱动算法，构造初始扫描方向低分辨率体素索引投影矩阵，具体方法为：取初始方向a=0，设定发射源的空间坐标为（-sSOD，0，0），探测板中心坐标为（sSDD-sSOD，0，0），按照x，y，z的顺序遍历降采样图像空间体素。如图2（a）所示，对于序号为(vx,vy,vz)的体素，记其编号jd＝vz×Xif×Yif+vy×Xif+vx，作点源S和编号jd体素中心点连线，其延长线与探测板相交点在探测平面的坐标为(pud，pvd)，作点源 S和编号jd体素各角点连线，其延长线与探测板相交点所包围矩形面积记为 其面积为sdj，若序号为(m，n)的投影像素点，其编号记为id＝pvd×U+pud，与 相交，且相交面积为sd′ij，，如图2（b）所示，则记低分辨率体素索引投影矩阵矩阵第id行，jd列的元素值aij=sd′ij/sdj； 

将所有具有非零值的元素的行号，列号以及元素值以三元组存放方式保存于内存设备中，得到的初始扫描方向低分辨率体素索引投影矩阵； 

4）针对步骤2）提供的降采样三维图像空间，构造各旋转角度下降采样三维图像空间旋转矩阵R具体方法为：首先对图像空间沿z轴（旋转轴）进行划分，根据图像旋转角度位置关系，利用线性插值算法构造降采样图像空间切面层旋转矩阵R0，该矩阵大小为，XifYif×XifYif为稀疏矩阵，对于所选切面层，可看作二维图像处理，以图像中心为原点，对于像素(row，col)，其2D坐标分别为(a0，b0)，其中 记其序号为j＝col×Xif+row。对应稀疏旋转矩阵R0的第j列。根据2D旋转坐标变换公式，经过旋转-a角度后该像素点的坐标位置为： 

a ′ = a 0 cos α + b 0 sin α b ′ = b 0 cos α - a 0 sin α ]]>

记 则在旋转后2D图像坐标中，像素(a'，b′)周围的4个像素点坐标分别为(a，b)，(a+1，b)，(a，b+1)，(a+1，b+1)，上述四点的序号分别为： 

i0＝b×Xif+a 

i1＝b×Xif+a+1 

i2＝(b+1)×Xif+a+1 

i3＝(b+1)×Xif+a+1 

上述序号分别对应在旋转矩阵R0的第j列，拥有非零元素的行号，矩阵元素的值分别为(a+1-a1)×(b+1-b1)，(a1-a)×(b+1-b1)，(a+1-a1)×(b1-b),(a1-a)×(b1-b)， 

将所有具有非零值的元素的行号，列号以及元素值以三元组存放方式保存于内存设备中，得到图像空间切面层旋转矩阵R0； 

根据降采样图像空间旋转矩阵R的分块对称性，利用公式 得到降采样图像空间旋转矩阵R； 

5）对步骤2）得到的降采样投影序列图像采用顶帽滤波方法进行滤波处理，对滤波处理结果二值化得到第一分割结果，具体方法为： 

对每一幅投影序列图像取反色，并进行顶帽滤波，滤波模板采用圆形掩模，根据降采样系数决定掩模直径r。在原始投影中，r一般为15-20像素，在降采样图像中，根据降采样系数，对r进行等比缩小。对顶帽变换结果进行阈值分割，分割阈值为最大灰度的1/10，得到第一分割结果； 

对步骤2）得到的降采样投影序列图像采用Frangi血管滤波方法进行滤波处理，该滤波方法为Alejandro Frangi于1998年首次在知名国际会议“Medical Image Computing and Computer-Assisted Interventation”所提出，发表于科技论文《Multiscale vessel enhancement filtering》，该方法适用于图像数据中管状结构的增强显示，被广泛应用于血管图像增强。对滤波处理结果二值化得到第二分割结果，具体方法为： 

利用Frangi滤波器对每一幅投影序列图像进行增强。增强参数为：ScaleRange=[1，r]；ScaleRatio=1.0；BetaOne=0.5，BetaTwo=1.0，其中，r等于步骤5-a顶帽变换的直径，对步骤Frangi滤波器增强结果进行阈值分割，分割阈值为最大灰度的15%，得到第二分割结果； 

对第一分割结果和第二分割结果求并集，得到低分辨率投影序列图像血管分割结果，所述步骤5）中顶帽滤波算法掩模形状与直径，分割阈值，Frangi滤波器增强参数，分割阈值参数，均为本方案的优选参数，采用其他参数仍可能得到类似结果； 

6）利用步骤3）所得到的三元组格式存储的初始扫描方向低分辨率体素索引投影矩阵和步骤4）所得到的各旋转角度降采样图像空间旋转矩阵，对步骤5）所得到的低分辨率投影序列图像血管分割结果进行反投影操作：具体方法为：对低分辨率投影图像血管分割结果，根据像素编号排列拉伸低分辨率投影图像血管分割结果，得到低分辨率投影图像血管分割结果向量，将所得向量同初始扫描方向低分辨率体素索引投影矩阵的转置矩阵相乘得到中间低分辨率反投影结果；将该中间低分辨率反投影结果同对应旋转角度降采样图像空间旋转矩阵的转置矩阵相乘，得到该旋转角度下的低分辨率反投影结果； 

7）利用步骤6）所得到的低分辨率反投影结果确定低分辨率三维图像空间血管掩模VOI，具体方法为：记旋转角度为α的低分辨率反投影结果为Uα，对于同一心动相位的少量角度投影，即投影数小于等于5的情况，掩模定义方式为： 

VOI = ∩ α U i ]]>

对于包含心脏运动的多个角度投影，即投影数大于5的情况，记总投影角度数为K，则首先定义反投影和VOI*： 

VOI * = Σ k = 1 K U k ]]>

随后用分割阈值τ对VOI*进行阈值分割，其中，分割阈值τ可以在[0，k]区间内自由确定，例如K/2，2K/3等。 

在本发明中，τ的一种优选参数为： 

τ = max ( 1 , 2 KΔα γ ) ]]>

其中，Δα为不同旋转角度之间的角度间隔，γ为完整扫描覆盖角度。 

8）对步骤7）得到的低分辨率三维图像空间血管掩模，根据步骤2）所设定的降采样倍数dif，  进行升采样操作，得到三维图像空间血管掩模，具体方法为：将低分辨率三维图像空间血管掩模的每一个低分辨率体素在x，y，z方向等分为dif份，拆分成为dif*dif*dif个高分辨率体素，将该低分辨率体素的值作为其拆分成的高分辨率体素的值，并按照拆分顺序排列高分辨率体素，当低分辨率三维图像空间血管掩模的所有低分辨率体素完成上述操作后，即得到三维图像空间血管掩模； 

9）根据步骤1）中所记录的射线源到探测板距离SDD，射线源到C臂旋转中心距离SOD，初始采样旋转角度，投影图像长度U宽度V，投影图像像素边长h，，设定的三维图像空间长度X、宽度Y、高度Z，三维图像空间体素边长l，以及步骤8）所得的到三维图像空间血管掩模，利用距离驱动算法计算各旋转角度下完整投影矩阵，计算方法同步骤3）类似，但遍历体素时仅仅考虑掩模VOI中的非零体素，此步骤不再计算旋转矩阵，通过遍历投影角度方式计算完整投影矩阵，具体方法为： 

对于每个投影角度α，点源的中心坐标为：(-sSOD cosα,-sSOD sinα,0)，探测板中 心坐标为((sSDD-sSOD)cosα,(sSDD-sSOD)sinα,0)，按照x，y，z的顺序遍历图像空间体素，对于序号为(vx,vy,vz)的体素，记其编号j＝vz×X×Y+vy×X+vx，如果步骤8）所得血管掩模在编号j的体素具有非零值，则作点源S和编号j图像空间体素中心点连线，其延长线与探测板相交点在探测平面的坐标为(pu，pv)，作点源S和编号j图像空间体素各角点连线，其延长线与探测板相交点所包围矩形面积记为 其面积为sj，若序号为(m，n)的投影像素点，其编号记为i＝pv×U+pu，与 相交，且相交面积为s′ji，，如图2（b）所示，则记矩阵第i行，j列的元素值aij=s′ji/sj； 

将所有具有非零值的元素的行号，列号以及元素值以三元组存放方式保存于内存设备中，得到投影角度α下的体素索引投影矩阵； 

10）根据步骤1）所记录的投影序列图像，根据步骤9）所得到的各旋转角度下的完整投影矩阵，完成三维血管结构重建，具体方法为： 

101）设定一个重建结果向量，所述重建结果向量长度为图像空间体素总个数X×Y×Z，重建结果向量元素值全部为0；设定一个步长向量，所述步长向量长度为图像空间体素总个数X×Y×Z，步长向量元素值全部为0；设定一个单元向量，所述单元向量长度为图像空间体素总个数X×Y×Z，单元向量元素值全部为1； 

102）对于每个投影角度，进行如下操作：将单元向量与当前投影角度下的完整投影矩阵相乘，将相乘结果同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前投影角度的单方向步长向量，将该单方向步长向量累加至步长向量；对所有投影角度完成上述操作后，将累加完成的步长向量记录保存为迭代步长向量； 

103）执行重建结果向量迭代更新步骤100至300次，完成指定次数的重建结果向量迭代更新步骤后，将重建结果向量保存为重建结果三维血管体数据，所述重建结果向量迭代更新的方法如下： 

对于每个投影角度，进行如下操作：将当前的重建结果向量同当前投影角度下的完整投影矩阵相乘，得到当前方向临时投影向量，用当前方向投影图像向量减去所述当前方向临时投影向量，得到当前方向临时投影差向量，将所述当前方向临时投影差向量同当前投影角度下的完整投影矩阵的转置矩阵相乘，得到当前方向差异反投影向量，用所述当前方向差异反投影向量点除以迭代步长向量，得到当前方向迭代步进向 量，将当前方向迭代步进向量累加至重建结果向量；对所有角度完成上述操作后，完成一次重建结果向量更新。 

步骤10）中，具体迭代次数需要根据数据条件和重建效果要求设定，投影数据对比度增加，则所需迭代次数减小，迭代次数增加，相应重建效果提升，本方案中的优选参数为200次。