CN114119801A

CN114119801A - 三维数字减影血管造影方法、装置、电子设备及存储介质

Info

Publication number: CN114119801A
Application number: CN202111448108.3A
Authority: CN
Inventors: 于留青; 刘肖琳
Original assignee: Individual
Current assignee: Individual
Priority date: 2021-11-30
Filing date: 2021-11-30
Publication date: 2022-03-01

Abstract

本申请实施例提供了一种三维数字减影血管造影方法、装置、电子设备及存储介质。首先提出一种用X射线成像***拍摄多帧带有标志物的被测物体图像，通过对图像中标志物的计算，得到成像***与被测物体之间高精度成像几何关系，得到被测物体三维图像的CT图像重建方法。用CT图像重建方法分别对掩模运行和充盈运行下拍摄的多帧图像进行CT图像重建，将充盈运行与掩模运行两种三维图像进行对应坐标像素灰度相减，得到减影血管造影图。由于坐标系建立在被测物体上，被测物体坐标系包括了被测物体的运动或变动，因此，被测物体三维图像结果是高精度匹配对齐的，自动有效消除CT运行状态的各种运动和误差，高精度得到三维数字减影血管造影图像。

Description

三维数字减影血管造影方法、装置、电子设备及存储介质

技术领域

本申请涉及图像处理技术领域，具体涉及一种三维数字减影血管造影方法、装置、电子设备及存储介质。

背景技术

三维数字减影血管造影技术(3D Digital Subtraction Angiography，3D-DSA)是近年来快速发展的一种新型血管造影技术，用于检测、监测人体心脏、大脑等部位血管的三维结构、尺寸、工况等。

它的主要原理是在掩模状态下，通过C形臂CT或螺旋CT等CT技术对人体的感兴趣部位进行CT三维图像重建，得到掩模状态下的三维CT图像；然后对人体注入造影剂，并进行三维图像重建，得到充盈状态下的三维CT图像。将两个三维CT图像进行配准、相减，得到三维数字减影血管造影图像。

目前，成像***与被测物体之间的成像几何关系的精度，以及图像减影的配准精度主要靠CT成像***的结构刚度、加工精度和运动同步控制精度等硬件的性能等来保证的，并通常要求被测物体保持稳定，这种方式成本高，且精度低。因此，如何提高三维DSA中的成像几何关系标定精度和后续图像减影的配准精度，是当前亟待解决的问题。

发明内容

本申请实施例提供了一种三维数字减影血管造影方法、装置、电子设备及存储介质，以放宽各类CT机加工与运动精度要求，提高成像***标定精度，并且方便可靠、高精度的得到三维数字减影血管造影图像。

第一方面，本申请实施例提供一种三维数字减影血管造影方法，包括：

获取多帧第一投影图像，所述多帧第一投影图像是血管处于掩膜状态下，CT成像***在多个第一成像角度下向被测物体发射X射线进行成像得到，所述被测物体上设置有多个标志物，所述CT成像***在预定轨道上运动；

获取至少一组第二投影图像，其中，每组第二投影图像中包括多帧第二投影图像，所述每组第二投影图像中的多帧第二投影图像是血管处于充盈状态下，所述CT成像***在多个第二成像角度下向被测物体发射X射线进行成像得到；

根据所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像对应的成像几何关系；

根据所述每组第二投影图像中的多帧第二投影图像中的标志物的像素坐标，确定所述每组第二投影图像中的多帧第二投影图像对应的成像几何关系；

根据所述多帧第一投影图像分别对应的成像几何关系，对所述多帧第一投影图像进行CT三维图像重建，得到第一CT三维图像；

根据所述每组第二投影图像中的多帧第二投影图像对应的成像几何关系，对所述每组第二投影图像中的多帧第二投影图像进行CT三维图像重建，得到至少一幅第二CT三维图像；

根据所述第一CT三维图像和所述至少一幅第二CT三维图像，得到至少一幅三维数字减影血管造影图像。

第二方面，本申请实施例提供一种三维数字减影血管造影装置，包括：

获取单元，用于获取多帧第一投影图像，所述多帧第一投影图像是血管处于掩膜状态下，CT成像***在多个第一成像角度下向被测物体发射X射线进行成像得到，所述被测物体上设置有多个标志物，所述CT成像***在预定轨道上运动；

处理单元，用于根据所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像对应的成像几何关系；

第三方面，本申请实施例提供一种电子设备，包括：处理器，所述处理器与存储器相连，所述存储器用于存储计算机程序和所述投影射线图像，所述处理器用于执行所述存储器中存储的计算机程序，以使得所述电子设备执行如第一方面所述的方法。

第四方面，本申请实施例提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序使得计算机执行如第一方面所述的方法。

第五方面，本申请实施例提供一种计算机程序产品，所述计算机程序产品包括存储了计算机程序的非瞬时性计算机可读存储介质，所述计算机可操作来使计算机执行如第一方面所述的方法。

实施本申请实施例，具有如下有益效果：

传统的各类CT设备中，被测物体与成像***之间的成像几何关系主要是靠机械装置或设备硬件来完成和保证精度的。本申请是通过摄影测量等算法获取成像***与被测物体之间的成像几何关系，并且高精度的优化获取到的成像几何关系。由于本申请的全局坐标系与被测物体固连，成像几何关系自动地包含了成像***的位置姿态，包括了实际工况下振动、晃动等引起的预定轨道的运动误差。特别是在数字减影血管造影运行中的掩模状态和充盈状态两种三维CT图像重建过程中，保持被测物体上的标志物与被测物体的相对位置不变，使多次投影成像环节具有相同的坐标系，自动消除由于被测物体运动变化等引起的成像几何关系的误差，从而提高了成像几何关系的标定精度，进而提高了CT图像重建精度。由于提高了CT成像重建精度，特别是提高了不同状态的体积图像匹配的对应性精度，从而有效提高数字减影血管造影图像的质量和精度。

附图说明

图1为本申请实施例提供的一种三维数字减影血管造影方法的流程示意图；

图2为本申请实施例提供的一种获取多帧第一投影图像的示意图；

图3为本申请实施例提供的一种成像***的内参数和外参数的示意图；

图4为本申请实施例提供的一种三维数字减影血管造影装置的功能单元组成框图；

图5为本申请实施例提供的一种电子设备的结构示意图。

具体实施方式

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

本申请的说明书和权利要求书及所述附图中的术语“第一”、“第二”、“第三”和“第四”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、***、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤或单元。

为了便于理解本申请的实施例，先对本申请涉及到的相关知识进行介绍。

在锥束CT(CBCT)成像中，其投影成像过程可以通过公式(1)表示：

p＝Sf 公式(1)

其中，S是反映成像***与被检体之间的成像几何关系的***矩阵，f为待重建的三维被检体，p为投影图像。

其中，CT图像重建过程可以通过公式(2)表示：

f＝S^-1p 公式(2)；

然而，由于有限角度重建具有病态性，所以用广义逆处理是病态性的，则待重建图像的估计

可以表示为：

从公式(3)式可以看到，成像***与被检物体之间的成像几何关系的获得是CT图像重建的前提，成像几何关系的精度是影响CT重建质量的关键因素。

目前，三维DSA的常规主要设备之一是C型臂CT，然而该CT需要围绕被测物体在预定轨道上高速运动，且CT设备重量重、运动快、轨道曲率半径较小，很容易引起振动、晃动等误差因素，这些因素降低了该CT设备的三维CT图像重建的质量，以及掩模运行与充盈运行两个状态下的三维CT图像的配准精度。

如何提高预定轨道CT的DSA过程中的成像几何关系的标定精度，以提高图像重建质量和后续图像减影配准对齐的精度，特别是如何提高由于被测物体运动引起的配准对齐的精度，是当前亟待解决的问题。

为了便于理解本申请的技术方案，先对本申请相关概念进行解释和说明。

本申请的CT成像***是一种满足中心透视投影的成像模型的成像***。其中，该成像***包括射线源和射线探测器，射线源是中心透视投影的焦点，并且在本申请中射线源和射线探测器之间的相对距离保持不变。

本申请的成像***可以为射线源和射线探测器之间相对距离保持不变的任意类型的CT成像***。包括但不限定：各类预定轨道的C型臂CT、螺旋CT、合成X线体层成像(Tomosynthesis)CT以及基于DR的锥形束计算机断层扫描成像(Cone beam ComputedTomography，CBCT)等。

本申请的坐标系包括全局坐标系、被测物体坐标系、成像***坐标系以及图像坐标系。其中，被测物体坐标系的坐标原点建立在被测物体上或与被测物体固连的承载装置上。为了便于计算，本申请中将被测物体坐标系作为全局坐标系。成像***坐标系的坐标原点建立在射线源中心，且z轴与成像***光轴重合。标志物的空间坐标为标志物在被测物体坐标系下的三维空间坐标，标志物的像素坐标为标志物在图像坐标系下的二维像素坐标。

本申请中被测物体在多个成像角度下的多帧投影图像，由于被测物体上设置有多个标志物，则每帧投影图像含有若干个标志物。用图像处理方法识别提取各投影图像上的标志物，得到标志物的像素坐标，并用图像处理数学方法将各不同投影图像间的同名标志物(同名点)进行匹配编号，例如RANSAC等算法。为了方便这些投影图像间同名标志物的自动匹配，可以将其中多个标志物做成不同的特定形状，例如，三角形，十字形等。对于各帧投影图像，采用图像处理和模式识别算法检测出这种特定形状的标志物时，可以将它们自动匹配配对编号，其余的标志物可以通过与该特定标志物的相对位置关系来匹配编号，这样就快速准确的识别出了每帧投影图像中所包含的标志物以及标志物的编号，从而快速得到不同的投影图像中的相同的标志物。

可选的，多个标志物中可以有两个以上标志物的空间位置关系或相对距离是已知的。

可选的，标志物的材质可以为铅、钢或骨质材料等对x射线敏感的任何材料，形状可以为点珠状或线状等，本申请不对标志物的材质和形状进行限定。

应理解，将标志物设置于被测物体上，可以是将标志物直接设置于被测物体上，也可以是将标志物间接设置于被测物体上，比如，可以将标志物设置于附件上，将附件设置于被测物体上，且该附件是可以自由拆卸的附件。因此，本申请不对标志物设置于被测物体上的方式进行限定。此外，该多个标志物之间的位置关系可以是任意的。比如，该多个标志物可以随机地固连安置在被测物体上，跟随被测物体运动；多个标志物的相互关系也可以是固定的，例如，将多个标志物安置在某个架子上，然后将架子固连在被测物体上。

为了精确地获取投影图像中的标志物的像素坐标，通常采用亚像素定位算法，例如，首先对投影图像中的标志物的灰度曲面或相关系数矩阵进行二次或三次曲面拟合，然后求极值，将极值点作为标志物的像素坐标。

为了减少标志物的射线阻挡效应影响，当使用图像处理算法确定出投影图像中的各标志物并且匹配编号对齐后，可以将该标志物的像素值用该标志物位置周边的像素值进行插值，取代该标志物的原有像素值，从而抑制掉该射线敏感标志物在射线图像上的影响。也可以将标志物的像素值置为空，即不参加三维图像重建的计算，减少标志物的射线敏感材料的扰动影响，从而不影响三维图像重建的质量。

参阅图1，图1为本申请实施例提供的一种三维数字减影血管造影方法的流程示意图。该方法包括以下步骤：

101、获取多帧第一投影图像。

其中，多帧第一投影图像是血管处于掩膜状态下(即人体中未注入造影剂的状态下)，CT成像***在多个第一成像角度下向被测物体发射X射线进行成像得到。CT成像***在预定轨道上运动，即CT成像***中的射线源和射线探测器之间相对距离保持不变。

其中，多帧第一投影图像与多个第一成像角度一一对应，多个第一成像角度可以是任意的成像角度。示例性的，在任意一个第一成像角度下CT成像***向被测物体发射X射线，X射线照射到被测物体，射线探测器在射线源的对应方位上记录一系列投影数据，基于该一系列投影数据生成与该成像角度对应的一帧第一投影图像，进而得到多帧第一投影图像。

进一步的，被测物体上设置有多个标志物，且每帧第一投影图像中包含有该多个标志物中的部分或全部。

具体的，控制CT成像***进行多次转动和/或平移，每次转动和/或平移得到一个成像角度，则得到一个第一成像角度；并在多次转动和/或平移被测物体后的多个第一成像角度下分别控制CT成像***向被测物体发射X射线，得到多帧第一投影图像。

如图2所示，将被测物体设置于承载装置上，则可控制CT成像***在预定轨道上绕被测物体旋转和/或平移，在每次转动和/或平移之后，控制射线源向被测物体发射X射线，射线探测器在射线源的对应方位上记录一系列投影数据，生成与每次转动和/或平移之后的投影图像，得到多帧第一投影图像。

需说明，本申请的成像***是在预定轨道上运动，当成像***的运动轨迹已知时，获取多帧第一投影图像的方式就是控制CT成像***在预先设定的轨道上运动，在到达多个预先设定的位置或第一成像角度时，向被测物体发射X射线，得到多帧第一投影图像。

在本申请的一个实施方式中，承载装置是可转动或者平移的装置，则可控制成像***和被检测物体同时进行多次转动和/或平移，并且在每次转动和/或平移后控制射线源向被检测物体发射X射线，得到多帧第一投影图像。

102：获取至少一组第二投影图像。

示例性的，与获取第一投影图像的方式类似，针对每组第二投影图像来说，在被测物体血管处于充盈状态(即注入造影剂)的情况下，控制CT成像***在多个第二成像角度下向被测物体发射X射线，得到该组第二投影图像中的多帧第二投影图像。应说明，获取任意一组第二投影图像中的多帧第二投影图像的方式，与上述获取多帧第一投影图像的方式类似，不再叙述。

应说明，任意两组第二投影图像所对应的多个第二成像角度之间可以相同，也可以不同；并且，任意一组投影图像对应的多个第二成像角度和多个第一成像角度可以相同，也可以不同。也就是说，每次进行三维CT图像重建过程，所使用的成像角度可以是相同的，也可以是不同的，本申请对此不做限定。

举例来说，在A、B、C、D四个位置分别获取了四帧第一投影图像，则在获取任意一组第二投影图像时，可以在该A、B、C、D四个位置进行成像，也可以不在该A、B、C、D四个位置中的任意一个位置处进行成像。

应说明，在不同次的投影成像环节上，需要保持标志物与被测物的相对几何关系不变，直到若干次投影成像环节结束，以便使多次投影成像环节具有相同的坐标系。

应说明，CT成像***在掩膜状态和充盈状态两种模式运行时，CT成像***围绕被测物体获取投影图像，借助相同的C形臂X射线设备成像***来采集记录数据组。可选的，CT成像***在轨道上还可以在掩模运行和充盈运行之间返回运行，以便可以在充盈运行中再次从与掩模运行时相同的起始位置启动轨道上运行。例如，通常利用C形臂X射线设备以锥形束几何形状的形式执行运转(锥形束CT运行)。在C形臂围绕被测物体旋转期间，在不同的旋转角度记录多个投影图像。可选的，至少旋转200°。

103：根据每帧第一投影图像中的标志物的像素坐标，确定每帧第一投影图像对应的成像几何关系。

其中，每帧第一投影图像对应的成像几何关系为获取每帧第一投影图像时被测物体与成像***之间的成像几何关系，即成像***坐标系相对于被测物体坐标系的姿态和位置，且后面所提到的第一投影图像的成像几何关系均与此处的解释类似，不再叙述。

示例性的，每帧第一投影图像对应的成像几何关系包括每帧第一投影图像对应的内参数和外参数。如图3所示，每帧投影的内参数包括获取每帧第一投影图像对应的焦距，即射线源C到射线探测器的距离D_SD和主点坐标(C_x，C_y)；其中，主点为射线源C到射线探测器成像面的垂直交点，射线源与主点的连线定义为光轴，主点坐标为每帧第一投影图像的中心像素点的像素坐标；外参数是获取每帧第一投影图像时成像***坐标系相对于被测物体坐标系的姿态和位置关系，且姿态通过旋转矩阵表示，位置关系通过平移矩阵表示。

进一步的，被测物体坐标系、成像***坐标系和图像坐标系之间的成像几何关系可以通过公式(4)表示：

其中，Z_C为一个常数因子，(X，Y，Z)指标志物P在被测物体坐标系中的空间坐标，

指标志物P经中心透视投影得到的像点p的像素坐标，即标志物P在投影图像中的像素坐标；du和dv分别为射线探测器像元尺寸；R₀为旋转矩阵，该旋转矩阵为一个3×3单位正交阵，它的元素是旋转角的三角函数组合，旋转角用于表示将被测物体坐标系的姿态变换到与成像***坐标系的姿态一致时分别绕三个坐标轴转过的欧拉角。T₀为平移矩阵，T₀＝(T_x,T_y,T_z)称为平移向量，T_x，T_y，T_z分别表示将被测物体坐标系的原点移至成像***坐标系的原点在三个坐标轴上的平移量，T₀也可以称为被测物体坐标系原点在成像***坐标系中的坐标，等价为成像***坐标系的原点在被测物体坐标系中的坐标。

上述可以看出本申请所要获取的每帧第一投影图像对应的成像几何关系就是获取每帧第一投影图像对应的主点坐标、射线源到射线探测器的距离、旋转矩阵以及平移矩阵。此外，本申请中将每帧第一投影图像对应的主点坐标(Cx，Cy)的初值设置为该帧第一投影图像的中心像素点的像素坐标，所以，每帧第一投影图像对应的主点坐标可以基于每帧第一投影图像直接确定出来。而且，不失一般性，在本申请实施例中，假设整个成像过程中射线源C到射线探测器的距离D_SD保持不变，因此任意两帧第一投影图像所对应的焦距相同，因此第一次确定出了射线源C到射线探测器的距离D_SD后，就无需再次确定D_SD；则该多帧第一投影图像均使用第一次确定出的D_SD。为了简化计算，也可以将CT***的标称D_SD作为计算模型中的D_SD的初值。

示例性的，根据CT成像***在预定轨道上的运动位置姿态(即CT成像***的标称参数)，确定每帧第一投影图像对应的初始成像几何关系；根据每帧第一投影图像的初始成像几何关系以及每帧第一投影图像中的标志物的像素坐标，确定每帧第一投影图像的成像几何关系。

下面介绍四种获取每帧第一投影图像的初始成像几何关系以及每帧第一投影图像中的标志物的空间坐标的方法。

方法1：成像***的标称参数是已知的，且每帧第一投影图像中的标志物的空间坐标全部已知时。

由于成像***的标称参数是已知的，也就是说成像***在每个第一成像角度下的初始成像几何关系是已知的，因此可根据成像***的标称参数直接确定出每帧第一投影图像对应的初始成像几何关系；

其中，所谓标志物的空间坐标是已知的，即标志物相对于被测物体坐标系的坐标原点之间的相对位置姿态关系是设定好的，相对位置姿态关系不会发生变化。因此，当建立好被测物体坐标系之后，可直接得到各标志物的空间坐标。示例性的，可以预先设置好多个标志物相对于被测物体坐标系的坐标之间的相对位置关系，就可使每帧第一投影图像中的标志物的空间坐标全部已知。

方法2：成像***的标称参数已知，但多帧第一投影图像中有些投影图像中存在未知空间坐标的标志物。

示例性的，多帧第一投影图像中的k帧第一投影图像中存在未知空间坐标的标志物，0＜k≤n(n为多帧第一投影图像的数量)，即部分或全部第一投影图像中存在空间坐标未知的标志物，对于存在空间坐标未知的标志物的第一投影图像来说，其上的所有标志物的空间坐标可以均未知，也可能是部分标志物的空间坐标未知，另外一部分标志物的空间坐标已知。

下面以第一投影图像A为例说明获取第一投影图像A中空间坐标未知的标志物的空间坐标的过程，第一投影图像A为k帧第一投影图像中的任意一帧。应理解，k帧第一投影图像中其他帧第一投影图像中的空间坐标未知的标志物的空间坐标的获取过程与此类似，不再叙述。

示例性的，从多帧第一投影图像中选取第一投影图像B，第一投影图像A与第一投影图像B之间的成像角度大于角度阈值，其中，第一投影图像B与第一投影图像A具有p个相同标志物，且该p个相同标志物中包含了第一投影图像A中空间坐标未知的标志物。

在本申请的一个实施方式中，由于成像***的标称参数是已知的，则可根据成像***的标称参数确定出第一投影图像A和第一投影图像B对应的初始成像几何关系；最后，根据p个相同标志物分别在第一投影图像A与第一投影图像B中的像素坐标，以及第一投影图像A和第一投影图像B对应的初始成像几何关系，确定出p个相同标志物的空间坐标；若投影图像中仍然还存在未知空间坐标的标志物，则可以按照上述选取第一投影图像A和第一投影图像B的方式，得到第一投影图像A中未知空间坐标的标志物的空间坐标。并且使已知空间坐标的标志物个数大于或等于第一阈值。

在本申请的另一个实施方式中，根据CT成像***在预定轨道上的运动轨迹位姿，确定第一投影图像A和第一投影图像B在成像时成像***之间的相对旋转平移矩阵；根据相对旋转平移矩阵、p个相同的标志物分别在第一投影图像A以及第一投影图像B中的像素坐标，得到p个标志物中的未知空间坐标的标志物的空间坐标，即进行三角交会，得到p个相同标志物的空间坐标。同样，若投影图像中仍然还存在未知空间坐标的标志物，则可以按照上述选取第一投影图像A和第一投影图像B的方式，获取未知空间坐标的标志物的空间坐标，可得到第一投影图像A中未知空间坐标的标志物的空间坐标。

对于其他帧第一投影图像来说，也可以按照上述的方法，求取其他帧第一投影图像中各个标志物的空间坐标。应说明，在求取任意一帧第一投影图像中的标志物的空间坐标时，若其中包含的某个标志物的空间坐标已知了(比如，在求取其他帧第一投影图像的标志物的三维坐标时，已经求取出了该标志物的空间坐标)，则不需要再次求取该标志物的空间坐标，只需要求取其他帧第一投影图像中未知空间坐标的标志物的空间坐标即可，这样每帧第一投影图像中的标志物的空间坐标均是已知的。

应说明，基于成像***之间的相对旋转平移矩阵得到的标志物的空间坐标是指定成像***坐标系下的坐标，根据应用需求，指定建立与被测物体固连的全局坐标系，得到各标志物的全局空间坐标。又由于标称参数已知，则每帧第一投影图像对应的初始成像几何关系也是知道的了。

方法3：成像***的标称参数未知，但每帧第一投影图像中的标志物的空间坐标是已知的。

由于成像***的标称参数未知，无法直接得到每帧第一投影图像的初始成像几何关系；但是，所有标志物的空间坐标是已知的，因此每帧第一投影图像对应的内外参数(即初始成像几何关系)可以根据公式(4)直接使用每帧第一投影图像中包含的标志物的空间坐标以及像素坐标构成该帧投影图像的投影矩阵元素的线性方程组，求解该方程组，再根据投影矩阵与成像***内外参数之间的关系，从投影矩阵分解得到该帧投影图像对应的初始成像几何关系。

求解该成像几何内外参数的方法，根据已知标志物的个数不同可有多种数学方法，例如PnP算法等。可选的，当特别设计满足一定约束条件的已知空间坐标的标志物的分布时，最少只需要3个已知空间坐标的标志物即可实现该帧投影图像的内外参数的求解。当直接采用***设备的射线源到探测器距离的标称参数作为内参数时，可直接用3个已知空间坐标的标志物即可实现该帧投影图像的外参数的求解。

方法4：成像***的标称参数未知，且每帧投影图像中已知空间坐标的标志物的数量均少于第一阈值时。

示例性的，从多帧第一投影图像中选取第一投影图像A和第一投影图像B组成第一投影图像对，第一投影图像A与第一投影图像B之间的成像角度差值大于角度阈值，其中，第一投影图像A和第一投影图像B具有w个相同的标志物，并将这w个相同的标志物作为w个第一标志物，w大于或等于第二阈值。

虽然成像***的标称参数未知，但是由于成像***在预定轨道上运动，因此，可以基于每帧第一投影图像的成像位置，以及成像***在预定轨道上的运动轨迹，确定第一投影图像A和第一投影图像B成像时所述成像***之间的相对旋转平移矩阵；然后，根据w个第一标志物分别在第一投影图像A和第一投影图像B中的像素坐标以及相对旋转平移矩阵，确定w个第一标志物的空间坐标；根据w个第一标志物分别在第一投影图像A中的像素坐标以及w个第一标志物的空间坐标，确定第一投影图像A的初始成像几何关系；根据w个第一标志物分别在第一投影图像B中的像素坐标以及w第一标志物的空间坐标，确定投影图像B的初始成像几何关系。

这样就可以基于三角交会的方法，求取出两帧第一投影图像的初始成像几何关系。有了初始成像几何关系的第一投影图像之后，可以从剩余未知成像几何关系的第一投影图像中选取一帧，以及从已知成像几何关系的第一投影图像中选取一帧，组成一个新的投影图像对，这两帧第一投影图像中具有相同的标志物，且相同标志物中已知空间坐标的标志物的数量大于或等于第一阈值。这样基于pnp方法，就可以直接求取出未知成像几何关系的第一投影图像的初始成像几何关系。若这两帧第一投影图像中相同的标志物中存在未知空间坐标的标志物，则可以基于这两帧第一投影图像的初始成像几何关系，以及未知空间坐标的标志物分别在两帧第一投影图像中的像素坐标，确定出该存在未知空间坐标的标志物的空间坐标；依次循环，可以得到每帧第一投影图像对应的初始成像几何关系，以及每帧第一投影图像中的标志物的空间坐标。

方法5：成像***的标称参数未知，且部分投影图像中已知空间坐标的标志物的数量大于或等于第一阈值。

示例性的，部分标志物的空间坐标已知，也就造成部分投影图像中包含有已知空间坐标的标志物的数量大于或等于第一阈值，另外一部分投影图像中包含有已知空间坐标的标志物的数量小于第一阈值；因此，针对部分已知的情况，获取多帧投影图像中的目标投影图像，其中，目标投影图像中包含已知空间坐标的标志物数量大于或者等于第一阈值；基于已知空间坐标的标志物的空间坐标和像素坐标以及上述的公式(4)，可得到目标投影图像的初始成像几何关系；

在确定出目标投影图像的初始成像几何关系之后，可以将目标投影图像移入到已知成像几何关系的投影图像集合中，将除目标投影图像之外的剩余投影图像移入到未知成像几何关系的投影图像集合，执行多次位姿估计操作，即，使用已知成像几何关系的投影图像集合中的投影图像与未知成像几何关系的投影图像重新组成一个投影图像对；基于2D-3D位姿估计方法，去估计未知成像几何关系的投影图像集合中的投影图像的成像几何关系，得到每帧第一投影图像的初始成像几何关系。

其中，上述多次位姿估计操作中的第i次位姿估计操作包括以下步骤内容：

从投影图像集合h_i-1中选取第一投影图像A，以及从投影图像集合w_i-1中选取第一投影图像B，其中，第一投影图像A和第一投影图像B具有p个相同的标志物，并将该p个相同的标志物作为p个第二标志物，其中，投影图像集合h_i-1和投影图像集合w_i-1为执行第i-1次位姿估计操作得到，其中，投影图像集合h_i是由多帧第一投影图像中已知成像几何关系的投影图像组成；投影图像集合w_i包括多帧投影图像中除投影图像集合h_i之外的所有投影图像，即投影图像集合w_i为多帧投影图像中成像几何关系未知(主要是外参数未知)的投影图像组成。当i＝1时，则投影图像集合h₀中包含的投影图像为目标投影图像，投影图像集合w₀为多帧第一投影图像中除目标投影图像之外的剩余投影图像；

当p个第二标志物中包含有已知空间坐标的第二标志物的数量大于或等于第三阈值时，比如，可以设置第三阈值为3，基于2D-3D位姿估计方法，根据成像***的内参数，p个第二标志物中已知空间坐标的第二标志物的空间坐标，以及已知空间坐标的第二标志物在第一投影图像B中的像素坐标，确定第一投影图像B对应的外参数。

至此获取到了第一投影图像B的初始成像几何关系，因此将第一投影图像B从投影图像集合w_i-1中移出，放入到投影图像集合h_i-1，分别得到与第i次位姿估计操作对应的投影图像集合h_i和投影图像集合w_i，即将已知成像几何关系的投影图像移到投影图像集合h_i中。应理解，当执行完N次位姿估计操作后，投影图像集合w_N为空集，即求取出每帧第一投影图像的初始成像几何关系。

进一步的，针对p个第二标志物中不知空间坐标的第二标志物，则可以基于第一投影图像A和第一投影图像B的初始成像几何关系，确定该第二标志物的空间坐标，从而确定出p个第二标志物的空间坐标。因此，执行每次位姿估计操作之后，如果参与位姿估计操作的两帧投影图像上相同标志物中存在未知空间坐标的标志物，则可基于位姿估计操作得到的成像几何关系，计算出未知空间坐标的标志物的空间坐标，这样随着位姿估计操作执行的次数增多，则每帧第一投影图像上会有更多的标志物的空间坐标被计算出来。因此，若任意一个标志物至少在两帧投影图像中出现，则执行完多次位姿估计操作之后，每帧第一投影图像上的标志物的空间坐标都可以计算出来。

示例性的，经过上述五种处理方法，可获取第一投影图像的初始成像几何关系以及多个标志物中各个标志物的空间坐标。进一步的，根据每帧第一投影图像对应的初始成像几何关系，以及每帧第一投影图像中的标志物的像素坐标，和各标志物的空间坐标进行优化处理，得到每帧第一投影图像的成像几何关系。

示例性，采用光束法平差优化方法，高精度地优化得到每帧第一投影图像的成像几何关系和被测物体上的标志物在被测物体坐标系下的空间坐标。

具体的，光束法平差优化方法中，根据每帧投影图的初始成像几何关系中的内参数、外参数以及各个标志物的空间坐标，按成像模型重新计算每帧第一投影图像中的各个标志物在每帧第一投影图像中的投影(标志物在投影图像中的像素坐标)，即重投影像点，并以各个标志物的重投影结果与实际投影结果之间的距离之和最小作为优化的目标函数。对于目标函数的最优求解，可通过泰勒展开对观测方程和约束方程进行线性化，再逐步迭代计算各平差参数的修正值。如果算法收敛，则迭代计算的修正值会逐步趋于零。当最后的修正值小于阈值时，确定成像关系已经被很好的满足，得到了优化后的成像几何关系，以及优化后的各个标志物的空间坐标。

具体的，x_ij为多个标志物中的第i个标志物在第j帧投影图像中的实际投影结果。为了标记上的方便，假定每一个标志物在所有的投影图像中可视。因此，可确定目标函数为最小化所有标志物的重投影像点和实际投影像点之间的距离平方和。应理解，先计算每帧第一投影图像中各个标志物的重投影像点和实际投影像点之间的距离平方和；然后，将多帧投影图像的距离平方和作为目标函数。当然，也可以计算每帧第一投影图像中各个标志物的重投影像点和实际投影像点之间的差值的绝对值；然后，将多帧投影图像的绝对值之和作为目标函数。本申请中主要以距离平方和作为目标函数为例进行说明。

其中，目标函数可以通过公式(5)表示：

其中，m为多个标志物的数量，n为多帧投影图像的数量，x_ij和

分别表示第i个标志物在第j帧投影图像中的实际投影像点和重投影像点，d(x，y)为非齐次坐标表示的两个像点x和y之间的欧氏距离，k表示投影图像对应的内参数组成的列向量。a_j为第j帧投影图像对应的外参数组成列向量，b_i为第i个标志物的空间坐标组成的列向量。

进一步的，由于多帧第一投影图像均为同一成像***所生成，则多帧第一投影图像上的内参数都相同，上述目标函数有特殊的参数划分形式，相应的正规方程有特殊的稀疏分块形式。可采用共内参多视图物点形式的稀疏LM光束法间接平差，优化每帧投影的内参数、外参数和各个标志物的空间坐标。将上述目标函数取值最小时每帧第一投影图像的内参数、外参数和各个标志物的空间坐标，作为相应每帧第一投影图像的成像几何关系和各个标志物的空间坐标。上述确定最小化所有标志物的重投影像点和实际投影像点之间的距离平方和为目标函数的过程既可以在像方进行也可以在物方进行。

实际优化计算过程中，参加优化迭代的参数可以是不同的，例如当内参数和/或标志物的空间坐标是准确的时，内参数和/或标志物的空间坐标可以不参加优化计算。这里的优化目标函数可以有不同的具体形式，例如，距离平方和可以用距离绝对值和来代替等。

应理解，上述给出了使用所有的第一投影图像进行优化的方法，在实际应用中，也可以使用多帧第一投影图像中的部分第一投影图像进行优化，其优化方式与使用所有的第一投影图像的优化方式类似，不再叙述。

应说明，当多个标志物组中有两个以上标志物的空间位置关系或相对距离是已知时，可以将这些已知的条件关系，作为优化目标函数的一项或几项联立约束条件，进行联合求解，可以使优化求解过程的最优解更加稳定，精度更高。可选的，CT成像***的对成像几何关系有比较精确的特殊状态，例如，CT成像***的初始或/和终止位置的状态值，将该特殊状态值作为优化目标函数的联立约束条件，也可以使优化求解过程更稳定，精度更高。

104：根据每组第二投影图像中的多帧第二投影图像中的标志物的像素坐标，确定每组第二投影图像中的多帧第二投影图像对应的成像几何关系。

其中，每组第二投影图像中任意一帧第二投影图像对应的成像几何关系为获取该帧第二投影图像时，被测物体与CT成像***之间的相对成像几何关系。

针对每组第二投影图像来说，获取该组第二投影图像中任意一帧第二投影图像的成像几何关系与上述获取每帧第一投影图像的成像几何关系的方式类似，不再详细叙述。实际操作时，由于在充盈状态下的标志物的空间坐标相对于掩膜状态下的标志物的空间坐标是没有发生变化的，因此，可以认为多个标志物的空间坐标都是已知的，也就是说，任意一组第二投影图像中的任意一帧第二投影图像中的标志物的空间坐标都是已知的。

因此，当成像***的标称参数是已知时，则可以通过上述的方法1获取任意一帧第二投影图像的初始成像几何关系，然后，再通过上述的优化方法，确定每帧第二投影图像的成像几何关系；当成像***的标称参数未知时，则可以通过上述的方法3获取任意一帧第二投影图像的初始成像几何关系，然后再通过上述的优化方法，确定每帧第二投影图像的成像几何关系。

应说明，即使任意一组第二投影图像中的任意一帧第二投影图像中的标志物的空间坐标全部都是已知的，也可以将它们当做未知的，按照上述方法2和方法4，确定任意一帧第二投影图像的初始成像几何关系，然后再通过上述的优化方法，确定每帧第二投影图像的成像几何关系。

105：根据多帧第一投影图像分别对应的成像几何关系，对多帧第一投影图像进行CT三维图像重建，得到第一CT三维图像。

示例性的，确定出每帧第一投影图像对应的成像几何关系之后，则可以将每帧第一投影图像以及每帧第一投影图像对应的成像几何关系作为CT重建算法的输入数据，通过CT重建算法对多帧第一投影图像进行处理，实现对被测物体的CT图像三维重建，得到第一CT三维图像。

106：根据每组第二投影图像中的多帧第二投影图像对应的成像几何关系，对每组第二投影图像中的多帧第二投影图像进行CT三维图像重建，得到至少一幅第二CT三维图像。

同样的，确定出每组第二投影图像中的每帧第二投影图像对应的成像几何关系之后，可以将每组第二投影图像中的每帧第二投影图像以及该帧第二投影图像的成像几何关系，作为重建算法的输入数据，通过CT重建算法对每组第二投影图像中的多帧第二投影图像进行处理，实现对被测物体的CT图像三维重建，得到与每组第二投影图像对应的一幅第二CT三维图像，进而得到至少一幅第二CT三维图像。

107：根据第一CT三维图像和至少一幅第二CT三维图像，得到至少一幅三维数字减影血管造影图像。

示例性的，获取第一CT三维图像中各个体素的空间坐标，其中，第一CT三维图像中各个体素的空间坐标是在以被测物体固连标志物为坐标系下的空间坐标，即第一CT三维图像中各个体素相对于被测物体坐标系原点的位置；获取每幅第二CT三维图像中各个体素的空间坐标，其中，每幅第二CT三维图像中各个体素的空间坐标是在以被测物体固连标志物为坐标系下的空间坐标，即该幅第二CT三维图像中各个体素相对于被测物体坐标系原点的位置；根据第一CT三维图像中各个体素的空间坐标，以及每幅第二CT三维图像中各个体素的空间坐标，确定第一CT三维图像和每幅第二CT三维图像中相对应的体素，即将空间坐标相同的体素作为相对应的体素；将每幅第二CT三维图像与第一CT三维图像中相对应的体素的灰度值相减，得到至少一幅三维数字减影血管造影图像。

可以理解，当只获取了一组第二投影图像时，则可以获取出一幅三维数字减影血管造影图像，则基于该幅三维数字减影血管造影图像可以观察血管的结构、形态、尺寸等；当获取了多组第二投影图像时，则可以获取出多幅三维数字减影血管造影图像，即获取了三维数字减影血管造影图像序列，基于该三维数字减影血管造影图像序列可以显示造影剂在血管中的流动状态和医疗介入性治疗的手术状态。

可以看出，在本申请实施例中，CT成像***拍摄带有标志物的被测物体，通过标志物在图像中的像素坐标和标志物的空间坐标，用摄像测量技术高精度的得到CT成像***与被测物体之间的成像几何关系。由于坐标系建立在被测物体上，该成像几何关系自动地包含了成像***的位置姿态状态，包括振动晃动等各种误差。只要在两个拍摄投影环节中所设被测物体标志物的相对关系保持不变，则得到的被测物体CT三维成像结果是同一个坐标系的同一物体，是高精度匹配对齐的，有效地消除在各个CT运行状态期间的刚性运动或变动，提高了掩膜状态下的CT三维图像和充盈状态下的三维CT图像的对齐精度。

在本申请的一个实施方式中，当CT成像***在掩膜状态和充盈状态两种模式下运行时，在相同的位置发射X射线，也就是多个第一成像角度和多个第二成像角度一一对应相同时，则可以理解当确定出每帧第一投影图像对应的成像几何关系之后，则无需再去确定与每组第二投影图像中的每帧第二投影图像的成像几何关系，可以将每帧第一投影图像的初始成像几何关系，作为与该帧第一投影图像对应的第二投影图像的初始成像几何关系，直接参加后面的优化计算。同样的，可以只确定一组第二投影图像中的每帧第二投影图像的初始成像几何关系，无需确定其他的投影图像的成像几何关系，可以将每帧第二投影图像的初始成像几何关系作为与其对应的投影图像的初始成像几何关系，直接参加后面的优化计算，确保高精度。

参阅图4，图4本申请实施例提供的一种三维数字减影血管造影装置的功能单元组成框图。三维数字减影血管造影装置400包括：获取单元401和处理单元402，其中：

获取单元401，用于获取多帧第一投影图像，所述多帧第一投影图像是血管处于掩膜状态下，CT成像***在多个第一成像角度下向被测物体发射X射线进行成像得到，所述被测物体上设置有多个标志物，所述CT成像***在预定轨道上运动；获取至少一组第二投影图像，其中，每组第二投影图像中包括多帧第二投影图像，所述每组第二投影图像中的多帧第二投影图像是血管处于充盈状态下，所述CT成像***在多个第二成像角度下向被测物体发射X射线进行成像得到；

处理单元402，用于根据所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像对应的成像几何关系；

在一些可能的实施方式中，在根据所述第一CT三维图像和所述至少一幅第二CT三维图像，得到至少一幅三维数字减影血管造影图像方面，处理单元402，具体用于：

获取所述第一CT三维图像中各个体素的空间坐标，其中，所述第一CT三维图像中各个体素的空间坐标是所述被测物体坐标系下的空间坐标；

获取所述每幅第二CT三维图像中各个体素的空间坐标，其中，所述每幅第二CT三维图像中各个体素的空间坐标是所述被测物体坐标系下的空间坐标；

根据所述第一CT三维图像中各个体素的空间坐标，以及所述每幅第二CT三维图像中各个体素的空间坐标，确定所述第一CT三维图像和所述每幅第二CT三维图像中相对应的体素；

将所述每幅第二CT三维图像以及在所述第一CT三维图像中相对应的体素的灰度值相减，得到所述至少一幅三维数字减影血管造影图像。

在一些可能的实施方式中，在根据所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像对应的成像几何关系方面，处理单元402，具体用于：

根据所述CT成像***在所述预定轨道上的运动位置姿态，确定所述每帧第一投影图像对应的初始成像几何关系；

根据所述每帧第一投影图像的初始成像几何关系，以及所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像的成像几何关系。

在一些可能的实施方式中，当所述成像***的标称参数已知，且所述多帧第一投影图像中的k帧第一投影图像中存在未知空间坐标的标志物；在根据所述每帧第一投影图像的初始成像几何关系，以及所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像的成像几何关系方面，处理单元402，具体用于：

针对第一投影图像A，从所述多帧第一投影图像中选取第一投影图像B，其中，所述第一投影图像A为所述k帧第一投影图像中的任意一帧，所述第一投影图像A和所述第一投影图像B为不同的投影图像，所述第一投影图像A和所述第一投影图像B具有p个相同的标志物，所述p个标志物中包含所述第一投影图像A中的未知空间坐标的标志物；

根据所述第一投影图像A的初始成像几何关系、所述第一投影图像B的初始成像几何关系、所述p个标志物分别在所述第一投影图像A以及所述第一投影图像B中的像素坐标，得到所述p个标志物中的未知空间坐标的标志物的空间坐标；

根据所述每帧第一投影图像对应的初始成像几何关系、所述每帧第一投影图像中的标志物的像素坐标，以及所述每帧第一投影图像中的标志物的空间坐标进行优化处理，得到所述每帧第一投影图像对应的成像几何关系。

针对第一投影图像A，从所述多帧第一投影图像中选取第一投影图像B，其中，所述第一投影图像A为所述k帧第一投影图像中的任意一帧，所述第一投影图像A和所述第一投影图像B为不同的投影图像，所述第一投影图像A和所述第一投影图像B具有p个相同的标志物，所述p个标志物中包含未知空间坐标的标志物；

根据所述成像***的运动轨迹位置姿态，确定所述第一投影图像A和所述第一投影图像B成像时所述成像***之间的相对旋转平移矩阵；

根据所述相对旋转平移矩阵、所述p个标志物分别在所述第一投影图像A以及所述第一投影图像B中的像素坐标，得到所述p个标志物中的未知空间坐标的标志物的空间坐标；

在一些可能的实施方式中，当所述成像***的标称参数未知，且所述每帧第一投影图像中已知空间坐标的标志物的数量均少于第一阈值时，在根据所述成像***在所述预定轨道上的运动位置姿态，确定所述每帧第一投影图像对应的初始成像几何关系方面，处理单元402，具体用于：

根据所述成像***的运动位置姿态，确定第一投影图像A和第一投影图像B在成像时所述成像***之间的相对旋转平移矩阵，其中，所述第一投影图像A和所述第一投影图像B为所述多帧第一投影图像中具有w个相同的标志物的两帧投影图像，且w大于或等于第二阈值；

将所述w个相同的标志物作为w个第一标志物；

根据所述w个第一标志物分别在所述第一投影图像A和所述第一投影图像B中的像素坐标以及所述相对旋转平移矩阵，确定所述w个第一标志物的空间坐标；

根据所述w个第一标志物分别在所述第一投影图像A中的像素坐标以及所述w个第一标志物的空间坐标，确定所述第一投影图像A的初始成像几何关系；

根据所述w个第一标志物分别在所述第一投影图像B中的像素坐标以及所述w个第一标志物的空间坐标，确定所述第一投影图像B的初始成像几何关系；

根据所述w个第一标志物的空间坐标，所述第一投影图像A的初始成像几何关系或所述第一投影图像B的初始成像几何关系，确定所述每帧第一投影图像的初始成像几何关系以及所述每帧第一投影图像中的标志物的空间坐标。

在一些可能的实施方式中，在根据所述每帧第一投影图像对应的初始成像几何关系，以及所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像的成像几何关系方面，处理单元402，具体用于：

参阅图5，图5为本申请实施例提供的一种三维数字减影血管造影装置的结构示意图。如图5所示，电子设备500包括发射器501、处理器502和存储器503。它们之间通过总线504连接。存储器503用于存储计算机程序和数据，并可以将存储器503存储的数据传输给处理器502。

处理器502用于读取存储器503中的计算机程序执行以下操作：

控制发射器501向被测物体发射X射线，以获取多帧第一投影图像，所述多帧第一投影图像是血管处于掩膜状态下，CT成像***在多个第一成像角度下向被测物体发射X射线进行成像得到，所述被测物体上设置有多个标志物，所述CT成像***在预定轨道上运动；

控制发射器501向被测物体发射X射线，获取至少一组第二投影图像，其中，每组第二投影图像中包括多帧第二投影图像，所述每组第二投影图像中的多帧第二投影图像是血管处于充盈状态下，所述CT成像***在多个第二成像角度下向被测物体发射X射线进行成像得到；

其中，处理器502的具体实现功能与上述处理单元402的功能相似，不再叙述。

本申请实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行以实现如上述方法实施例中记载的任何一种三维数字减影血管造影方法的部分或全部步骤。

本申请实施例还提供一种计算机程序产品，所述计算机程序产品包括存储了计算机程序的非瞬时性计算机可读存储介质，所述计算机程序可操作来使计算机执行如上述方法实施例中记载的任何一种三维数字减影血管造影方法的部分或全部步骤。

需要说明的是，对于前述的各方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本申请并不受所描述的动作顺序的限制，因为依据本申请，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于可选实施例，所涉及的动作和模块并不一定是本申请所必须的。

在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述的部分，可以参见其他实施例的相关描述。

在本申请所提供的几个实施例中，应该理解到，所揭露的装置，可通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个***，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本申请各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件程序模块的形式实现。

所述集成的单元如果以软件程序模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储器中。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储器中，包括若干指令用以使得一台计算机设备(可为个人计算机、服务器或者网络设备等)执行本申请各个实施例所述方法的全部或部分步骤。而前述的存储器包括：U盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。

本领域普通技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储器中，存储器可以包括：闪存盘、只读存储器(英文：Read-Only Memory，简称：ROM)、随机存取器(英文：Random Access Memory，简称：RAM)、磁盘或光盘等。

以上对本申请实施例进行了详细介绍，本文中应用了具体个例对本申请的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本申请的方法及其核心思想；同时，对于本领域的一般技术人员，依据本申请的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本申请的限制。

Claims

1.一种三维数字减影血管造影方法，其特征在于，包括：

2.根据权利要求1所述的方法，其特征在于，所述根据所述第一CT三维图像和所述至少一幅第二CT三维图像，得到至少一幅三维数字减影血管造影图像，包括：

3.根据权利要求1或2所述的方法，其特征在于，所述根据所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像对应的成像几何关系，包括：

4.根据权利要求3所述的方法，其特征在于，

当所述成像***的标称参数已知，且所述多帧第一投影图像中的k帧第一投影图像中存在未知空间坐标的标志物；所述根据所述每帧第一投影图像的初始成像几何关系，以及所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像的成像几何关系，包括：

5.根据权利要求3所述的方法，其特征在于，

6.根据权利要求3所述的方法，其特征在于，

当所述成像***的标称参数未知，且所述每帧第一投影图像中已知空间坐标的标志物的数量均少于第一阈值时，所述根据所述成像***在所述预定轨道上的运动位置姿态，确定所述每帧第一投影图像对应的初始成像几何关系，包括：

将所述w个相同的标志物作为w个第一标志物；

7.根据权利要求6所述的方法，其特征在于，所述根据所述每帧第一投影图像对应的初始成像几何关系，以及所述每帧第一投影图像中的标志物的像素坐标，确定所述每帧第一投影图像的成像几何关系，包括：

8.一种三维数字减影血管造影装置，其特征在于，包括：

9.一种电子设备，其特征在于，包括：处理器和存储器，所述处理器与所述存储器相连，所述存储器用于存储计算机程序和所述多帧投影图像，所述处理器用于执行所述存储器中存储的计算机程序，以执行如权利要求1-7中任一项所述的方法。

10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序和所述多帧投影图像，所述计算机程序被处理器执行以实现如权利要求1-7中任一项所述的方法。