CN112989081B - 一种数字重建影像库的构建方法和装置 - Google Patents

一种数字重建影像库的构建方法和装置 Download PDF

Info

Publication number
CN112989081B
CN112989081B CN202110548584.6A CN202110548584A CN112989081B CN 112989081 B CN112989081 B CN 112989081B CN 202110548584 A CN202110548584 A CN 202110548584A CN 112989081 B CN112989081 B CN 112989081B
Authority
CN
China
Prior art keywords
ray
projection
image
images
drr
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.)
Expired - Fee Related
Application number
CN202110548584.6A
Other languages
English (en)
Other versions
CN112989081A (zh
Inventor
舒丽霞
李萌
蔺嫦燕
陈彧
濮欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
BEIJING INSTITUTE OF HEART LUNG AND BLOOD VESSEL DISEASES
Beijing Anzhen Hospital
Original Assignee
BEIJING INSTITUTE OF HEART LUNG AND BLOOD VESSEL DISEASES
Beijing Anzhen Hospital
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by BEIJING INSTITUTE OF HEART LUNG AND BLOOD VESSEL DISEASES, Beijing Anzhen Hospital filed Critical BEIJING INSTITUTE OF HEART LUNG AND BLOOD VESSEL DISEASES
Priority to CN202110548584.6A priority Critical patent/CN112989081B/zh
Publication of CN112989081A publication Critical patent/CN112989081A/zh
Application granted granted Critical
Publication of CN112989081B publication Critical patent/CN112989081B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/50—Information retrieval; Database structures therefor; File system structures therefor of still image data
    • G06F16/51—Indexing; Data structures therefor; Storage structures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/37—Determination of transform parameters for the alignment of images, i.e. image registration using transform domain methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10116—X-ray image
    • G06T2207/10124—Digitally reconstructed radiograph [DRR]

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种数字重建影像库的构建方法和装置，包括：获取一对术中X线图像与术前CT图像，X线图像由术中X光机采集，CT图像是重建的3D CT图像；获取X光机投影人体生成X线图像时的投影定位部分信息，以此计算得到生成X线图像时的X线投影向量场；采样X线投影向量场的未知投影参数，组合X线投影向量场与采样的未知投影参数，得到系列投影向量场；模拟生成X线图像的过程，将CT图像中的CT值转换成衰减系数，以系列投影向量场逐一投影CT图像，获得系列DRR图像，实现DRR图像库的构建。本发明显著缩小了DRR图像库，在加快建库的同时，还能改善X线与CT图像配准性能，并与现有单次快速DRR技术有效兼容。

Description

一种数字重建影像库的构建方法和装置
技术领域
本发明涉及2D/3D医学图像配准领域，尤其涉及术中X线图像与术前CT图像配准过程中数字重建影像库的构建方法和装置。
背景技术
在介入、放疗、内窥镜等手术过程中，医生借助X线照射病人生成的X线图像来实施手术，然而术中X线图像比较模糊，噪声丰富，而且一些组织如血管等，不能显示在其中，影响手术精确操作；因此，经常需要将术中2D X线图像与术前重建的3D CT图像叠加在一起，即配准2D X线图像与3D CT图像，以帮助医生观察病灶，导航手术器械，实现精准治疗。配准过程通常如下：首先，利用投影变换和数字重建影像（Digital ReconstructedRadiograph, DRR）方法，将重建的3D CT图像投影到二维空间，一次投影得到一张2D DRR图像；遍历投影空间的所有投影变换，得到一个完整数字重建影像库（简称DRR图像库）；然后，在DRR图像库中搜索，寻找与X线图像最为匹配的DRR图像；最后，叠加X线图像与最佳匹配DRR图像，实现X线图像与CT图像的配准。投影变换与DRR图像一一对应，在DRR图像库中搜索最佳匹配DRR图像，也就是搜索配准X线与CT图像的投影变换。可见，构建DRR图像库是配准X线与CT图像一个至关重要的环节。
然而，构建DRR图像库极为耗时。采用最为常用的光线追踪DRR方法，生成一张DRR图像需要几秒钟；构建一套完整的DRR图像库，则需要花费数百乃至数千个小时，令人难以承受。细密地采样投影空间，生成庞大的DRR图像库，可以从中搜索到与X线图像最佳匹配的DRR图像；如果降低投影空间采样频率，缩小DRR图像库，那么只能从中搜索到与X线图像次佳匹配的DRR图像，继而导致X线与CT图像配准精度下降。为了平衡配准精度和建库时间，现有技术主要从如下两个方向加快库的构建：1、使用硬件加速，或者改进DRR算法，以加快单次DRR投影过程，缩短单张DRR 图像的生成时间；此类技术需遍历投影空间内的全部投影变换，所生成的图像库没有缩小。2、将投影参数分为平面内参数和平面外参数，先对平面外参数采样投影，生成部分DRR图像库，再对其中每张DRR进行平面内形变处理，以替代平面内参数采样投影生成的DRR图像，从而获得完整的图像库；因为仅对平面外参数采样投影DRR，建库时间大幅缩减，但是平面外参数采样后投影时X线穿过的3D CT区域与平面内参数采样后投影时X线穿过的3D CT区域不同，生成的DRR图像内容也有所不同，前者经平面内形变处理后，并不能有效替代后者，反而损失和混淆DRR图像信息，导致配准精度降低；此类技术需间接地遍历投影空间内的全部投影变换，所生成的图像库也没有缩小。
发明内容
我们发现，术中X光机及其采集的术中X线图像dicom数据提供有生成术中X线图像时的投影定位部分信息，充分利用这些信息，可以显著缩小DRR图像库，大幅加快DRR图像库的构建，并且一定程度提升配准精度和速度。
我们将X线与CT图像的配准过程视为术中人体与术前3D CT图像的重合过程。假设术中人体与术前3D CT在空间位置上完全重合，如果已知术中投影人体生成X线图像的投影变换，以此变换投影3D CT图像，所生成的DRR图像即是与X线图像最佳匹配的DRR图像，此投影变换就是配准X线与CT图像的投影变换；如果只知道生成X线图像时的投影定位部分信息，那么结合这些已知信息，遍历投影空间内未知的投影变换，在如此构建的DRR图像库中搜索，亦可以搜寻到与X线图像最佳匹配的DRR图像，生成它的投影变换，即配准X线与CT图像的投影变换。第二种情况下，构建DRR图像库时，结合生成X线图像时的投影定位部分信息，只需遍历投影空间内未知的投影变换，相对遍历投影空间内全部的投影变换，所构建的DRR图像库明显缩小，建库时间可以大幅缩短。
有鉴于此，本发明提供了一种数字重建影像库的构建方法，包括：
获取一对术中X线图像与术前CT图像，所述术中X线图像由术中X光机采集，所述CT图像是对CT切片重建而成的3D CT图像；
获取所述X光机投影人体生成所述术中X线图像时的投影定位部分信息；依据这些已知信息，计算得到所述X光机投影生成所述术中X线图像时的X线投影向量场，由此确定投影空间内的已知部分；
在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；
模拟X线投影人体生成X线图像的过程，将所述CT图像中的CT值转换成衰减系数，以所述系列投影向量场逐一投影CT图像，获得系列DRR 图像，实现DRR 图像库的构建。
优选的，所述X光机为C型臂X光机。
基于上述构建方法，本发明还提供了一种数字重建影像库的构建装置，包括：
数据获取单元，用于：获取一对术中X线图像与术前CT图像，所述术中X线图像由术中X光机采集，所述CT图像是对CT切片重建而成的3D CT图像；
X线投影向量场计算单元，用于：获取所述X光机投影人体生成所述术中X线图像时的投影定位部分信息，依据这些已知信息，计算得到所述X光机投影生成所述术中X线图像时的X线投影向量场，由此确定投影空间内的已知部分；
DRR图像库生成单元，用于：在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；模拟X线投影人体生成X线图像的过程，将所述CT图像中的CT值转换成衰减系数，以所述系列投影向量场逐一投影CT图像，获得系列DRR 图像，实现DRR图像库的构建。
优选的，所述X光机为C型臂X光机。
本发明充分利用能获取的X光机生成术中X线图像的投影定位部分信息，只需遍历投影空间内未知的投影变换，使得生成的DRR图像库大幅缩小，建库时间显著加快。如前所述，现有两类技术所生成的图像库并没有缩小，相对而言，本发明具有以下有益效果：
（1）本发明显著缩小了图像库，在加快DRR 建库的同时，还能加快X线与CT图像配准速度，并提升配准精度。
首先，现有技术需直接或间接遍历投影空间内的全部投影变换，图像库中图像数量并没有减少，而本发明仅遍历投影空间内未知的投影变换，图像库中图像数量大幅减少，也就是说，相对现有技术，本发明大幅减少了配准时需要搜索的DRR 图像数量，因此可以明显加快配准速度。
其次，现有技术需直接或间接遍历投影空间内的全部投影变换，因而投影空间内的所有投影参数都会对配准精度产生影响；相对而言，本发明利用投影定位已知信息得到的投影参数不会产生配准误差，仅余下未知的投影参数影响配准精度，因此，虽然缩小了DRR图像库，反而一定程度可以提升配准精度。
再者，相较于现有基于平面内外参数的建库技术，本发明提升配准精度的优势更为突出。前者通过平面内形变处理平面外参数采样投影生成的DRR图像，来替代平面内参数采样投影生成的DRR图像，以加快图像库的构建，但是如前所述，这种替代损失和混淆了DRR图像信息，导致配准精度降低；本发明结合投影定位已知信息，将投影空间缩小在未知参数空间内，然后逐一对其中每一个投影变换进行DRR，所生成的DRR图像没有信息损失和混淆。
（2）本发明可以有效兼容现有的单次快速DRR技术，在保障X线与CT图像配准性能得到一定改善的前提下，进一步加快DRR图像库的构建。
本发明着眼于缩小DRR 图像库，现有单次快速DRR技术着眼于缩短单张DRR 图像的生成时间，两者互为补充，彼此有效结合，可以在保障X线与CT图像配准性能得到一定改善的前提下，进一步加快DRR图像库构建。
相对而言，现有基于平面内外参数的建库技术虽然也可与现有单次快速DRR技术相结合，但并不能保障X线与CT图像配准性能得到改善。此类技术最终获得的图像库中图像数量并没有减少，也就是说，配准时需要搜索的DRR图像数量并未减小，配准速度得不到提升。而且，如前所述，平面内形变处理平面外参数采样投影生成的DRR图像，并不能有效替代平面内参数采样投影生成的DRR图像，反而损失和混淆DRR图像信息，导致配准精度降低。
附图说明
为了清楚地理解本发明的技术方案，下面对描述本发明具体实施方式时用到的附图进行简要说明。显而易见，这些附图仅是本发明的一部分附图，本领域普通技术人员在不付出创造性劳动的前提下，还可以获得其它附图。
图1是本发明实施例1提供的一种数字重建影像库的构建方法流程示意图；
图2是本发明实施例1提供的X线投影向量场计算流程示意图；
图3是本发明实施例1提供的X光机投影定位部分信息示意图
图4是本发明实施例1提供的DRR图像库生成流程示意图；
图5是本发明实施例1提供的系列X线投影向量场计算流程示意图
图6是本发明实施例2提供的一种数字重建影像库的构建装置结构示意图；
图7是采用本发明实施例3和传统DRR方法分别配准一组术中X线与术前CT图像对的配准效果图。
具体实施方式
为使本发明的目的、技术手段和有益效果更加清楚完整，下面结合附图对本发明的具体实施方式进行描述。
实施例1
本发明实施例1提供了一种数字重建影像库的构建方法，其流程示意图如图1所示，包括以下步骤：
S1、获取一对术中X线图像与术前CT图像，所述术中X线图像由术中X光机采集，所述CT图像是对CT切片重建而成的3D CT图像。
S2、获取所述X光机投影人体生成所述术中X线图像时的投影定位部分信息，依据这些信息，计算得到所述X光机投影生成所述术中X线图像时的X线投影向量场，由此确定生成X线图像时投影空间内的已知部分。
图2显示了步骤S2的实现过程，具体包括以下步骤：
S21、获取投影定位部分信息：
如图3所示，从所述术中X线图像的dicom数据及所述术中X光机说明书中，读取X光机生成所述术中X线图像时的部分参数，这些参数包括：
S211、等中心点距地面高度，从X光机说明书中获取；
S212、治疗床距地面高度，先从X光机说明书中获取治疗床高度初始值，再从X线图像dicom数据标签中读取治疗床高度变化值，两者之和即治疗床距地面高度；
S213、点光源S沿X线投影方向到等中心点C的距离，从X光机说明书中获取；
S214、点光源S沿X线投影方向到平板探测器中心点Dc的距离，从X线图像dicom数据标签中读取；
以及从X线图像dicom数据标签中读取的机架带动X光机围绕等中心点C绕三个轴向旋转的旋转角 α、β、γ，和从X光机说明书中获取的平板探测器尺寸和X线图像分辨率。
上述多为大型C型臂X光机定位参数的获取方式；对于中小型C型臂X光机，等中心点距地面高度会在手术中进行调整，则多从X线图像dicom数据标签中读取。
S22、确定等中心点C坐标(xc, yc, zc)中的x c ：
如图3所示，坐标系以垂直地面向上为x轴方向，以人体头部指向足端为z轴方向，采用右手法则确定y轴方向，并令X线照射区域人体最上端平行于治疗床的切面为平面yoz，S215代表X线照射区域人体厚度；等中心点C的坐标(xc, yc, zc)中，x c = H c - (H t + H CT)，其中H c表示等中心点C距地面高度，H t为治疗床距地面高度，H CT表示CT切片的高度，也即X线照射区域人体厚度；yc和zc为未知参数；
X光机在投影人体时，给出了人体所在高度信息H c 、H t 、H CT，没有给出人体在yoz平面的确切位置，所以只能确定等中心点坐标中的x c ，而不能确定yc和zc；下面在S23至S26的计算步骤中，先假设参数yc和zc已知。
S23、计算X光机围绕等中心点绕三个轴向旋转的变换矩阵：
首先平移坐标系，将原点移至等中心点C，计算平移矩阵T t ；然后计算X光机依次绕x、y、z三个轴向旋转α、β、γ角后所得的旋转矩阵R；最后再平移坐标系，将原点从等中心点C移回原坐标系位置，平移矩阵为T t 的逆矩阵T t -1 ；变换矩阵T的计算公式如下：
Figure 30420DEST_PATH_IMAGE001
其中
Figure 160050DEST_PATH_IMAGE002
，
Figure 505580DEST_PATH_IMAGE003
即x轴，
Figure 328043DEST_PATH_IMAGE004
为绕
Figure 126235DEST_PATH_IMAGE003
旋转α角所对应的旋转矩阵；向量
Figure 121872DEST_PATH_IMAGE005
，是坐标系绕x轴旋转α角后y轴所对应向量，
Figure 259593DEST_PATH_IMAGE006
表示绕
Figure 315273DEST_PATH_IMAGE007
旋转β角所对应的旋转矩阵；向量
Figure 804024DEST_PATH_IMAGE008
，是坐标系先绕x轴旋转α角、再绕y轴旋转β角后z轴所对应向量，
Figure 337773DEST_PATH_IMAGE009
表示绕
Figure 595579DEST_PATH_IMAGE010
旋转γ角所对应的旋转矩阵。旋转矩阵
Figure 822161DEST_PATH_IMAGE011
依照如下公式计算：
Figure 532628DEST_PATH_IMAGE012
这里
Figure 870069DEST_PATH_IMAGE013
表示绕向量
Figure 247960DEST_PATH_IMAGE014
旋转θ角的旋转矩阵，将旋转矩阵
Figure 379864DEST_PATH_IMAGE011
的参数与
Figure 843207DEST_PATH_IMAGE013
的参数一一对应，即可计算得到各个矩阵值；
无论实际操作中X光机绕三个轴向旋转的先后次序如何，最终得到的旋转矩阵与按上述方法计算得到的旋转矩阵R相同。
S24、计算X光机点光源的坐标：
点光源S坐标(xs, ys, zs)的计算公式如下：
Figure 718759DEST_PATH_IMAGE015
其中，
Figure 951157DEST_PATH_IMAGE016
表示点光源S沿X线投影方向到等中心点C的距离，默认点光源S初始位置在治疗床下方，X线
Figure 253962DEST_PATH_IMAGE017
初始位置垂直于治疗床。
S25、计算平板探测器所有点的坐标：
平板探测器所有点Dmn坐标
Figure 204601DEST_PATH_IMAGE018
的计算方法如下：
首先根据公式（6）计算平板探测器中心点Dc的坐标：
Figure 883844DEST_PATH_IMAGE019
其中
Figure 970749DEST_PATH_IMAGE020
表示点光源S沿X线投影方向到平板探测器中心点Dc的距离，然后根据公式（7）计算平板探测器所有点Dmn的坐标：
Figure 721753DEST_PATH_IMAGE022
其中，
Figure 894109DEST_PATH_IMAGE023
表示点Dmn的初始位置坐标，
Figure 314726DEST_PATH_IMAGE024
，
Figure 584033DEST_PATH_IMAGE025
，SL*SW为平板探测器尺寸，RL*RW是X线图像分辨率；当
Figure 166324DEST_PATH_IMAGE026
，
Figure 888292DEST_PATH_IMAGE027
时，
Figure 847021DEST_PATH_IMAGE028
,由此推算出
Figure 970835DEST_PATH_IMAGE029
；将
Figure 989606DEST_PATH_IMAGE029
代入式（7），得到平板探测器所有点Dmn的坐标
Figure 198871DEST_PATH_IMAGE018
。
S26、计算X光机点光源投影到平板探测器所有点的X线投影向量场：
X线投影向量场
Figure 695711DEST_PATH_IMAGE030
的计算公式如下：
Figure 877294DEST_PATH_IMAGE031
其中
Figure 129284DEST_PATH_IMAGE032
。
在上述S22到S26计算生成X线图像投影向量场的全部过程中，参数yc和zc未知；参数
Figure 763527DEST_PATH_IMAGE033
在步骤S21中已经获得，为已知参数；参数H CT = RCT-h*Ph，RCT-h和Ph分别表示CT切片的高度分辨率和像素高度，均从CT切片的dicom数据中读取，因而参数H CT亦为已知参数。
假设已知参数yc和zc，通过上述步骤，便可准确计算得到生成X线图像时的投影向量场，以此投影向量场投影3D CT图像，就可得到与X线图像最佳匹配的DRR图像；换句话说，一旦确定参数yc和zc，就可以实现X线与CT图像配准。然而，如前所述，实际上参数yc和zc未知，为了实现X线与CT图像配准，就需要确定参数yc和zc的值。考虑到DRR图像与投影变换一一对应，如果寻找到与X线图像最佳匹配的DRR图像，其对应的投影变换，就是配准X线与CT图像的投影变换，也即投影生成X线图像的投影向量场；知道投影向量场，便可推出参数yc和zc的值。所以，为了确定参数yc和zc，实现X线与CT图像配准，需要在X线照射区域遍历（yc,zc）的所有组合，得到所有可能的投影向量场，逐一对3D CT图像进行投影DRR，获得DRR图像库；库中与X线图像最佳匹配的DRR图像所对应的（yc, zc）组合，就是参数yc和zc的值。
因此，经步骤S2确定投影空间内的已知部分后，下一步，就是遍历（yc, zc）的所有可能组合，生成DRR图像库。
S3、在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；模拟X线投影人体生成X线图像的过程，将所述CT图像中的CT值转换成衰减系数，以所述系列投影向量场逐一投影CT图像，获得系列DRR 图像，实现DRR 图像库的构建。
图4显示了步骤S3的具体实现过程，包括以下三个步骤：
S31、在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换。
如图5所示，步骤S31又包括如下步骤：
S311、指定坐标系原点：
令CT dicom数据中第一个被读取的体素点为坐标系原点；
S312、采样投影向量场的未知参数：
在投影向量场
Figure 126376DEST_PATH_IMAGE030
未知参数yc和zc所在的yoz平面内，高频率、等间隔、全覆盖地采样两个参数，得到系列（yci, zcj）值；其中
Figure 896886DEST_PATH_IMAGE034
，A表示3D CT在yoz平面覆盖的区域，也即X线在yoz平面的照射区域；i = 0,1,…,I-1，j = 0,1,…,J-1，I和J分别表示在y轴和z轴上的采样点数；
S313、获得系列投影向量场：
将采样值（yci ,zcj）逐一代入公式（8），获得系列投影向量场
Figure 319777DEST_PATH_IMAGE035
，其中i和j对应采样值（yci ,zcj），i = 0,1,…,I-1，j = 0,1,…,J-1；m和n对应平板探测器点Dmn，m = 0,1,…，RL-1，n = 0,1,…，RW-1；例如，采样值(yc0 ,zc0) 对应得到X光机点光源S到平板探测器所有点Dmn的投影向量场
Figure 441316DEST_PATH_IMAGE036
，采样值(yc2 ,zc3) 对应得到点光源S到探测器所有点Dmn的投影向量场
Figure 342276DEST_PATH_IMAGE037
；遍历(yci ,zcj)所有组合，也就遍历了投影空间内未知的所有投影变换。
S32、将CT图像中的CT值转化成X射线衰减系数：
为了模拟X线投影人体生成术中X线图像的过程，首先将所述3D CT图像中每一个体素的CT值转化为相应组织对X线的衰减系数，转化公式如下：
Figure 232872DEST_PATH_IMAGE038
其中，k表示X线穿过人体的第k类组织，
Figure 826664DEST_PATH_IMAGE039
表示第k类组织对X线的衰减系数，
Figure 435500DEST_PATH_IMAGE040
表示水对X线的衰减系数，HU表示CT值。
S33、以所述系列投影向量场逐一投影CT图像到平板探测器，获得系列DRR 图像，实现DRR图像库的构建：
模拟X线投影人体生成术中X线图像的过程，X线从点光源S出发，以所述系列投影向量场
Figure 343413DEST_PATH_IMAGE035
（i = 0,1,…,I-1，j = 0,1,…,J-1）逐一投影3D CT图像到平板探测器所有点Dmn（m = 0,1,…，RL-1，n = 0,1,…，RW-1），采用光线追踪法生成系列DRR 图像；点Dmn的灰度值DRR(m,n)计算如下：
Figure 150832DEST_PATH_IMAGE041
其中，I 0表示X线光源强度，可以指定为1；k表示X线穿过人体的第k类组织，
Figure 853209DEST_PATH_IMAGE039
表示第k类组织的衰减系数；
Figure 746079DEST_PATH_IMAGE042
表示X线
Figure 926524DEST_PATH_IMAGE035
穿过第k类组织的长度，通过计算X线穿入和穿出第k类组织的两点间距得到；
投影向量场
Figure 854029DEST_PATH_IMAGE035
遍历m和n，得到对应采样参数（yci ,zcj）投影生成的DRR图像；例如，投影向量场
Figure 727307DEST_PATH_IMAGE036
遍历m和n，得到采样参数(yc0 ,zc0) 投影生成的DRR图像；投影向量场
Figure 107473DEST_PATH_IMAGE037
遍历m和n，得到采样参数(yc2 ,zc3) 投影生成的DRR图像；以系列投影向量场
Figure 91609DEST_PATH_IMAGE035
（i =0,1,…,I-1，j = 0,1,…,J-1）逐一投影3D CT图像，获得系列DRR 图像，实现DRR图像库的构建。
以上即为本发明实施例1提供的一种数字重建影像库构建方法的具体实施方式。相对于现有技术，此方法具有如下优势：（1）显著缩小了图像库，在加快DRR 建库的同时，X线与CT图像配准搜索的DRR 图像量亦大幅减少，能明显加快配准速度；而且，配准不受投影空间内已知参数影响，仅受投影空间内未知参数影响，因而一定程度可以提升配准精度；相较于存在信息损失和混淆的平面内外参数建库技术，本发明提升配准精度的优势更为突出。（2）本发明着眼于缩小DRR 图像库，现有单次快速DRR技术着眼于缩短单张DRR 图像的生成时间，两者互为补充，彼此有效结合，可以在保障X线与CT图像配准性能得到一定改善的前提下，进一步加快DRR图像库构建。
实施例2
基于实施例1提供的一种数字重建影像库构建方法，本发明实施例2还提供了一种数字重建影像库的构建装置，其结构示意图如图6所示，该装置包括以下单元：
数据获取单元41，用于：获取一对术中X线图像与术前CT图像，所述术中X线图像由术中X光机采集，所述CT图像是对CT切片重建而成的3D CT图像；
X线投影向量场计算单元42，用于：获取所述X光机投影人体生成所述术中X线图像时的投影定位部分信息，依据这些信息，计算得到所述X光机投影生成所述术中X线图像时的X线投影向量场，由此确定生成X线图像时投影空间内的已知部分；
DRR图像库生成单元43，用于：在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；模拟X线投影人体生成X线图像的过程，将所述CT图像中的CT值转换成衰减系数，以所述系列投影向量场逐一投影CT图像，获得系列DRR 图像，实现DRR 图像库的构建。
相对于现有技术，实施例2提供的一种数字重建影像库构建装置具有如下优势：（1）显著缩小了图像库，在加快DRR 建库的同时，X线与CT图像配准搜索的DRR 图像量亦大幅减少，能明显加快配准速度；而且，配准不受投影空间内已知参数影响，仅受投影空间内未知参数影响，因而能一定程度提升配准精度；相较于存在信息损失和混淆的平面内外参数建库技术，本发明提升配准精度的优势更为突出。（2）本发明着眼于缩小DRR 图像库，现有单次快速DRR技术着眼于缩短单张DRR 图像的生成时间，两者互为补充，彼此有效结合，可以在保障X线与CT图像配准性能得到一定改善的前提下，进一步加快DRR图像库构建。
所述X线投影向量场计算单元42包括：
获取投影信息子单元421，用于：从X光机说明书中获取等中心点距地面高度、点光源S沿X线投影方向到等中心点C的距离、平板探测器尺寸和X线图像分辨率；从X线图像dicom数据标签中读取点光源S沿X线投影方向到平板探测器中心点Dc的距离、机架带动X光机围绕等中心点C绕三个轴向旋转的旋转角α、β、γ；以及先从X光机说明书中获取治疗床高度初始值，再从X线图像dicom数据标签中读取治疗床高度变化值，两者之和即治疗床距地面高度。
上述多为大型C型臂X光机定位参数的获取方式；对于中小型C型臂X光机，等中心点距地面高度会在手术中进行调整，则多从X线图像dicom数据标签中读取。
计算X线投影向量场子单元422，用于：计算生成所述术中X线图像时的X线投影向量场，具体又包括：确定等中心点子单元4221、计算变换矩阵子单元4222、计算点光源坐标子单元4223、计算探测器坐标子单元4224，和计算投影向量场子单元4225。
所述确定等中心点子单元4221，用于：确定等中心点C坐标(xc, yc, zc)中的x c ；
坐标系以垂直地面向上为x轴方向，以人体头部指向足端为z轴方向，采用右手法则确定y轴方向，并令X线照射区域人体最上端平行于治疗床的切面为平面yoz；等中心点C的坐标(xc, yc, zc)中，x c = H c - (H t + H CT)，其中H c表示等中心点C距地面高度，H t为治疗床距地面高度，H CT表示CT切片的高度，也即X线照射区域人体厚度；yc和zc为未知参数。
所述计算变换矩阵子单元4222，用于：计算X光机以等中心点为中心、围绕三个轴向旋转的变换矩阵T；
首先平移坐标系，将原点移至等中心点C，计算平移矩阵T t ；然后计算X光机依次绕x、y、z三个轴向旋转α、β、γ角后所得的旋转矩阵R；最后再平移坐标系，将原点从等中心点C移回原坐标系位置，平移矩阵为T t 的逆矩阵T t -1 ；变换矩阵T的计算公式如下：
Figure 873621DEST_PATH_IMAGE001
其中
Figure 183379DEST_PATH_IMAGE002
，
Figure 988524DEST_PATH_IMAGE003
即x轴，
Figure 838669DEST_PATH_IMAGE004
为绕
Figure 412869DEST_PATH_IMAGE003
旋转α角所对应的旋转矩阵；向量
Figure 690267DEST_PATH_IMAGE005
，是坐标系绕x轴旋转α角后y轴所对应向量，
Figure 717129DEST_PATH_IMAGE006
表示绕
Figure 902122DEST_PATH_IMAGE007
旋转β角所对应的旋转矩阵；向量
Figure 330830DEST_PATH_IMAGE008
，是坐标系先绕x轴旋转α角、再绕y轴旋转β角后z轴所对应向量，
Figure 779129DEST_PATH_IMAGE009
表示绕
Figure 558866DEST_PATH_IMAGE010
旋转γ角所对应的旋转矩阵。旋转矩阵
Figure 750813DEST_PATH_IMAGE011
依照如下公式计算：
Figure 34027DEST_PATH_IMAGE012
这里
Figure 122068DEST_PATH_IMAGE013
表示绕向量
Figure 185839DEST_PATH_IMAGE014
旋转θ角的旋转矩阵，将旋转矩阵
Figure 853581DEST_PATH_IMAGE011
的参数与
Figure 319197DEST_PATH_IMAGE013
的参数一一对应，即可计算得到各个矩阵值；
无论实际操作中X光机绕三个轴向旋转的先后次序如何，最终得到的旋转矩阵与按上述方法计算得到的旋转矩阵R相同。
所述计算点光源坐标子单元4223，用于：计算X光机点光源S的坐标(xs, ys, zs)，计算公式如下：
Figure 312561DEST_PATH_IMAGE015
其中，
Figure 598049DEST_PATH_IMAGE016
表示点光源S沿X线投影方向到等中心点C的距离，默认点光源S初始位置在治疗床下方，X线
Figure 335061DEST_PATH_IMAGE017
初始位置垂直于治疗床。
所述计算探测器坐标子单元4224，用于：计算平板探测器所有点Dmn的坐标
Figure 655184DEST_PATH_IMAGE018
，计算方式如下：
首先根据公式（6）计算平板探测器中心点Dc的坐标：
Figure 553870DEST_PATH_IMAGE019
其中
Figure 592233DEST_PATH_IMAGE020
表示点光源S沿X线投影方向到平板探测器中心点Dc的距离，然后根据公式（7）计算平板探测器所有点Dmn的坐标：
Figure 132936DEST_PATH_IMAGE043
其中，
Figure 245248DEST_PATH_IMAGE023
表示点Dmn的初始位置坐标，
Figure 642731DEST_PATH_IMAGE024
，
Figure 840494DEST_PATH_IMAGE025
，SL*SW为平板探测器尺寸，RL*RW是X线图像分辨率；当
Figure 981626DEST_PATH_IMAGE026
，
Figure 214024DEST_PATH_IMAGE027
时，
Figure 516829DEST_PATH_IMAGE028
,由此推算出
Figure 467468DEST_PATH_IMAGE029
；将
Figure 146711DEST_PATH_IMAGE029
代入式（7），得到平板探测器所有点Dmn的坐标
Figure 499195DEST_PATH_IMAGE018
。
计算投影向量场子单元4225，用于：计算X光机点光源S投影到平板探测器所有点Dmn的X线投影向量场
Figure 707322DEST_PATH_IMAGE030
，计算公式如下：
Figure 145257DEST_PATH_IMAGE031
其中
Figure 565874DEST_PATH_IMAGE032
。
在上述计算X线投影向量场子单元422的全部过程中，参数yc和zc未知；参数
Figure 569602DEST_PATH_IMAGE033
在子单元421中已经获得，为已知参数；参数H CT = RCT-h*Ph，RCT-h和Ph分别表示CT切片的高度分辨率和像素高度，均从CT切片的dicom数据中读取，因而参数H CT亦为已知参数。
为了确定参数yc和zc，实现X线与CT图像配准，需要在X线照射区域遍历（yc, zc）的所有组合，得到所有可能的投影向量场，逐一对3D CT图像进行投影DRR，获得DRR图像库；库中与X线图像最佳匹配的DRR图像所对应的（yc, zc）组合，就是参数yc和zc的值。因此，确定投影空间内的已知部分后，下一单元，是遍历（yc, zc）的所有可能组合，生成DRR图像库。
所述DRR图像库生成单元43包括：系列投影向量场计算子单元431，CT图像预处理子单元432，和生成DRR图像库子单元433。
所述系列投影向量场计算子单元431，用于：在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换，具体又包括：
确定坐标系原点子单元4311，用于：令CT dicom数据中第一个被读取的体素点为坐标系原点；
采样未知参数子单元4312，用于：在投影向量场
Figure 417472DEST_PATH_IMAGE030
未知参数yc和zc所在的yoz平面内，高频率、等间隔、全覆盖地采样两个参数，得到系列（yci, zcj）值；其中
Figure 139441DEST_PATH_IMAGE034
，A表示3D CT在yoz平面覆盖的区域，也即X线在yoz平面的照射区域；i = 0,1,…,I-1，j = 0,1,…,J-1，I和J分别表示在y轴和z轴上的采样点数；
计算系列投影向量场子单元4313，用于：将采样值（yci ,zcj）逐一代入公式（8），获得系列投影向量场
Figure 98169DEST_PATH_IMAGE035
，其中i和j对应采样值（yci ,zcj），i = 0,1,…,I-1，j = 0,1,…,J-1；m和n对应平板探测器点Dmn，m = 0,1,…，RL-1，n = 0,1,…，RW-1；例如，采样值(yc0 ,zc0)对应得到X光机点光源S到平板探测器所有点Dmn的投影向量场
Figure 221983DEST_PATH_IMAGE036
，采样值(yc2 ,zc3) 对应得到点光源S到探测器所有点Dmn的投影向量场
Figure 240755DEST_PATH_IMAGE037
；遍历(yci ,zcj)所有组合，也就遍历了投影空间内未知的所有投影变换。
所述CT图像预处理子单元432，用于：为了模拟X线投影人体生成术中X线图像的过程，首先将所述3D CT图像中每一个体素的CT值转化为相应组织对X线的衰减系数，转化公式如下：
Figure 184440DEST_PATH_IMAGE038
其中，k表示X线穿过人体的第k类组织，
Figure 946860DEST_PATH_IMAGE039
表示第k类组织对X线的衰减系数，
Figure 925180DEST_PATH_IMAGE040
表示水对X线的衰减系数，HU表示CT值。
所述生成DRR图像库子单元433，用于：模拟X线投影人体生成术中X线图像的过程，X线从点光源S出发，以所述系列投影向量场
Figure 380432DEST_PATH_IMAGE035
（i = 0,1,…,I-1，j = 0,1,…,J-1）逐一投影3D CT图像到平板探测器所有点Dmn（m = 0,1,…，RL-1，n = 0,1,…，RW-1），采用光线追踪法生成系列DRR 图像；点Dmn的灰度值DRR(m,n)计算如下：
Figure 14676DEST_PATH_IMAGE044
其中，I 0表示X线光源强度，可以指定为1；k表示X线穿过人体的第k类组织，
Figure 377524DEST_PATH_IMAGE039
表示第k类组织的衰减系数；
Figure 148034DEST_PATH_IMAGE042
表示X线
Figure 570925DEST_PATH_IMAGE035
穿过第k类组织的长度，通过计算X线穿入和穿出第k类组织的两点间距得到；
投影向量场
Figure 426885DEST_PATH_IMAGE035
遍历m和n，得到对应采样参数（yci ,zcj）投影生成的DRR图像；例如，投影向量场
Figure 593425DEST_PATH_IMAGE036
遍历m和n，得到采样参数(yc0 ,zc0) 投影生成的DRR图像；投影向量场
Figure 484020DEST_PATH_IMAGE037
遍历m和n，得到采样参数(yc2 ,zc3) 投影生成的DRR图像；以系列投影向量场
Figure 77813DEST_PATH_IMAGE035
（i =0,1,…,I-1，j = 0,1,…,J-1）逐一投影3D CT图像，获得系列DRR 图像，实现DRR图像库的构建。
实施例3
采用实施例1和2提供的数字重建影像库构建方法和装置构建DRR图像库，以互信息为相似性度量，Powell作为优化搜索，本发明实施例3配准了行胸主动脉腔内修复术患者的一组术前CT血管造影与术中X线图像对；采用传统6自由度光线追踪DRR方法构建DRR图像库，以互信息为相似性度量，Powell作为优化搜索，我们也配准上述CT血管造影与X线图像对；我们将两种方法的实验结果进行了比较。两种配准方法唯一差异在于DRR图像库构建方法不同：实施例3采用的是实施例1和2提供的数字重建影像库构建方法和装置，而后者采用的是传统6自由度光线追踪DRR图像库构建方法。
我们在一台配置为intel i9-9820X CPU和112G内存的电脑上进行实验，实验环境为matlab。两种DRR图像库构建方法对平移的采样间隔均为3mm，对角度的采样间隔为1°。待配准的CT血管造影与X线图像对中，术前 CT切片分辨率为512*512，像素大小为0.6836mm*0.6836mm，H CT为 350mm；术中X线图像由GE Innova4100-IQ C型臂采集，从说明书及dicom数据中获得的X光机生成术中X线图像时的投影部分信息如表1所示。
表1. X光机生成术中X线图像时的投影部分信息
Figure DEST_PATH_IMAGE045
表2. 实验结果比较
Figure DEST_PATH_IMAGE046
表2对两种配准方法的实验结果进行了定量比较，图7则对两种方法的配准效果进行目视比较。图7上、中、下三张子图分别对应待配准的术中X线图像、传统DRR方法配准图像对后DRR与X线图像的叠加效果图（叠加后的DRR与X线图像以棋盘格交替显示）、以及实施例3配准图像对后DRR与X线图像的叠加效果图。
从表2可以看出，在配准上述CT血管造影与X线图像对时，相对于传统6自由度光线追踪DRR方法，实施例3采用本发明提出的数字重建影像库构建方法和装置，使得（1）DRR图像库中的图像数量大幅减少，建库时间从上千个小时骤减至不到三个小时；（2）配准速度加快了三倍多；（3）配准后DRR图像与X线图像的相似度得到一定程度提高。目视比较图7中图以及下图中箭头所指部位，特别观察DRR图像中肋骨边缘与X线图像中肋骨边缘的吻合程度，也可以看到相对于传统DRR方法，实施例3采用本发明提出的数字重建影像库构建方法和装置，能使肋骨部位的吻合程度得到提升，也即使得配准精度得到改善。
以上所述仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制。虽然本发明以较佳实施例揭露如上，然而并非用以限定本发明。任何熟悉本领域的技术人员，在不脱离本发明技术方案范围情况下，都可利用上述揭示的方法和技术内容对本发明技术方案做出许多可能的变动和修饰，或修改为等同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。

Claims (8)

1.一种数字重建影像库的构建方法，其特征在于，包括：
获取一对术中X线图像与术前CT图像，所述术中X线图像由术中X光机采集，所述CT图像是对CT切片重建而成的3D CT图像；
获取所述X光机投影人体生成所述术中X线图像时的投影定位部分信息；依据这些已知信息，计算得到所述X光机投影生成所述术中X线图像时的X线投影向量场，由此确定投影空间内的已知部分；
在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；
模拟X线投影人体生成X线图像的过程，将所述CT图像中的CT值转换成衰减系数，以所述系列投影向量场逐一投影CT图像，获得系列数字重建影像，即系列DRR 图像，实现DRR 图像库的构建。
2.根据权利要求1所述的方法，其特征在于，所述X光机为C型臂X光机。
3.根据权利要求2所述的方法，其特征在于，获取所述X光机投影人体生成所述术中X线图像时的投影定位部分信息；依据这些已知信息，计算得到所述X光机投影生成所述术中X线图像时的X线投影向量场，由此确定投影空间内的已知部分，具体包括：
步骤A、获取投影定位部分信息；
从所述术中X光机说明书及术中X线图像的dicom数据中，读取X光机生成所述术中X线图像时的部分参数，包括：等中心点距地面高度、治疗床距地面高度、点光源沿X线投影方向到等中心点的距离、点光源沿X线投影方向到平板探测器中心点的距离、机架带动X光机围绕等中心点绕三个轴向旋转的旋转角、平板探测器尺寸、X线图像分辨率；
步骤B、确定等中心点C坐标(xc, yc, zc)中的x c ；
坐标系以垂直地面向上为x轴方向，以人体头部指向足端为z轴方向，采用右手法则确定y轴方向，并令X线照射区域人体最上端平行于治疗床的切面为平面yoz；等中心点C的坐标(xc, yc, zc)中，x c = H c- (H t+ H CT)，其中H c表示等中心点C距地面高度，H t为治疗床距地面高度，H CT表示CT切片的高度，也即X线照射区域人体厚度（S215）；yc和zc为未知参数；
步骤C、计算X光机围绕等中心点绕三个轴向旋转的变换矩阵T；
首先平移坐标系，将原点移至等中心点C，计算平移矩阵T t ；然后计算X光机依次绕x、y、z三个轴向旋转α、β、γ角后所得的旋转矩阵R；最后再平移坐标系，将原点从等中心点C移回原坐标系位置，平移矩阵为T t 的逆矩阵T t -1 ；变换矩阵T的计算公式如下：
Figure 857691DEST_PATH_IMAGE001
Figure 688244DEST_PATH_IMAGE002
其中
Figure 744056DEST_PATH_IMAGE003
，
Figure 343665DEST_PATH_IMAGE004
即x轴，
Figure 216943DEST_PATH_IMAGE005
为绕
Figure 659425DEST_PATH_IMAGE004
旋转α角所对应的旋转矩阵；向量
Figure 643562DEST_PATH_IMAGE006
，是坐标系绕x轴旋转α角后y轴所对应向量，
Figure 363256DEST_PATH_IMAGE007
表示绕
Figure 407436DEST_PATH_IMAGE008
旋转β角所对应的旋转矩阵；向量
Figure 87947DEST_PATH_IMAGE009
，是坐标系先绕x轴旋转α角，再绕y轴旋转β角后z轴所对应向量，
Figure 875774DEST_PATH_IMAGE010
表示绕
Figure 449975DEST_PATH_IMAGE011
旋转γ角所对应的旋转矩阵；旋转矩阵
Figure 665056DEST_PATH_IMAGE012
依照如下公式计算：
Figure 816551DEST_PATH_IMAGE013
这里
Figure 673649DEST_PATH_IMAGE014
表示绕向量
Figure 102356DEST_PATH_IMAGE015
旋转θ角的旋转矩阵，将旋转矩阵
Figure 488338DEST_PATH_IMAGE012
的参数与
Figure 268075DEST_PATH_IMAGE014
的参数一一对应，即可计算得到各个矩阵值；
无论实际操作中X光机绕三个轴向旋转的先后次序如何，最终得到的旋转矩阵与按上述方法计算得到的旋转矩阵R相同；
步骤D、计算X光机点光源S的坐标(xs, ys, zs)，计算公式如下：
Figure 273072DEST_PATH_IMAGE016
其中，
Figure 556285DEST_PATH_IMAGE017
表示点光源S沿X线投影方向到等中心点C的距离，默认点光源S初始位置在治疗床下方，X线
Figure 378748DEST_PATH_IMAGE018
初始位置垂直于治疗床；
步骤E、计算平板探测器所有点Dmn的坐标
Figure 380202DEST_PATH_IMAGE019
，计算方法如下：
首先根据公式（6）计算平板探测器中心点Dc的坐标：
Figure 438157DEST_PATH_IMAGE020
其中
Figure 575877DEST_PATH_IMAGE021
表示点光源S沿X线投影方向到平板探测器中心点Dc的距离，然后根据公式（7）计算平板探测器所有点Dmn的坐标：
Figure 569241DEST_PATH_IMAGE022
其中，
Figure 323570DEST_PATH_IMAGE023
表示点Dmn的初始位置坐标，
Figure 795003DEST_PATH_IMAGE024
，
Figure 925245DEST_PATH_IMAGE025
，SL*SW为平板探测器尺寸，RL*RW是X线图像分辨率；当
Figure 89510DEST_PATH_IMAGE026
，
Figure 799977DEST_PATH_IMAGE027
时，
Figure 75101DEST_PATH_IMAGE028
，由此推算出
Figure 577626DEST_PATH_IMAGE029
；将
Figure 647214DEST_PATH_IMAGE030
代入式（7），得到平板探测器所有点Dmn的坐标
Figure 110556DEST_PATH_IMAGE019
;
步骤F、计算X光机点光源S投影到平板探测器所有点Dmn的X线投影向量场
Figure 189370DEST_PATH_IMAGE031
，计算公式如下：
Figure 421769DEST_PATH_IMAGE032
其中
Figure 537623DEST_PATH_IMAGE033
；
在步骤B到步骤F计算生成X线图像投影向量场的全部过程中，参数yc和zc未知，参数
Figure 488262DEST_PATH_IMAGE034
为已知参数。
4.根据权利要求3所述的方法，其特征在于，在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；模拟X线投影人体生成X线图像的过程，将所述CT图像中的CT值转换成衰减系数，以所述系列投影向量场逐一投影CT图像，获得系列DRR图像，实现DRR 图像库的构建，具体包括：
步骤G、在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；具体包括：
步骤G1、指定坐标系原点；
令CT dicom数据中第一个被读取的体素点为坐标系原点；
步骤G2、采样投影向量场的未知参数；
在投影向量场
Figure 105188DEST_PATH_IMAGE031
未知参数yc和zc所在的yoz平面内，高频率、等间隔、全覆盖地采样两个参数，得到系列（yci, zcj）值；其中
Figure 192093DEST_PATH_IMAGE035
，A表示3D CT在yoz平面覆盖的区域，也即X线在yoz平面的照射区域；i = 0,1,…,I-1，j = 0,1,…,J-1，I和J分别表示在y轴和z轴上的采样点数；
步骤G3、获得系列投影向量场；
将采样值（yci ,zcj）逐一代入公式（8），获得系列投影向量场
Figure 993695DEST_PATH_IMAGE036
，其中i和j对应采样值（yci ,zcj），i = 0,1,…,I-1，j = 0,1,…,J-1；m和n对应平板探测器点Dmn，m = 0,1,…，RL-1，n = 0,1,…，RW-1；遍历(yci ,zcj)所有组合，也就遍历了投影空间内未知的所有投影变换；
步骤H、将CT图像中的CT值转化成X射线衰减系数：
为了模拟X线投影人体生成术中X线图像的过程，首先将所述3D CT图像中每一个体素的CT值转化为相应组织对X线的衰减系数，转化公式如下：
Figure 166051DEST_PATH_IMAGE037
其中，k表示X线穿过人体的第k类组织，
Figure 321089DEST_PATH_IMAGE038
表示第k类组织对X线的衰减系数，
Figure 793658DEST_PATH_IMAGE039
表示水对X线的衰减系数，HU表示CT值；
步骤I、以所述系列投影向量场逐一投影CT图像到平板探测器，获得系列DRR 图像，实现DRR图像库的构建；
模拟X线投影人体生成术中X线图像的过程，X线从点光源S出发，以所述系列投影向量场
Figure 375949DEST_PATH_IMAGE036
，i = 0,1,…,I-1，j = 0,1,…,J-1，逐一投影3D CT图像到平板探测器所有点Dmn，m= 0,1,…，RL-1，n = 0,1,…，RW-1，采用光线追踪法生成系列DRR 图像；点Dmn的灰度值DRR(m,n)计算如下：
Figure 910967DEST_PATH_IMAGE040
其中，I 0表示X线光源强度，指定值为1；k表示X线穿过人体的第k类组织，
Figure 869696DEST_PATH_IMAGE038
表示第k类组织的衰减系数；
Figure 931193DEST_PATH_IMAGE041
表示X线
Figure 949964DEST_PATH_IMAGE036
穿过第k类组织的长度，通过计算X线穿入和穿出第k类组织的两点间距得到；
投影向量场
Figure 221546DEST_PATH_IMAGE036
遍历m和n，得到对应采样参数（yci ,zcj）投影生成的DRR图像；以系列投影向量场
Figure 718386DEST_PATH_IMAGE036
，i = 0,1,…,I-1，j = 0,1,…,J-1，逐一投影3D CT图像，获得系列DRR 图像，实现DRR图像库的构建。
5.一种数字重建影像库的构建装置，其特征在于，包括：
数据获取单元，用于：获取一对术中X线图像与术前CT图像，所述术中X线图像由术中X光机采集，所述CT图像是对CT切片重建而成的3D CT图像；
X线投影向量场计算单元，用于：获取所述X光机投影人体生成所述术中X线图像时的投影定位部分信息，依据这些信息，计算得到所述X光机投影生成所述术中X线图像时的X线投影向量场，由此确定投影空间内的已知部分；
DRR图像库生成单元，用于：在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换；模拟X线投影人体生成X线图像的过程，将所述CT图像中的CT值转换成衰减系数，以所述系列投影向量场逐一投影CT图像，获得系列DRR 图像，实现DRR图像库的构建。
6.根据权利要求5所述的装置，其特征在于，所述X光机为C型臂X光机。
7.根据权利要求6所述的装置，其特征在于，所述X线投影向量场计算单元包括：获取投影信息子单元和计算X线投影向量场子单元，其中：
所述获取投影信息子单元，用于：从所述术中X光机说明书及术中X线图像的dicom数据中，读取X光机生成所述术中X线图像时的部分参数，包括：等中心点距地面高度、治疗床距地面高度、点光源沿X线投影方向到等中心点的距离、点光源沿X线投影方向到平板探测器中心点的距离、机架带动X光机围绕等中心点绕三个轴向旋转的旋转角、平板探测器尺寸、X线图像分辨率；
所述计算X线投影向量场子单元，用于：计算生成所述术中X线图像时的X线投影向量场，由此确定投影空间内的已知部分，具体包括：确定等中心点子单元、计算变换矩阵子单元、计算点光源坐标子单元、计算探测器坐标子单元，和计算投影向量场子单元，其中：
所述确定等中心点子单元，用于：确定等中心点C坐标(xc, yc, zc)中的x c ；
坐标系以垂直地面向上为x轴方向，以人体头部指向足端为z轴方向，采用右手法则确定y轴方向，并令X线照射区域人体最上端平行于治疗床的切面为平面yoz；等中心点C的坐标(xc, yc, zc)中，x c = H c- (H t+ H CT)，其中H c表示等中心点C距地面高度，H t为治疗床距地面高度，H CT表示CT切片的高度，也即X线照射区域人体厚度（S215）；yc和zc为未知参数；
所述计算变换矩阵子单元，用于：计算X光机以等中心点为中心、围绕三个轴向旋转的变换矩阵T；
首先平移坐标系，将原点移至等中心点C，计算平移矩阵T t ；然后计算X光机依次绕x、y、z三个轴向旋转α、β、γ角后所得的旋转矩阵R；最后再平移坐标系，将原点从等中心点C移回原坐标系位置，平移矩阵为T t 的逆矩阵T t -1 ；变换矩阵T的计算公式如下：
Figure 634389DEST_PATH_IMAGE001
Figure 824062DEST_PATH_IMAGE042
其中
Figure 723885DEST_PATH_IMAGE003
，
Figure 899783DEST_PATH_IMAGE004
即x轴，
Figure 670293DEST_PATH_IMAGE043
为绕
Figure 30867DEST_PATH_IMAGE004
旋转α角所对应的旋转矩阵；向量
Figure 152407DEST_PATH_IMAGE006
，是坐标系绕x轴旋转α角后y轴所对应向量，
Figure 115683DEST_PATH_IMAGE044
表示绕
Figure 6279DEST_PATH_IMAGE008
旋转β角所对应的旋转矩阵；向量
Figure 537754DEST_PATH_IMAGE009
，是坐标系先绕x轴旋转α角、再绕y轴旋转β角后z轴所对应向量，
Figure 881011DEST_PATH_IMAGE010
表示绕
Figure 673079DEST_PATH_IMAGE011
旋转γ角所对应的旋转矩阵；旋转矩阵
Figure 418182DEST_PATH_IMAGE012
依照如下公式计算：
Figure 120558DEST_PATH_IMAGE013
这里
Figure 216690DEST_PATH_IMAGE014
表示绕向量
Figure 397136DEST_PATH_IMAGE015
旋转θ角的旋转矩阵，将旋转矩阵
Figure 386958DEST_PATH_IMAGE012
的参数与
Figure 260236DEST_PATH_IMAGE014
的参数一一对应，即可计算得到各个矩阵值；
无论实际操作中X光机绕三个轴向旋转的先后次序如何，最终得到的旋转矩阵与按上述方法计算得到的旋转矩阵R相同；
所述计算点光源坐标子单元，用于：计算X光机点光源S的坐标(xs, ys, zs)，计算公式如下：
Figure 578085DEST_PATH_IMAGE016
其中，
Figure 562221DEST_PATH_IMAGE045
表示点光源S沿X线投影方向到等中心点C的距离，默认点光源S初始位置在治疗床下方，X线
Figure 157282DEST_PATH_IMAGE018
初始位置垂直于治疗床；
所述计算探测器坐标子单元，用于：计算平板探测器所有点Dmn的坐标
Figure 201461DEST_PATH_IMAGE019
，计算方式如下：
首先根据公式（6）计算平板探测器中心点Dc的坐标：
Figure 6606DEST_PATH_IMAGE020
其中
Figure 60013DEST_PATH_IMAGE021
表示点光源S沿X线投影方向到平板探测器中心点Dc的距离，然后根据公式（7）计算平板探测器所有点Dmn的坐标：
Figure 634213DEST_PATH_IMAGE022
其中，
Figure 973928DEST_PATH_IMAGE023
表示点Dmn的初始位置坐标，
Figure 790DEST_PATH_IMAGE024
，
Figure 592308DEST_PATH_IMAGE025
，SL*SW为平板探测器尺寸，RL*RW是X线图像分辨率；当
Figure 21015DEST_PATH_IMAGE026
，
Figure 282364DEST_PATH_IMAGE027
时，
Figure 62101DEST_PATH_IMAGE046
，由此推算出
Figure 191731DEST_PATH_IMAGE029
；将
Figure 474945DEST_PATH_IMAGE030
代入式（7），得到平板探测器所有点Dmn的坐标
Figure 562986DEST_PATH_IMAGE019
;
所述计算投影向量场子单元，用于：计算X光机点光源S投影到平板探测器所有点Dmn的X线投影向量场
Figure 689074DEST_PATH_IMAGE031
，计算公式如下：
Figure 356816DEST_PATH_IMAGE032
其中
Figure 760115DEST_PATH_IMAGE033
；
在上述计算X线投影向量场子单元的全部过程中，参数yc和zc为未知参数，参数
Figure 753479DEST_PATH_IMAGE034
为已知参数。
8.根据权利要求7所述的装置，其特征在于，所述DRR图像库生成单元包括：系列投影向量场计算子单元，CT图像预处理子单元，和生成DRR图像库子单元，其中：
所述系列投影向量场计算子单元，用于：在所述X线投影向量场未知的参数空间内采样未知投影参数，组合所述X线投影向量场与所述采样的未知投影参数，得到系列投影向量场，由此遍历投影空间内未知的投影变换，具体又包括：
确定坐标系原点子单元，用于：令CT dicom数据中第一个被读取的体素点为坐标系原点；
采样未知参数子单元，用于：在投影向量场
Figure 852016DEST_PATH_IMAGE031
未知参数yc和zc所在的yoz平面内，高频率、等间隔、全覆盖地采样两个参数，得到系列（yci, zcj）值；其中
Figure 589028DEST_PATH_IMAGE035
，A表示3D CT在yoz平面覆盖的区域，即X线在yoz平面的照射区域；i = 0,1,…,I-1，j = 0,1,…,J-1，I和J分别表示在y轴和z轴上的采样点数；
计算系列投影向量场子单元，用于：将采样值（yci ,zcj）逐一代入公式（8），获得系列投影向量场
Figure 846834DEST_PATH_IMAGE036
，其中i和j对应采样值（yci ,zcj），i = 0,1,…,I-1，j = 0,1,…,J-1；m和n对应平板探测器点Dmn，m = 0,1,…，RL-1，n = 0,1,…，RW-1；遍历(yci ,zcj)所有组合，也就遍历了投影空间内未知的所有投影变换；
所述CT图像预处理子单元，用于：为了模拟X线投影人体生成术中X线图像的过程，首先将所述3D CT图像中每一个体素的CT值转化为相应组织对X线的衰减系数，转化公式如下：
Figure 745520DEST_PATH_IMAGE037
其中，k表示X线穿过人体的第k类组织，
Figure 846200DEST_PATH_IMAGE038
表示第k类组织对X线的衰减系数，
Figure 121324DEST_PATH_IMAGE039
表示水对X线的衰减系数，HU表示CT值；
所述生成DRR图像库子单元，用于：模拟X线投影人体生成术中X线图像的过程，X线从点光源S出发，以所述系列投影向量场
Figure 233636DEST_PATH_IMAGE036
，i = 0,1,…,I-1，j = 0,1,…,J-1，逐一投影3DCT图像到平板探测器所有点Dmn，m = 0,1,…，RL-1，n = 0,1,…，RW-1，采用光线追踪法生成系列DRR 图像；点Dmn的灰度值DRR(m,n)计算如下：
Figure 568802DEST_PATH_IMAGE040
其中，I 0表示X线光源强度，指定值为1；k表示X线穿过人体的第k类组织，
Figure 639002DEST_PATH_IMAGE038
表示第k类组织的衰减系数；
Figure 717817DEST_PATH_IMAGE041
表示X线
Figure 950215DEST_PATH_IMAGE036
穿过第k类组织的长度，通过计算X线穿入和穿出第k类组织的两点间距得到；
投影向量场
Figure 190703DEST_PATH_IMAGE036
遍历m和n，得到对应采样参数（yci ,zcj）投影生成的DRR图像；以系列投影向量场
Figure 406921DEST_PATH_IMAGE036
，i = 0,1,…,I-1，j = 0,1,…,J-1，逐一投影3D CT图像，获得系列DRR 图像，实现DRR图像库的构建。
CN202110548584.6A 2021-05-20 2021-05-20 一种数字重建影像库的构建方法和装置 Expired - Fee Related CN112989081B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110548584.6A CN112989081B (zh) 2021-05-20 2021-05-20 一种数字重建影像库的构建方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110548584.6A CN112989081B (zh) 2021-05-20 2021-05-20 一种数字重建影像库的构建方法和装置

Publications (2)

Publication Number Publication Date
CN112989081A CN112989081A (zh) 2021-06-18
CN112989081B true CN112989081B (zh) 2021-08-27

Family

ID=76337030

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110548584.6A Expired - Fee Related CN112989081B (zh) 2021-05-20 2021-05-20 一种数字重建影像库的构建方法和装置

Country Status (1)

Country Link
CN (1) CN112989081B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115205417B (zh) * 2022-09-14 2023-03-14 首都医科大学附属北京安贞医院 一种投影变换的计算方法、装置、设备及存储介质
CN116433476B (zh) * 2023-06-09 2023-09-08 有方（合肥）医疗科技有限公司 Ct图像处理方法及装置

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7366278B2 (en) * 2004-06-30 2008-04-29 Accuray, Inc. DRR generation using a non-linear attenuation model
CN106875376B (zh) * 2016-12-29 2019-10-22 中国科学院自动化研究所 腰椎配准先验模型的构建方法以及腰椎配准方法
US10434335B2 (en) * 2017-03-30 2019-10-08 Shimadzu Corporation Positioning apparatus and method of positioning by generation of DRR image from X-ray CT image data
WO2019012686A1 (ja) * 2017-07-14 2019-01-17 株式会社日立製作所 粒子線治療装置およびｄｒｒ画像作成方法
CN112614169B (zh) * 2020-12-24 2022-03-25 电子科技大学 基于深度学习网络的2d/3d脊椎ct层级配准方法

Also Published As

Publication number Publication date
CN112989081A (zh) 2021-06-18

Similar Documents

Publication Publication Date Title
US11707241B2 (en) System and method for local three dimensional volume reconstruction using a standard fluoroscope
JP5345947B2 (ja) 対象物を撮像する撮像システム及び撮像方法
CN105916444B (zh) 由二维x射线图像来重建三维图像的方法
Penney et al. A comparison of similarity measures for use in 2-D-3-D medical image registration
US11559266B2 (en) System and method for local three dimensional volume reconstruction using a standard fluoroscope
CN112989081B (zh) 一种数字重建影像库的构建方法和装置
US10433797B2 (en) Systems and methods for ultra low dose CT fluoroscopy
JP6637781B2 (ja) 放射線撮像装置及び画像処理プログラム
JP6349278B2 (ja) 放射線撮像装置、画像処理方法及びプログラム
US11127153B2 (en) Radiation imaging device, image processing method, and image processing program
CN109350059B (zh) 用于肘部自动对准的组合的转向引擎和界标引擎
US8358874B2 (en) Method for displaying computed-tomography scans, and a computed-tomography system or computed-tomography system assembly for carrying out this method
TWI836493B (zh) 註冊二維影像資料組與感興趣部位的三維影像資料組的方法及導航系統
CN104700377B (zh) 获得对计算机断层扫描数据进行射束硬化校正的射束硬化校正系数的方法和装置
WO2006085253A2 (en) Computer tomography apparatus, method of examining an object of interest with a computer tomography apparatus, computer-readable medium and program element
JP6703470B2 (ja) データ処理装置及びデータ処理方法
US7116808B2 (en) Method for producing an image sequence from volume datasets
US10682184B2 (en) Tissue sampling system
Zheng et al. A unifying MAP-MRF framework for deriving new point similarity measures for intensity-based 2D-3D registration
JP6505462B2 (ja) 医用画像処理装置、ｘ線コンピュータ断層撮影装置
Shu et al. Isocentric fixed angle irradiation-based DRR: a novel approach to enhance x-ray and CT image registration
CN117726672A (zh) 一种介入设备实时三维成像方法、***和装置
CN116704068A (zh) 一种用于双源ct设备的图像校正方法、成像方法和***

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for 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: 20210827

CF01 Termination of patent right due to non-payment of annual fee