CN115205417B - 一种投影变换的计算方法、装置、设备及存储介质 - Google Patents

一种投影变换的计算方法、装置、设备及存储介质 Download PDF

Info

Publication number
CN115205417B
CN115205417B CN202211113427.3A CN202211113427A CN115205417B CN 115205417 B CN115205417 B CN 115205417B CN 202211113427 A CN202211113427 A CN 202211113427A CN 115205417 B CN115205417 B CN 115205417B
Authority
CN
China
Prior art keywords
image
coordinate system
ray
point
projection
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.)
Active
Application number
CN202211113427.3A
Other languages
English (en)
Other versions
CN115205417A (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 CN202211113427.3A priority Critical patent/CN115205417B/zh
Publication of CN115205417A publication Critical patent/CN115205417A/zh
Application granted granted Critical
Publication of CN115205417B publication Critical patent/CN115205417B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • A61B6/5229—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data combining image data of a patient, e.g. combining a functional image with an anatomical image
    • A61B6/5235—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data combining image data of a patient, e.g. combining a functional image with an anatomical image combining images from the same or different ionising radiation imaging techniques, e.g. PET and CT
    • 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/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biophysics (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明提供了一种投影变换的计算方法、装置、设备及存储介质。采用多个标记点对患者不同部位进行体表标记，获取患者的一对术前CT图像与术中X线图像；依据采集CT图像时的扫描参数，计算得到CT图像上所有标记点3D体素坐标的第一3D空间位置；依据术中X光机投影患者生成X线图像时的投影已知参数，计算得到X线图像上所有标记点2D像素坐标的第二3D空间位置；依据所有标记点的第一3D空间位置和第二3D空间位置，计算术中X光机投影患者生成X线图像时的投影未知参数，得到配准术前CT图像与术中X线图像的投影变换。本发明有效改善了现有投影变换计算方法的不足，显著提升了投影变换计算的准确性。

Description

一种投影变换的计算方法、装置、设备及存储介质
技术领域
本发明涉及CT和X线图像配准领域，具体涉及一种投影变换的计算方法、装置、设备及存储介质。
背景技术
配准CT和X线图像的正确投影变换，是训练和评价有监督的CT和X线图像深度学习配准网络模型的重要依据。因此，通过两图中标注的同名点对计算得到配准CT和X线图像的投影变换，其准确与否，对于训练和评价有监督的CT和X线图像配准网络模型至关重要，直接关系到配准网络模型的性能优劣。
配准CT和X线图像的过程，我们可以视为将术前采集重建的3D CT与术中患者身体重叠的过程。如果两者对应位置完全重叠，那么以术前3D CT替代术中人体，模拟X线照射人体，以生成术中X线图像的投影变换对3D CT进行投影，在探测器上最后生成的投影图像，就是与术中X线图像配准的3D CT的投影图像；反之，与术中X线图像配准的3D CT的投影图像，其所对应的投影变换，就是配准CT和X线图像的投影变换。
现有投影变换计算方法所采用的模拟投影方式有两种：1）以3D CT中心为投影变换空间的原点；2）以点光源为投影变换空间的原点。两种方式皆以点光源穿过3D CT中心的X线作为X线光束中心线，并以此为基础确定光源到等中心的距离（source-axis distance,SAD）和光源到探测器中心的距离（source-isocenter distance, SID）；两种方式均以原点为旋转中心，绕X、Y、Z三个轴向旋转；三个旋转角度未知，SID和SAD为估计值。
事实上，这两种模拟投影方式均不符合实际。真实情况下，医生在C型臂下实施手术，术中C型臂X光机采用等中心定角照射的方式投影患者生成X线图像。与遵循等中心定角照射的真实X线投影相比较，现有两种模拟的投影方式存在如下不同：
1）真实X线投影以C型臂等中心点为旋转中心，而现有模拟X线投影以3D CT中心或者点光源为旋转中心；
2）真实X线投影以L臂、枢轴、C臂旋转轴为三个旋转轴，L臂独立于枢轴与C臂旋转轴，三者并未始终保持相互垂直的关系，而现有模拟X线投影以坐标系X、Y、Z轴为旋转轴，三者始终相互垂直；
3）真实X线投影以从点光源出发穿过等中心点的X线作为X线光束中心线，而现有模拟X线投影皆以从点光源出发穿过3D CT中心的X线作为X线光束中心线；等中心点不同于3D CT中心，因此两者对SID和SAD的设定不一致；
4）真实X线投影的投影变量为：L臂旋转角α，枢轴旋转角β，C臂旋转角γ，点光源S到等中心点C的距离dSAD、点光源S 到平板探测器中心点DC的距离dSID，以及等中心点坐标(xc,yc,zc)共8个参数；而现有模拟X线投影的投影变量为：绕坐标系X、Y、Z三个轴旋转的三个旋转角，以及沿X、Y、Z三个轴向的平移，共6个参数；
5）真实X线投影提供有已知生成X线图像时的L臂旋转角α，枢轴旋转角β，C臂旋转角γ，点光源S到等中心点C的距离dSAD、点光源S 到平板探测器中心点Dc的距离dSID共五个参数的值；而现有模拟X线投影没有生成X线图像时的准确投影参数值。
因为与真实投影存在上述明显不同，所以采用现有模拟投影方式，通过两图中标注的同名点对计算得到的投影变换不可能是准确的投影变换，以此训练和评价深度学习的CT与X线配准网络，所得到的配准网络模型性能堪忧。
发明内容
鉴于上述原因，本发明提出了一种基于等中心定角照射的投影变换的计算方法、装置、设备及存储介质，有效克服了现有投影变换计算方法存在的不足，使得计算得到投影变换更加逼近真实的X线投影变换，所获得的投影参数更为准确。
第一方面，本发明提供了一种投影变换的计算方法，所述方法包括：
采用多个标记点对患者不同部位进行体表标记，获取患者的一对术前CT与术中X线图像，所述CT图像是由CT切片重建而成的3D图像，所述X线图像则由术中X光机采集；
获取采集所述CT图像时的扫描参数，以及所述CT图像上所有标记点的3D体素坐标，计算得到所述CT图像上所有标记点3D体素坐标的第一3D空间位置；
获取所述术中X光机投影患者生成所述X线图像时的投影已知参数，以及所述X线图像所有标记点的2D像素坐标，计算得到所述X线图像上所有标记点2D像素坐标的第二3D空间位置；
依据所述所有标记点的第一3D空间位置和第二3D空间位置，计算所述术中X光机投影患者生成所述X线图像时的投影未知参数，得到配准术前CT图像与术中X线图像的投影变换。
可选地，所述术中X光机为C型臂X光机。
第二方面，基于上述投影变换的计算方法，本发明还提供了一种投影变换的计算装置，包括：
图像获取模块，用于采用多个标记点对患者不同部位进行体表标记，获取患者的一对术前CT与术中X线图像，所述CT图像是由CT切片重建而成的3D图像，所述X线图像则由术中X光机采集；
第一计算模块，用于获取采集所述CT图像时的扫描参数，以及所述CT图像上所有标记点的3D体素坐标，计算得到所述CT图像上所有标记点3D体素坐标的第一3D空间位置；
第二计算模块，用于获取所述术中X光机投影患者生成所述X线图像时的投影已知参数，以及所述X线图像所有标记点的2D像素坐标，计算得到所述X线图像上所有标记点2D像素坐标的第二3D空间位置；
第三计算模块，用于依据所述所有标记点的第一3D空间位置和第二3D空间位置，计算所述术中X光机投影患者生成所述X线图像时的投影未知参数，得到配准术前CT图像与术中X线图像的投影变换。
第三方面，基于上述投影变换的计算方法，本发明还提供了一种计算机设备，包括：存储器和处理器，所述存储器和所述处理器之间互相通信连接，所述存储器中存储有计算机指令，所述处理器通过执行所述计算机指令，从而执行第一方面，或者第一方面任意一种可选实施方式中所述的方法。
第四方面，基于上述投影变换的计算方法，本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于使所述计算机执行第一方面，或者第一方面任意一种可选实施方式中所述的方法。
相对于现有的传统投影参数计算技术，本发明可以显著提升投影变换计算的准确性，原因在于：
1）本发明完全遵循术中X光机投影患者生成X线图像时所采用的等中心定角照射原则，投影以C型臂等中心点为旋转中心，以L臂及与之独立的枢轴和C臂旋转轴为三个旋转轴，而现有投影变换计算所采用的投影则以3D CT中心或者点光源为旋转中心，以坐标系X、Y、Z轴为旋转轴，投影的运动形式与实际不符；
2）本发明采用与真实X线投影完全相同的投影变量，包括L臂旋转角α、枢轴旋转角β，C臂旋转角γ，点光源S到等中心点C的距离dSAD、点光源S 到平板探测器中心点DC的距离dSID，以及等中心点的空间位置(xc,yc,zc)共8个参数；而现有投影变换计算所采用的投影变量则为：绕坐标系X、Y、Z三个轴旋转的三个旋转角，以及沿X、Y、Z三个轴向的平移，共6个参数，参数的设置和数量均与实际不符；
3）本发明在计算投影变换时，先从X线图像数据中读取L臂旋转角α，枢轴旋转角β，C臂旋转角γ，点光源S到等中心点C的距离dSAD、点光源S 到平板探测器中心点Dc的距离dSID共五个参数的值用作投影变换的已知参数，投影变换仅需计算等中心点坐标(xc,yc,zc)三个未知投影参数；而现有投影变换计算并没有利用X线图像所提供的已知投影参数，需要计算投影变换全部6个参数的值。
简言之，本发明将现有投影变换计算中不符合实际的投影运动形式、不符合实际的投影参数设置和数量全部纠正为符合实际情况；同时充分利用X线图像所提供的已知投影参数值，大幅减少了需要计算的投影参数数量；因而能有效改善现有投影变换计算的不足，显著提升投影变换计算的准确性。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例1提供的一种投影变换计算方法的计算流程示意图；
图2为本发明实施例1提供的CT图像中标记点第一3D空间位置的计算流程示意图；
图3为本发明实施例1提供的X线图像中标记点第二3D空间位置的计算流程示意图；
图4为本发明实施例1提供的C型臂及其投影方式和投影参数示意图；
图5为本发明实施例1提供的3D等中心坐标系中术中X光机相关各点的初始3D空间位置的计算流程示意图；
图6为本发明实施例1提供的3D患者坐标系中X线图像中标记点第二3D空间位置的计算流程示意图；
图7为本发明实施例1提供的依据标记点第一和第二3D空间位置计算投影未知参数的流程示意图；
图8为本发明实施例1提供的3D患者坐标系下CT与X线图像中标记点对的对应关系示意图；
图9为本发明实施例2提供的一种投影变换的计算装置的结构示意图；
图10为本发明实施例3提供的计算机设备的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1
本发明实施例提供了一种投影变换的计算方法，如图1所示，包括以下步骤：
步骤S1：采用多个标记点对患者不同部位进行体表标记，获取带有体表标记患者的一对术前CT与术中X线图像，CT图像是由CT切片重建而成的3D图像，X线图像则由术中X光机采集。示例性地，术中X光机为C型臂X光机，但不以此为限。
步骤S2：获取采集CT图像时的扫描参数，以及CT图像上所有标记点的3D体素坐标，计算得到CT图像上所有标记点3D体素坐标的第一3D空间位置；如图2所示，具体步骤包括：
步骤S21：从术中CT图像的dicom数据中，读取采集CT图像时的扫描参数，扫描参数包括：CT切片的分辨率RCL*RCW、扫描层间距dd、扫描孔径da；同时读取CT图像上所有标记点的3D体素坐标Qi3(xi3,yi3,zi3)，i=1,…,N；
步骤S22：确定3D患者坐标系：
3D患者坐标系中，指定以3D CT图像对应患者头端左侧背部第一个体素为坐标原点，以垂直地面向上为x轴正方向，以患者头部指向足端为z轴正方向，采用右手法则确定y轴正方向；
步骤S23：计算3D患者坐标系中CT图像标记点第一3D空间位置：
在3D患者坐标系下，计算得到CT图像上所有标记点的第一3D空间位置为：
Figure 979853DEST_PATH_IMAGE001
Figure 10000230688111
i=1,…,N。
步骤S3：获取术中X光机投影患者生成X线图像时的投影已知参数，以及X线图像所有标记点的2D像素坐标，计算得到X线图像上所有标记点2D像素坐标的第二3D空间位置。如图3所示，具体步骤包括：
步骤S31：获取术中X光机投影患者生成X线图像时的投影已知参数：
从术中X线图像的dicom数据及术中X光机说明书中，读取投影已知参数，包括： L臂旋转角α，枢轴旋转角β，C臂旋转角γ，点光源S到等中心点C的距离dSAD、点光源S 到平板探测器中心点DC的距离dSID、平板探测器尺寸SL*SW、X线图像分辨率RL*RW，具体参数如图4所示；
大型C型臂X光机与中小型C型臂X光机所提供的投影参数有所不同，但上述参数都能从X线图像dicom数据标签及C型臂说明书中读取。
步骤S32：确定等中心点、点光源及X线图像标记点的初始3D空间位置，包括如下步骤：
步骤S321：确定3D等中心坐标系和X线图像2D坐标系：
3D等中心坐标系中，指定以等中心点C为原点，以垂直地面向上为x轴正方向，以患者头部指向足端为z轴正方向，采用右手法则确定y轴正方向；
X线图像2D坐标系中，指定以X线图像左上角为原点，以患者头部指向足端为x轴正方向，以X线图像从左到右为y轴正方向；
步骤S322：计算3D等中心坐标系中术中X光机相关各点的初始3D空间位置，如图5所示，具体步骤包括：
步骤S3221：预设C型臂X光机的初始位置为：点光源S在治疗床下方，探测器平面平行于水平面，点光源到探测器平面中心点的连线SDc始终垂直于探测器平面；
步骤S3222：获取X线图像上所有标记点的2D像素坐标Di20(xi2,yi2 )，i=1,…,N；
步骤S3223：在3D等中心坐标系下，根据等中心C初始3D空间位置C30’(0,0,0)和点光源S到等中心点C的距离dSAD，计算得到点光源S的初始3D空间位置为S30’(-dSAD,0,0)；
步骤S3234：在3D等中心坐标系下，根据点光源S 到平板探测器中心点Dc的距离dSID，计算得到探测器中心点Dc的初始3D空间位置为Dc30’(dSID-dSAD,0,0)，以及X线图像上所有标记点的初始3D空间位置为：
Figure 236150DEST_PATH_IMAGE004
。
步骤S33：计算3D患者坐标系中X线图像上所有标记点的第二3D空间位置，如图6所示，具体步骤包括：
步骤S331：在3D等中心坐标系下计算旋转矩阵：
3D等中心坐标系下，从C型臂初始位置开始，机架带动术中X光机和探测器围绕等中心点C旋转，依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角，由此计算得到旋转矩阵为：
Figure 936122DEST_PATH_IMAGE005
其中L臂、枢轴、C臂旋转轴的初始位置分别为
Figure 272425DEST_PATH_IMAGE006
，
Figure DEST_PATH_IMAGE007
为术中X光机绕L臂旋转α角所对应的旋转矩阵；绕L臂旋转α角后，枢轴从初始位置
Figure 856991DEST_PATH_IMAGE008
旋转到
Figure 837716DEST_PATH_IMAGE009
，
Figure 583955DEST_PATH_IMAGE010
为术中X光机绕枢轴旋转β角所对应的旋转矩阵；C臂旋转轴先绕L臂旋转α角，再绕枢轴旋转β角后，C臂旋转轴从初始位置
Figure 141975DEST_PATH_IMAGE011
旋转到
Figure 264652DEST_PATH_IMAGE012
，
Figure 847687DEST_PATH_IMAGE013
为术中X光机绕C臂旋转轴
Figure 764827DEST_PATH_IMAGE014
旋转γ角所对应的旋转矩阵；旋转矩阵
Figure 75723DEST_PATH_IMAGE015
依照如下公式计算：
Figure 736511DEST_PATH_IMAGE016
这里
Figure 426250DEST_PATH_IMAGE017
表示绕向量
Figure 779871DEST_PATH_IMAGE018
旋转θ角的旋转矩阵，将旋转矩阵
Figure 312483DEST_PATH_IMAGE019
的参数与
Figure 776963DEST_PATH_IMAGE020
的参数一一对应，计算各个矩阵值，得到旋转矩阵，需要说明的是，无论术中X光机绕机臂旋转的先后次序如何，得到的旋转矩阵结果均一致；
步骤S332： 在3D患者坐标系下计算术中X光机相关各点的第二3D空间位置，包括如下步骤：
步骤S3321：在3D患者坐标系中，预设等中心点C的第二3D空间位置为C3(xc,yc,zc)；
步骤S3322：在3D患者坐标系中，经旋转变换后，逐一计算得到：
点光源S的第二3D空间位置变换为S3=RS30’+C3，
探测器中心点Dc的第二3D空间位置变换为Dc3=RDc30’+C3，
X线图像上所有标记点的第二3D空间位置为Di3=RDi30’+C3，i=1,…，N，
其中，R为旋转矩阵，S30’、Dc30’、Di30’分别表示点光源S、探测器中心点DC和X线图像第i个标记点在3D等中心坐标系下的初始3D空间位置。
需要说明的是，术中X光机包括C型臂X光机，G型臂X光机和O型臂X光机三个机型。其中，C型臂X光机是最常见的术中X光机，采用单术中X光机***；G型臂X光机和O型臂X光机是C型臂X光机的改进机型，均采用双术中X光机***。就G型臂X光机或O型臂X光机中的单组术中X光机而言，其运动模式和投影方式与C型臂X光机大同小异，因此上述计算过程同样适用。
步骤S4：依据所有标记点的第一3D空间位置和第二3D空间位置，计算术中X光机投影患者生成X线图像时的投影未知参数，得到配准术前CT图像与术中X线图像的投影变换；如图7所示，具体步骤包括：
步骤S41：根据图8所示3D患者坐标系下CT与X线图像中标记点对的空间对应关系，计算3D患者坐标系中点光源到CT图像上第i个标记点Qi3的3D空间向量：
Figure 570475DEST_PATH_IMAGE021
同时，计算3D患者坐标系中点光源到X线图像上第i个标记点Di3的3D空间向量：
Figure 94998DEST_PATH_IMAGE022
其中，R为旋转矩阵，i=1,…，N，S30’和Di30’分别表示点光源S和X线图像上第i个标记点在3D等中心坐标系下的初始3D空间位置， S3、C3、Di3分别表示点光源S、等中心点C、X线图像上第i个标记点在3D患者坐标系下的第二3D空间位置，Qi3则表示CT图像上第i个标记点在3D患者坐标系下的第一3D空间位置，i=1,…,N；
步骤S42：计算得到一组投影未知参数xc,yc,zc的值：
CT图像上第i个标记点Qi3和X线图像上第i个标记点Di3为同名点对，根据向量
Figure DEST_PATH_IMAGE023
与
Figure 521431DEST_PATH_IMAGE024
共线，得到：
Figure 524022DEST_PATH_IMAGE025
其中，
Figure 312986DEST_PATH_IMAGE026
R为旋转矩阵，i=1,…，N，S30’和Di30’分别表示点光源S和X线图像上第i个标记点在3D等中心坐标系下的初始3D空间位置，Qi3表示CT图像上第i个标记点在3D患者坐标系下的第一3D空间位置；
步骤S43：计算得到最终投影未知参数xc,yc,zc的值，包括：
任意组合两个标记点在CT与X线图像中的第一和第二3D 空间位置(Qi3，Di3)和(Qj3，Dj3)，i=1,…,N，j=1,…,N，根据所述上述步骤S41和S42计算得到一组xc,yc,zc的值；组合多个两两标记点，得到多组xc,yc,zc的值；取所有组xc,yc,zc的平均值，作为最终投影未知参数xc,yc,zc的值。
相对于现有投影变换计算技术，一方面，针对其中不符合实际的投影运动形式，以及不符合实际的投影参数设置和数量，本发明将其全部转换为符合实际情况；另一方面，针对现有技术忽视X线图像所提供的已知投影参数，本发明充分利用已知投影参数值，大幅减少需要计算的投影参数量；从两个方面着手，本发明有效改善现有投影变换计算的不足，显著提升了投影变换计算的准确性。
实施例2
本发明实施例提供一种投影变换的计算装置，如图9所示，包括：
图像获取模块1，用于采用多个标记点对患者不同部位进行体表标记，获取患者的一对术前CT与术中X线图像，CT图像是由CT切片重建而成的3D图像，X线图像则由术中X光机采集。详细内容参见上述实施例1中步骤S1的描述，在此不再赘述。
第一计算模块2，用于获取采集CT图像时的扫描参数，以及CT图像上所有标记点的3D体素坐标，计算得到CT图像上所有标记点3D体素坐标的第一3D空间位置。详细内容参见上述实施例1中步骤S2的描述，在此不再赘述。
第二计算模块3，用于获取术中X光机投影患者生成X线图像时的投影已知参数，以及X线图像所有标记点的2D像素坐标，计算得到X线图像上所有标记点2D像素坐标的第二3D空间位置。详细内容参见上述实施例1中步骤S3的描述，在此不再赘述。
第三计算模块4，用于依据所有标记点的第一3D空间位置和第二3D空间位置，计算术中X光机投影患者生成X线图像时的投影未知参数，得到配准术前CT图像与术中X线图像的投影变换。详细内容参见上述实施例1中步骤S4的描述，在此不再赘述。
相对于现有投影变换计算技术，一方面，针对其中不符合实际的投影运动形式，以及不符合实际的投影参数设置和数量，本发明将其全部转换为符合实际情况；另一方面，针对现有技术忽视X线图像所提供的已知投影参数，本发明充分利用已知投影参数值，大幅减少需要计算的投影参数量；从两个方面着手，本发明有效改善现有投影变换计算的不足，显著提升了投影变换计算的准确性。
实施例3
图10示出了本发明实施例中计算机设备的结构示意图，如图10所示，包括：处理器1001和存储器1002，其中，处理器1001和存储器1002可以通过总线或者其他方式连接，图10中以通过总线连接为例。
处理器1001可以为中央处理器（Central Processing Unit，CPU）。处理器1001还可以为其他通用处理器、数字信号处理器（Digital Signal Processor，DSP）、专用集成电路（Application Specific Integrated Circuit，ASIC）、现场可编程门阵列（Field-Programmable Gate Array，FPGA）或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等芯片，或者上述各类芯片的组合。
存储器1002作为一种非暂态计算机可读存储介质，可用于存储非暂态软件程序、非暂态计算机可执行程序以及模块，如上述方法实施例中的方法所对应的程序指令/模块。处理器1001通过运行存储在存储器1002中的非暂态软件程序、指令以及模块，从而执行处理器的各种功能应用以及数据处理，即实现上述实施例1中的方法。
存储器1002可以包括存储程序区和存储数据区，其中，存储程序区可存储操作***、至少一个功能所需要的应用程序；存储数据区可存储处理器1001所创建的数据等。此外，存储器1002可以包括高速随机存取存储器，还可以包括非暂态存储器，例如至少一个磁盘存储器件、闪存器件、或其他非暂态固态存储器件。在一些实施例中，存储器1002可选包括相对于处理器1001远程设置的存储器，这些远程存储器可以通过网络连接至处理器1001。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。
一个或者多个模块存储在存储器1002中，当被处理器1001执行时，执行上述实施例1中的方法。
上述计算机设备具体细节可以对应参阅上述方法实施例中对应的相关描述和效果进行理解，此处不再赘述。
本领域技术人员可以理解，实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，实现的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述实施例1的全部流程。其中，存储介质可为磁碟、光盘、只读存储记忆体（Read-Only Memory，ROM）、随机存储记忆体（Random Access Memory，RAM）、快闪存储器（Flash Memory）、硬盘（Hard Disk Drive，缩写：HDD）或固态硬盘（Solid-StateDrive，SSD)等；存储介质还可以包括上述种类的存储器的组合。
需要说明的是，在本发明的描述中，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，不应理解为对本发明的限制；术语“第一”、“第二”仅用于描述目的，不应理解为指示或暗示相对重要性；而术语“安装”、“相连”、“连接”也应做广义理解。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
上述的不同实施方式中所涉及的技术特征，只要彼此之间未构成冲突，就可以相互结合。本领域技术人员在不脱离本发明的精神和范围的情况下作出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。

Claims (7)

1.一种投影变换的计算方法，其特征在于，所述方法包括：
采用多个标记点对患者不同部位进行体表标记，获取患者的一对术前CT图像与术中X线图像，所述CT图像是由CT切片重建而成的3D图像，所述X线图像则由术中X光机采集；
获取采集所述CT图像时的扫描参数，以及所述CT图像上所有标记点的3D体素坐标，计算得到所述CT图像上所有标记点3D体素坐标的第一3D空间位置；
获取所述术中X光机投影患者生成所述X线图像时的投影已知参数，以及所述X线图像所有标记点的2D像素坐标，计算得到所述X线图像上所有标记点2D像素坐标的第二3D空间位置，具体包括：
步骤31：获取所述术中X光机投影患者生成所述X线图像时的投影已知参数：
从所述术中X线图像的dicom数据及所述术中X光机说明书中，读取投影已知参数，所述投影已知参数包括： L臂旋转角α，枢轴旋转角β，C臂旋转角γ，点光源S到等中心点C的距离dSAD、点光源S 到平板探测器中心点DC的距离dSID、平板探测器尺寸SL*SW、X线图像分辨率RL*RW；
步骤32：确定等中心点、点光源及X线图像标记点的初始3D空间位置，包括：
步骤321：确定3D等中心坐标系和X线图像2D坐标系；
步骤322：计算3D等中心坐标系中X光机相关各点的初始3D空间位置；
步骤33：计算3D患者坐标系中X线图像上所有标记点的第二3D空间位置：
步骤331：在3D等中心坐标系下计算旋转矩阵；
步骤332：在3D患者坐标系下计算X光机相关各点的第二3D空间位置；
依据所述所有标记点的第一3D空间位置和第二3D空间位置，计算所述术中X光机投影患者生成所述X线图像时的投影未知参数，得到配准术前CT图像与术中X线图像的投影变换，具体包括：
步骤41：计算3D患者坐标系中点光源S到CT图像上第i个标记点Qi3的3D空间向量：
Figure QLYQS_1
同时计算3D患者坐标系中点光源S到X线图像上第i个标记点Di3的3D空间向量：
Figure QLYQS_2
其中，R为依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角的旋转矩阵，i=1,…，N，S30’和Di30’分别表示所述点光源S和X线图像上第i个标记点在3D等中心坐标系下的初始3D空间位置，S3、C3、Di3分别表示点光源S、等中心点C、X线图像上第i个标记点在3D患者坐标系下的第二3D空间位置，Qi3则表示CT图像上第i个标记点在3D患者坐标系下的第一3D空间位置，i=1,…,N；
步骤42：计算得到一组投影未知参数xc,yc,zc的值，包括：
CT图像上第i个标记点Qi3和X线图像上第i个标记点Di3为同名点对，根据向量
Figure QLYQS_3
与
Figure QLYQS_4
共线，得到：
Figure QLYQS_5
其中，
Figure QLYQS_6
R为依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角的旋转矩阵，i=1,…,N，S30’和Di30’分别表示所述点光源S和X线图像上第i个标记点在3D等中心坐标系下的初始3D空间位置，Qi3表示CT图像上第i个标记点在3D患者坐标系下的第一3D空间位置；
由此计算得到一组投影未知参数xc,yc,zc的值；
步骤43：计算得到最终投影未知参数xc,yc,zc的值，包括：
任意组合两个标记点在CT与X线图像中的第一和第二3D 空间位置(Qi3，Di3)和(Qj3，Dj3)，i=1,…,N，j=1,…,N，根据所述步骤41和步骤42计算得到一组xc,yc,zc的值；组合多个两两标记点，得到多组xc,yc,zc的值；取所有组xc,yc,zc的平均值，作为最终投影未知参数xc,yc,zc的值。
2.根据权利要求1所述的投影变换的计算方法，其特征在于，获取采集所述CT图像时的扫描参数，以及所述CT图像上所有标记点的3D体素坐标，计算得到所述CT图像上所有标记点3D体素坐标的第一3D空间位置，具体包括：
步骤21：从所述术前CT图像的dicom数据中，读取采集所述CT图像时的扫描参数，所述扫描参数包括：CT切片的分辨率RCL*RCW、扫描层间距dd、扫描孔径da；同时读取CT图像上所有标记点的3D体素坐标Qi3(xi3,yi3,zi3)，i=1,…,N；
步骤22：确定3D患者坐标系，包括：
3D患者坐标系中，指定以3D CT图像对应患者头端左侧背部第一个体素为坐标原点，以垂直地面向上为x轴方向，以患者头部指向足端为z轴方向，采用右手法则确定y轴方向；
步骤23：计算3D患者坐标系中CT图像标记点第一3D空间位置，包括：
在3D患者坐标系下，计算得到CT图像上所有标记点的第一3D空间位置为：
Figure QLYQS_7
，其中，i=1,…,N。
3.根据权利要求2所述投影变换的计算方法，其特征在于，所述术中X光机为C型臂X光机。
4.根据权利要求3所述投影变换的计算方法，其特征在于，
步骤321：确定3D等中心坐标系和X线图像2D坐标系：
3D等中心坐标系中，指定以等中心点C为原点，以垂直地面向上为x轴正方向，以患者头部指向足端为z轴正方向，采用右手法则确定y轴正方向；
X线图像2D坐标系中，指定以X线图像左上角为原点，以患者头部指向足端为x轴正方向, 以X线图像从左到右为y轴正方向；
步骤322：计算3D等中心坐标系中术中X光机相关各点的初始3D空间位置，包括：
步骤3221：预设C型臂X光机的初始位置为：点光源S在治疗床下方，探测器平面平行于水平面，点光源到探测器平面中心点的连线SDc始终垂直于探测器平面；
步骤3222：获取X线图像上所有标记点的2D像素坐标Di20(xi2,yi2)，i=1,…,N；
步骤3223：在3D等中心坐标系下，根据等中心C初始3D空间位置C30’(0,0,0)和点光源S到等中心点C的距离dSAD，计算得到点光源S的初始3D空间位置为S30’(-dSAD,0,0)；
步骤3224：在3D等中心坐标系下，根据点光源S 到平板探测器中心点Dc的距离dSID，计算得到探测器中心点Dc的初始3D空间位置为Dc30’(dSID-dSAD,0,0)，以及X线图像上所有标记点的初始3D空间位置为：
Figure QLYQS_8
步骤331：在3D等中心坐标系下计算旋转矩阵；
在3D等中心坐标系下，从C型臂初始位置开始，机架带动术中X光机和探测器围绕等中心点C旋转，依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角，由此计算得到旋转矩阵为：
Figure QLYQS_9
其中L臂、枢轴、C臂旋转轴的初始位置分别为
Figure QLYQS_12
，
Figure QLYQS_13
为术中X光机绕L臂旋转α角所对应的旋转矩阵；绕L臂旋转α角后，枢轴从初始位置
Figure QLYQS_16
旋转到
Figure QLYQS_11
，
Figure QLYQS_14
为术中X光机绕枢轴
Figure QLYQS_18
旋转β角所对应的旋转矩阵；C臂旋转轴先绕L臂旋转α角，再绕枢轴旋转β角后，C臂旋转轴从初始位置
Figure QLYQS_20
旋转到
Figure QLYQS_10
，
Figure QLYQS_15
为术中X光机绕C臂旋转轴
Figure QLYQS_17
旋转γ角所对应的旋转矩阵；旋转矩阵
Figure QLYQS_19
依照如下公式计算：
Figure QLYQS_21
这里
Figure QLYQS_22
表示绕向量
Figure QLYQS_23
旋转θ角的旋转矩阵，将旋转矩阵
Figure QLYQS_24
的参数与
Figure QLYQS_25
的参数一一对应，计算各个矩阵值，得到依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角的旋转矩阵；
步骤332：在3D患者坐标系下计算术中X光机相关各点的第二3D空间位置，包括：
步骤3321：预设3D患者坐标系中等中心点C的第二3D空间位置为C3，其坐标表示为(xc,yc,zc)；
步骤3322：在3D患者坐标系中，经旋转变换后，逐一计算得到：
点光源S的第二3D空间位置变换为S3=RS30’+C3，
探测器中心点Dc的第二3D空间位置变换为Dc3=RDc30’+C3，
X线图像上所有标记点的第二3D空间位置为Di3=RDi30’+C3，i=1,…，N，
其中R为依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角的旋转矩阵，S30’、Dc30’、Di30’分别表示所述点光源S、探测器中心点Dc和X线图像第i个标记点在3D等中心坐标系下的初始3D空间位置。
5.一种投影变换的计算装置，其特征在于，包括：
图像获取模块，用于采用多个标记点对患者不同部位进行体表标记，获取患者的一对术前CT图像与术中X线图像，所述CT图像是由CT切片重建而成的3D图像，所述X线图像则由术中X光机采集；
第一计算模块，用于获取采集所述CT图像时的扫描参数，以及所述CT图像上所有标记点的3D体素坐标，计算得到所述CT图像上所有标记点3D体素坐标的第一3D空间位置；
第二计算模块，用于获取所述术中X光机投影患者生成所述X线图像时的投影已知参数，以及所述X线图像所有标记点的2D像素坐标，计算得到所述X线图像上所有标记点2D像素坐标的第二3D空间位置，具体包括：
步骤31：获取所述术中X光机投影患者生成所述X线图像时的投影已知参数：
从所述术中X线图像的dicom数据及所述术中X光机说明书中，读取投影已知参数，所述投影已知参数包括： L臂旋转角α，枢轴旋转角β，C臂旋转角γ，点光源S到等中心点C的距离dSAD、点光源S 到平板探测器中心点DC的距离dSID、平板探测器尺寸SL*SW、X线图像分辨率RL*RW；
步骤32：确定等中心点、点光源及X线图像标记点的初始3D空间位置，包括：
步骤321：确定3D等中心坐标系和X线图像2D坐标系；
步骤322：计算3D等中心坐标系中X光机相关各点的初始3D空间位置；
步骤33：计算3D患者坐标系中X线图像上所有标记点的第二3D空间位置：
步骤331：在3D等中心坐标系下计算旋转矩阵；
步骤332：在3D患者坐标系下计算X光机相关各点的第二3D空间位置；
第三计算模块，用于依据所述所有标记点的第一3D空间位置和第二3D空间位置，计算所述术中X光机投影患者生成所述X线图像时的投影未知参数，得到配准术前CT图像与术中X线图像的投影变换，具体包括：
步骤41：计算3D患者坐标系中点光源S到CT图像上第i个标记点Qi3的3D空间向量：
Figure QLYQS_26
同时计算3D患者坐标系中点光源S到X线图像上第i个标记点Di3的3D空间向量：
Figure QLYQS_27
其中，R为依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角的旋转矩阵，i=1,…，N，S30’和Di30’分别表示所述点光源S和X线图像上第i个标记点在3D等中心坐标系下的初始3D空间位置，S3、C3、Di3分别表示点光源S、等中心点C、X线图像上第i个标记点在3D患者坐标系下的第二3D空间位置，Qi3则表示CT图像上第i个标记点在3D患者坐标系下的第一3D空间位置，i=1,…,N；
步骤42：计算得到一组投影未知参数xc,yc,zc的值，包括：
CT图像上第i个标记点Qi3和X线图像上第i个标记点Di3为同名点对，根据向量
Figure QLYQS_28
与
Figure QLYQS_29
共线，得到：
Figure QLYQS_30
其中，
Figure QLYQS_31
R为依次绕L臂、枢轴、C臂旋转轴旋转α角、β角和γ角的旋转矩阵，i=1,…,N，S30’和Di30’分别表示所述点光源S和X线图像上第i个标记点在3D等中心坐标系下的初始3D空间位置，Qi3表示CT图像上第i个标记点在3D患者坐标系下的第一3D空间位置；
由此计算得到一组投影未知参数xc,yc,zc的值；
步骤43：计算得到最终投影未知参数xc,yc,zc的值，包括：
任意组合两个标记点在CT与X线图像中的第一和第二3D 空间位置(Qi3，Di3)和(Qj3，Dj3)，i=1,…,N，j=1,…,N，根据所述步骤41和步骤42计算得到一组xc,yc,zc的值；组合多个两两标记点，得到多组xc,yc,zc的值；取所有组xc,yc,zc的平均值，作为最终投影未知参数xc,yc,zc的值。
6.一种计算机设备，其特征在于，包括：存储器和处理器，所述存储器和所述处理器之间互相通信连接，所述存储器中存储有计算机指令，所述处理器通过执行所述计算机指令，从而执行如权利要求1-4任一项所述的方法。
7.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于使计算机执行如权利要求1-4任一项所述的方法。
CN202211113427.3A 2022-09-14 2022-09-14 一种投影变换的计算方法、装置、设备及存储介质 Active CN115205417B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211113427.3A CN115205417B (zh) 2022-09-14 2022-09-14 一种投影变换的计算方法、装置、设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211113427.3A CN115205417B (zh) 2022-09-14 2022-09-14 一种投影变换的计算方法、装置、设备及存储介质

Publications (2)

Publication Number Publication Date
CN115205417A CN115205417A (zh) 2022-10-18
CN115205417B true CN115205417B (zh) 2023-03-14

Family

ID=83573308

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211113427.3A Active CN115205417B (zh) 2022-09-14 2022-09-14 一种投影变换的计算方法、装置、设备及存储介质

Country Status (1)

Country Link
CN (1) CN115205417B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116327228B (zh) * 2023-03-30 2024-04-30 杭州邦杰星医疗科技有限公司 一种2d-3d图像初始值计算的方法
CN117474906B (zh) * 2023-12-26 2024-03-26 合肥吉麦智能装备有限公司 基于脊柱x光图像匹配的术中x光机复位方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112989081A (zh) * 2021-05-20 2021-06-18 首都医科大学附属北京安贞医院 一种数字重建影像库的构建方法和装置
CN113870331A (zh) * 2021-10-07 2021-12-31 浙江大学 一种基于深度学习的胸部ct与x光实时配准算法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20220054199A1 (en) * 2018-12-19 2022-02-24 INDIAN INSTITUTE OF TECHNOLOGY MADRAS (IIT Madras) Robotic surgery systems and surgical guidance methods thereof
CN110533597B (zh) * 2019-03-26 2022-06-28 北京东软医疗设备有限公司 伪影处理及旋转中心确定方法、装置与设备、存储介质
EP3946062A4 (en) * 2019-04-04 2023-06-14 Centerline Biomedical, Inc. SPATIAL REGISTRATION OF A TRACKING SYSTEM WITH AN IMAGE USING TWO-DIMENSIONAL IMAGE PROJECTIONS
CN110148160A (zh) * 2019-05-22 2019-08-20 合肥中科离子医学技术装备有限公司 一种正交x射线影像快速2d-3d医学图像配准方法
CN114821031A (zh) * 2022-02-25 2022-07-29 上海极睿医疗科技有限公司 基于c臂机的术中图像匹配方法、装置和***

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112989081A (zh) * 2021-05-20 2021-06-18 首都医科大学附属北京安贞医院 一种数字重建影像库的构建方法和装置
CN113870331A (zh) * 2021-10-07 2021-12-31 浙江大学 一种基于深度学习的胸部ct与x光实时配准算法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
胸主动脉腔内修复术前CT血管造影与术中X射线图像的配准算法;贾瑞明 等;《中国组织工程研究》;20190919;5658-5663 *

Also Published As

Publication number Publication date
CN115205417A (zh) 2022-10-18

Similar Documents

Publication Publication Date Title
CN115205417B (zh) 一种投影变换的计算方法、装置、设备及存储介质
JP7443462B2 (ja) 治療計画画像、断片内３ｄ画像、及び断片内２ｄｘ線画像の画像見当合わせ
CN107789001B (zh) 一种用于成像扫描的摆位方法和***
US10424118B2 (en) Perspective representation of a virtual scene component
US11576645B2 (en) Systems and methods for scanning a patient in an imaging system
CN106572829B (zh) 医学成像***的扫描区域的定位的控制
JP5497436B2 (ja) 回転ｘ線走査計画システム
US20200268251A1 (en) System and method for patient positioning
CN110533597B (zh) 伪影处理及旋转中心确定方法、装置与设备、存储介质
JP2009022754A (ja) 放射線画像の位置揃えを補正する方法
US11127153B2 (en) Radiation imaging device, image processing method, and image processing program
US20180325472A1 (en) Patient position monitoring system based on 3d surface acquisition technique
CN112168357B (zh) C臂机空间定位模型构建***及方法
CN113298745B (zh) Cta三维重建镜像数据图像投影方法、图像处理方法及装置
US9020215B2 (en) Systems and methods for detecting and visualizing correspondence corridors on two-dimensional and volumetric medical images
TWI657802B (zh) 自定位牙齒電腦斷層掃描方法及裝置
CN114073579B (zh) 手术导航方法、装置、电子设备及存储介质
CN116612161A (zh) 用于医学图像配准的数据集的构建方法和装置
US20210295542A1 (en) Method, radiotherapy device, and computer-readable storage medium for image registration
CN110063741A (zh) 用于x射线成像的x射线成像单元
US20220301199A1 (en) Imaging System
Gopalakrishnan et al. Intraoperative 2d/3d image registration via differentiable x-ray rendering
CN113963056B (zh) Ct图像重建方法、装置、电子设备以及存储介质
US11922645B2 (en) Imaging system
JP7469929B2 (ja) 医用画像処理装置及びｘ線撮影装置

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