CN117392230A - Cbct旋转半径的测量方法、装置、设备及存储介质 - Google Patents

Cbct旋转半径的测量方法、装置、设备及存储介质 Download PDF

Info

Publication number
CN117392230A
CN117392230A CN202210770183.XA CN202210770183A CN117392230A CN 117392230 A CN117392230 A CN 117392230A CN 202210770183 A CN202210770183 A CN 202210770183A CN 117392230 A CN117392230 A CN 117392230A
Authority
CN
China
Prior art keywords
coordinate system
cbct
coordinates
loss function
rotation
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.)
Pending
Application number
CN202210770183.XA
Other languages
English (en)
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.)
Tinavi Medical Technologies Co Ltd
Original Assignee
Tinavi Medical Technologies Co Ltd
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 Tinavi Medical Technologies Co Ltd filed Critical Tinavi Medical Technologies Co Ltd
Priority to CN202210770183.XA priority Critical patent/CN117392230A/zh
Publication of CN117392230A publication Critical patent/CN117392230A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
    • 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/40—Arrangements for generating radiation specially adapted for radiation diagnosis
    • A61B6/4064—Arrangements for generating radiation specially adapted for radiation diagnosis specially adapted for producing a particular type of beam
    • A61B6/4085—Cone-beams
    • 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
    • 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/58—Testing, adjusting or calibrating thereof
    • 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]

Landscapes

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

Abstract

本申请提供了一种CBCT旋转半径的测量方法、装置、设备及存储介质。该方法包括：利用CBCT拍摄标定板得到多个二维图像，获取每个二维图像中的标定点在世界坐标系和像素坐标系下的坐标；建立标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系；基于每个二维图像中的标定点对应的第一转换关系建立第一损失函数，得到内参矩阵和外参矩阵；利用外参矩阵对应的逆矩阵，将每个二维图像对应的焦点位置转换到世界坐标系下，得到焦点位置在世界坐标系下的坐标；基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，得到CBCT旋转半径。本申请可以快速、高效、便捷地获取CBCT旋转半径，提高CBCT旋转半径的测量精度。

Description

CBCT旋转半径的测量方法、装置、设备及存储介质
技术领域
本申请涉及医疗设备技术领域，尤其涉及一种CBCT旋转半径的测量方法、装置、设备及存储介质。
背景技术
锥形束电子计算机断层扫描(Cone-beam-Computed-Tomography，简称CBCT)，由于其摆位灵活、操作方便、比螺旋CT的辐射剂量小，因此被广泛应用于多种手术中，比如辅助引导质子肿瘤治疗、影像导航机器人手术等。在将术前锥形束CT引入到手术过程中，确定CBCT的内参是至关重要的，然而受CBCT的使用过程以及其机械结构的特殊化影响，CBCT的内参并不是固定不变的，并且CBCT拍摄图像容易存在畸变。
目前在将术前用于病情诊断的锥形束CT图像用于规划手术方案并引入到手术导航过程中时，虽然CBCT内参和图像畸变参数可以通过相机标定的方式来获取，但是对于CBCT旋转中心和CBCT旋转半径却没有已知方法来获得。然而，该参数对于CBCT和其他模态融合来说却十分重要。因此，亟需一种能够快速、高效、便捷地获取CBCT旋转半径的方案。
发明内容
有鉴于此，本申请实施例提供了一种CBCT旋转半径的测量方法、装置、设备及存储介质，以解决现有技术存在的无法快速、高效、便捷地获得CBCT旋转半径的问题。
本申请实施例的第一方面，提供了一种CBCT旋转半径的测量方法，包括：利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，并获取每个二维图像中的标定点分别在世界坐标系和像素坐标系下的坐标；建立每个二维图像中的标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系；基于每个二维图像中的标定点对应的第一转换关系建立第一损失函数，基于第一损失函数得到内参矩阵和外参矩阵；利用外参矩阵对应的逆矩阵，将每个二维图像对应的焦点位置转换到世界坐标系下，得到焦点位置在世界坐标系下的坐标；基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于第二损失函数得到CBCT旋转半径。
本申请实施例的第二方面，提供了一种CBCT旋转半径的测量装置，包括：拍摄模块，被配置为利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，并获取每个二维图像中的标定点分别在世界坐标系和像素坐标系下的坐标；建立模块，被配置为建立每个二维图像中的标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系；第一拟合模块，被配置为基于每个二维图像中的标定点对应的第一转换关系建立第一损失函数，基于第一损失函数得到内参矩阵和外参矩阵；转换模块，被配置为利用第一损失函数最小化时的外参矩阵对应的逆矩阵，将每个二维图像对应的焦点位置转换到世界坐标系下，得到焦点位置在世界坐标系下的坐标；第二拟合模块，被配置为基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于第二损失函数得到CBCT旋转半径。
本申请实施例的第三方面，提供了一种电子设备，包括存储器，处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行程序时实现上述方法的步骤。
本申请实施例的第四方面，提供了一种计算机可读存储介质，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现上述方法的步骤。
本申请实施例采用的上述至少一个技术方案能够达到以下有益效果：
通过利利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，并获取每个二维图像中的标定点分别在世界坐标系和像素坐标系下的坐标；建立每个二维图像中的标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系；基于每个二维图像中的标定点对应的第一转换关系建立第一损失函数，基于第一损失函数得到内参矩阵和外参矩阵；利用外参矩阵对应的逆矩阵，将每个二维图像对应的焦点位置转换到世界坐标系下，得到焦点位置在世界坐标系下的坐标；基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于第二损失函数得到CBCT旋转半径。本申请利用CBCT的内参和外参标定过程来直接测量得到CBCT的旋转半径，不仅可以快速、高效、便捷地获得CBCT旋转半径，而且经测量得到的CBCT旋转半径精度高，测量结果的准确性高。
附图说明
为了更清楚地说明本申请实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
图1是本申请实施例提供的CBCT旋转半径的测量方法的流程示意图；
图2是本申请实施例提供的被拍摄物体在CBCT中的成像示意图；
图3是本申请实施例提供的不同焦点位置下的放射源对被拍摄物体进行扫描的示意图；
图4是本申请实施例提供的CBCT旋转半径的测量装置的结构示意图；
图5是本申请实施例提供的电子设备的结构示意图。
具体实施方式
以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本申请实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本申请。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本申请的描述。
随着计算机技术和计算机视觉的发展，影像导航技术越来越多的被应用到外科手术中。影像导航技术基于影像(比如X光、CT等)规划手术方案，并于术中将规划方案利用配准的方式带到术中的患者患处。基于影像导航技术的骨科手术，具有病灶定位准、术中创伤小、手术成功率高等优点。
在微创术式中，一般需要术中在患区固定标记物的情况下，采集患区图像。患区图像不但用于建立图像和患区坐标系之间的映射关系，还用于手术规划。受空间条件和经济条件限制，术中一般采取成像质量不佳的CBCT(即锥形束CT)来采集图像，并且为了尽量减少患者和医生所受辐射剂量，一般拍摄二维片(即二维图像或二维影像)。这无疑加大了医生的规划难度、降低了规划的准确性。因此，将术前用于病情诊断的锥形束CT图像用于规划手术方案并引入到术中导航，不但可以降低医生的规划难度，还能够提高规划精度。
本申请实施例中，CBCT全称为锥形束CT，其原理是利用X线发生器以较低的射线量围绕投照体做环形数字式投照，将围绕投照体多次数字投照后所获得的数据在计算机中重组，进而获得高清三维图像。锥形束CT采用一个X点光源，照射至一块大面积平板探测器之上，整个X光束因而形成一个锥形。
进一步地，在利用CBCT成像时，放射源围绕被扫描物体旋转，其旋转轨迹所处平面称为旋转平面，放射源到CBCT旋转中心的距离即旋转半径，从锥形束的几何关系可见，CBCT成像视野所对应的最大锥角主要由旋转半径来决定。
正如背景技术中描述的内容，在将术前锥形束CT引入到手术过程中，确定CBCT的内参是至关重要的，然而受CBCT用途以及其机械结构的特殊化影响，比如CBCT在使用过程中受自身重力、磨损或外力等影响，CBCT的相机位置及物理形态等会发生变化，因此CBCT的内参并不是固定不变的，并且CBCT拍摄图像容易存在畸变。
目前在将术前用于病情诊断的锥形束CT图像用于规划手术方案并引入到手术导航过程中时，虽然CBCT内参和图像畸变参数可以通过相机标定的方式来获取，但是对于CBCT旋转中心和CBCT旋转半径却没有已知方法来获得。然而，该参数对于CBCT和其他模态融合来说却十分重要，比如在将CBCT二维片和其他模态图像进行配准时可能会使用到CBCT旋转半径的值。因此，亟需一种能够快速、高效、便捷地获取CBCT旋转半径的方案。
鉴于以上现有技术中存在的问题，本申请提供了一种快速、高效、便捷地获取CBCT旋转半径的方法，本申请利用标定板对CBCT的拍摄相机进行标定，并且在相机标定过程中通过建立标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的转换关系，建立与上述转换关系相对应的损失函数；获取损失函数最小化时的内参矩阵和外参矩阵，并利用外参矩阵的逆矩阵，确定每个二维图像所对应焦点位置在世界坐标系下的坐标，最后利用另一损失函数计算最小化时的CBCT旋转半径，得到最终的测量结果。利用本申请实施例可以获得准确的CBCT旋转半径，提升后续利用CBCT旋转半径进行图像配准时的效果。
下面结合附图以及具体实施例对本申请技术方案进行详细说明。
图1是本申请实施例提供的CBCT旋转半径的测量方法的流程示意图。如图1所示，该CBCT旋转半径的测量方法具体可以包括：
S101，利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，并获取每个二维图像中的标定点分别在世界坐标系和像素坐标系下的坐标；
S102，建立每个二维图像中的标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系；
S103，基于每个二维图像中的标定点对应的第一转换关系建立第一损失函数，基于第一损失函数得到内参矩阵和外参矩阵；
S104，利用第一损失函数最小化时的外参矩阵对应的逆矩阵，将每个二维图像对应的焦点位置转换到世界坐标系下，得到焦点位置在世界坐标系下的坐标；
S105，基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于第二损失函数得到CBCT旋转半径。
具体地，本申请实施例中CBCT的相机标定过程是基于2D平面靶标的摄影机标定，通过拍摄多张标定板的图片，再通过多个实际中的标定点(世界坐标中的点)和图片上的标定点(像素坐标中的点)一一对应，即可求出世界坐标和像素坐标的对应关系。
进一步地，本申请实施例的标定板是通过在一个材料密度较低的板子上规则镶嵌一定数目的标定物得到的模具，镶嵌在标定板上的标定物可以采用钢珠等密度大于板子材质的物体。在实际应用中，标定板可以采用棋盘式的布局结构，例如通过在标定板上按照4x4的排列方式布置16个钢珠作为标定物，钢珠在标定板上的覆盖范围可以根据CBCT成像视野进行自由调整，标定物也称为标定点。
进一步地，本申请实施例的二维图像是通过CBCT在不同扫描角度(即不同拍摄角度)下对被拍摄物体(标定板)进行多次拍摄得到的图像，因此，二维图像也可以认为是CBCT所拍摄的二维片或者二维影像。由于CBCT在每次拍摄标定板时相机光心对应的焦点位置均不相同，因此在相机坐标系下，每个二维图像内的标定板在不同拍摄角度下的位姿也会不同。
需要说明的是，本申请实施例中CBCT的相机标定不限于采用张正友标定法，标定板上的标定物也不限于钢珠，标定物在标定板上的排列方式以及数量也可根据实际场景进行调整；应当理解的是，任何能够实现对CBCT的相机标定的方法均适用于本申请，并且其他密度大于标定板所用板材密度的物体均可作为本申请中的标定物，以上标定物在标定板上的排列方式以及标定物的数量也不构成对本申请技术方案的限定。
在一些实施例中，利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，包括：将标定板放置在CBCT的有效转动角度的成像区域内，利用CBCT围绕标定板进行旋转，并以不同的扫描角度拍摄一定数量的二维图像，其中，不同二维图像之间对应的CBCT的射线源的位置不同，标定板中按照等间距设置多个标定物，每个标定物对应一个标定点。
具体地，将模具(即标定板)放置在CBCT有效转动角度的成像区域内，利用CBCT按照不同的拍摄角度采集一定数量的二维图像，不同二维图像之间对应的相机拍摄角度不同，即对应的射线源的位置不同，因此光心的位置也不相同；在实际应用中，CBCT拍摄的二维图像的数量不少于10张，二维图像的数量越多，CBCT旋转半径的测量结果越准确。
进一步地，以空间内的任意点W为原点建立世界坐标系，由于标定板中每个标定物之间的距离是已知的，因此可以获得二维图像中的每个标定点在世界坐标系下的坐标；另外，由于钢珠跟板子的密度差别较大，X光穿过钢珠和板子时的损失不一样，所以感应器探测的信号也就不同，因此钢珠和板子在CBCT成像中出现明显差异，利用阈值分割的方式，可以轻松获取钢珠在每个二维图像的坐标(即在像素坐标系中的坐标值)，得到每个二维图像中钢珠对应的标定点在世界坐标系的坐标以及在像素坐标系的坐标。
在一些实施例中，建立每个二维图像中的标定点在世界坐标系下的坐标与在像素坐标系下的坐标之间的第一转换关系，包括：基于标定点在世界坐标系下的坐标，建立世界坐标系与相机坐标系之间的第二转换关系；利用投影法将标定点在相机坐标系中的坐标转变为图像坐标系中的坐标，以建立相机坐标系与图像坐标系之间的第三转换关系；基于图像坐标系的原点在像素坐标系下的坐标以及像素的尺寸，建立图像坐标系与像素坐标系之间的第四转换关系；将第二转换关系、第三转换关系和第四转换关系进行整合，得到标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系。
具体地，本申请实施例的世界坐标系Ow是以空间内的任意点W为原点建立的坐标系，将上述模具中钢珠在世界坐标系Ow的位置作为标定点在世界坐标系Ow下的坐标；相机坐标系Oc是在相机的角度上衡量物体的坐标系，相机坐标系Oc的原点在相机的光心上，z轴与相机光轴平行；图像坐标系O1用于表征从相机坐标系向图像坐标系的透视投影关系，特点为连续，原点位于相机光轴与成像平面的交点上，图像坐标系O1也可称为成像平面坐标系；像素坐标系Oo是从二维影像中直接得到的信息，二维图像被离散为一个个的像素进行表示，像素坐标系Oo的原点一般位于二维图像的左上角，比如可以将二维图像的左上角中的某个像素点作为原点建立像素坐标系Oo，像素坐标系Oo也可称为屏幕坐标系Oo。
下面在上述建立的坐标系基础上结合附图对本申请实施例中各坐标系之间的转换关系进行详细说明。图2是本申请实施例提供的被拍摄物体在CBCT中的成像示意图；如图2所示，该被拍摄物体在CBCT中的成像具体可以包括：
图2中的立方体表示被拍摄物体，Oc1表示X光的光源发射点，也称为放射源或者发射源等，左上角的示意图表示的是立方体在世界坐标系Ow下的位置，右下角的示意图表示利用放射源对被拍摄物体进行扫描时，被拍摄物体从世界坐标系Ow转换到图像坐标系O1对应的成像平面。中间穿过被拍摄物体的直线表示放射源Oc1与图像坐标系O1的成像平面之间的连线，放射源Oc1围绕该连线中的某一个点(即旋转中心)进行旋转。
进一步地，世界坐标系Ow的坐标可以通过刚体变换(即旋转和平移转换)的方式，转换为相机坐标系Oc的坐标，通过建立世界坐标系与相机坐标系之间的转换关系(即第二转换关系)，将被拍摄物体由世界坐标系Ow变换到相机坐标系Oc，第二转换关系可以表示为：
其中，[X,Y,Z]表示标定点在世界坐标系下的坐标，[Xc,Yc,Zc]表示标定点在相机坐标系下的坐标。
进一步地，将相机坐标系Oc中的坐标通过投影法转换到图像坐标系O1，即通过投影法建立相机坐标系Oc与图像坐标系O1之间的转换关系(即第三转换关系)，第三转换关系可以表示为：
其中，x、y表示标定点投影在图像坐标系中的坐标，f表示焦距，即焦点到平面的距离。
在实际应用中，投影法可以采用透视投影法，透视投影法也称为中心投影法，即利用中心投影法将被拍摄物体投射到投影面上，从而将被拍摄物体在相机坐标系Oc下的坐标变换为在图像坐标系O1下的坐标。
进一步地，像素坐标系u-v的原点为O0，图像坐标系x-y的原点O1在像素坐标系u-v的坐标为(u0，v0)，dx和dy分别表示每个像素在横轴x和纵轴y的物理尺寸；因此基于图像坐标系的原点在像素坐标系下的坐标以及像素的尺寸，可以建立图像坐标系O1与像素坐标系Oo之间的转换关系(第四转换关系)，第四转换关系可以表示为：
在一些实施例中，将第二转换关系、第三转换关系和第四转换关系进行整合得到第一转换关系，包括：利用第三转换关系和第四转换关系生成CBCT的拍摄相机对应的内参矩阵，并获取第二转换关系中CBCT的拍摄相机对应的外参矩阵；基于CBCT的拍摄相机的内参矩阵、外参矩阵以及标定点分别在世界坐标系和像素坐标系下的坐标进行整合得到第一转换关系。
具体地，在将世界坐标系Ow通过刚体变换到相机坐标系Oc，相机坐标系Oc再通过透视投影法转换到图像坐标系O1，图像坐标系O1离散化可得到像素坐标系Oo，通过将上述坐标系之间的转换关系进行整理可得到世界坐标系与像素坐标系之间的转换关系(即第一转换关系)，第一转换关系可以表示为：
其中，第一转换关系中的外参矩阵是从上述第二转换关系中直接得到的参数，第一转换关系中的内参矩阵是通过将第三转换关系与第四转换关系进行相加得到的参数。
在一些实施例中，经过整合得到第一转换关系之后，需要将第一转换关系进一步整合到损失函数(即第一损失函数)里面去，假设模具中有N个钢珠，CBCT总共拍摄M张二维片，则根据第一转换关系可以建立上述的第一损失函数。将基于第一转换关系建立的第一损失函数表示为：
其中，M表示二维图像的个数，N表示标定点的个数，mmn表示标定点在像素坐标系中的坐标，A表示内参矩阵，Rm、Tm构成外参矩阵，Mn表示标定点在世界坐标系中的坐标。
具体地，标定点在像素坐标系中的坐标以及标定点在世界坐标系中的坐标是可以直接获得的信息，通过对第一损失函数执行拟合操作，在参数空间内寻找使上述第一损失函数最小化时的一组参数(即内参和外参)，从而得到第一损失函数最小化时的内参矩阵和外参矩阵。
进一步地，通过以上实施例的相机标定过程，本申请实施例可以获得CBCT的内参以及对应每个拍摄角度的外参(即通过使上述损失函数最小化，得到的内参和外参)；其中，外参描述的是被拍摄物体由世界坐标系到相机坐标系的变换过程，在相机坐标系下，被拍摄物体在不同拍摄角度下的位姿是不同的。
进一步地，利用第一损失函数最小化时的外参矩阵所对应的逆矩阵，将拍摄每个二维图像时对应的焦点位置变换到世界坐标系下，这样就能够知道在世界坐标系内每个焦点相对于被拍摄物体的空间位置，每个焦点与各自成像平面之间的连线所形成的交点就是CBCT旋转中心。
下面结合附图以及具体实施例，对本申请实施例中利用CBCT在不同焦点下对被拍摄物体进行扫描的结果进行说明。图3是本申请实施例提供的不同焦点位置下的放射源对被拍摄物体进行扫描的示意图；如图3所示，该不同焦点位置下的放射源对被拍摄物体进行扫描的结果具体可以包括：
图3中的Oc1、Oc2和Oc3分别对应三个不同的放射源位置(即对应三个不同的焦点位置)，每个放射源与各自成像平面之间具有一条连线，三条连线对应的交点位置即CBCT旋转中心在世界坐标系下的位置，焦点到CBCT旋转中心之间的距离即CBCT旋转半径。由图3可以看出，在利用CBCT对被拍摄物体进行扫描过程中，焦点转动、被拍摄物体不动，且每个焦点位置与各自成像平面之间的连线长度相等。经过对拍摄过程中焦点的转动模拟，发现所有焦点位于一个以旋转中心为球心、以旋转半径为半径的球面上转动。在实际应用中，焦点的位置也可以理解为相机光心的位置，即放射源对应的位置。
在一些实施例中，将基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立的第二损失函数表示为：
ei(x0，y0，z0，r)＝(xi-x0)2+(yi-y0)2+(zi-z0)2-r2
其中，(x0,y0,z0)表示CBCT旋转中心的坐标，r表示CBCT旋转半径，(xi,yi,zi)表示不同扫描角度下焦点位置在世界坐标系的坐标。
具体地，在世界坐标系下，将不同拍摄位姿下焦点的坐标记为[xi,yi,zi]，将旋转中心的坐标记为[x0,y0,z0]，将旋转半径记为r。例如，以图3所示的实际场景中的扫描结果为例，包含三个不同拍摄位姿下的焦点，即Oc1、Oc2和Oc3，可以将每个焦点的坐标记为[x1,y1,z1]、[x2,y2,z2]、[x3,y3,z3]。
在一些实施例中，基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于第二损失函数得到CBCT旋转半径，包括：将第二损失函数进行最小化，利用最小二乘拟合的方式，在参数空间内寻找使第二损失函数的值最小化时的CBCT旋转中心的坐标以及CBCT旋转半径，并将第二损失函数最小化时对应的CBCT旋转半径作为CBCT旋转半径的测量结果。
具体地，在得到每个焦点在世界坐标系下的坐标之后，可以建立上述第二损失函数，利用最小二乘拟合的方式对第二损失函数进行拟合，在参数空间内寻找使上述第二损失函数最小化时的另一组参数(即旋转中心和旋转半径)，从而得到满足第二损失函数最小化时的旋转中心[x0,y0,z0]和旋转半径r。将第二损失函数最小化时的旋转半径r作为最终得到的CBCT旋转半径，从而完成对CBCT旋转半径的测量工作。
在本申请的另一些实施例中，第二损失函数除了采用以上公式进行表示外，还可以通过对上述公式进行变形得到以下公式：
其中，N表示焦点的个数，该公式为上述第二损失函数公式的变形，因此也可以用于拟合得到CBCT旋转半径，并且通过对上述公式中第二损失函数的变量进行整理，可以形成如下所示的式子：
需要说明的是，通过对上述第二损失函数公式进行变形以及变量整理等，可以得到不同表达形式的式子。但是应当理解的是，这些式子仅仅是表达形式上的区别，不同式子之间所表达的含义以及作用是相同的。通过对以上式子进行最小二乘拟合，当式子等号两边的值越接近时，说明旋转半径越接近实际值，因此可以将损失函数最小化时对应的旋转半径作为最终的测量结果。
根据本申请实施例提供的技术方案，本申请通过CBCT的相机标定过程，建立标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的转换关系，并建立与上述转换关系对应的第一损失函数；通过获取第一损失函数最小化时的内参矩阵和外参矩阵，并利用外参矩阵的逆矩阵，确定每个二维图像所对应焦点位置在世界坐标系下的坐标；最后，根据焦点位置的坐标与CBCT旋转中心及CBCT旋转半径之间建立起第二损失函数，对第二损失函数进行最小二乘拟合，求满足第二损失函数最小化时的CBCT旋转半径，从而得到最终的测量结果。通过本申请实施例不仅可以快速、高效、便捷地获取CBCT旋转半径，而且对CBCT旋转半径的测量精度高，提升后续利用CBCT旋转半径进行图像配准时的效果。
下述为本申请装置实施例，可以用于执行本申请方法实施例。对于本申请装置实施例中未披露的细节，请参照本申请方法实施例。
图4是本申请实施例提供的CBCT旋转半径的测量装置的结构示意图。如图4所示，该CBCT旋转半径的测量装置包括：
拍摄模块401，被配置为利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，并获取每个二维图像中的标定点分别在世界坐标系和像素坐标系下的坐标；
建立模块402，被配置为建立每个二维图像中的标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系；
第一拟合模块403，被配置为基于每个二维图像中的标定点对应的第一转换关系建立第一损失函数，基于第一损失函数得到内参矩阵和外参矩阵；
转换模块404，被配置为利用第一损失函数最小化时的外参矩阵对应的逆矩阵，将每个二维图像对应的焦点位置转换到世界坐标系下，得到焦点位置在世界坐标系下的坐标；
第二拟合模块405，被配置为基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于第二损失函数得到CBCT旋转半径。
在一些实施例中，图4的拍摄模块401将标定板放置在CBCT的有效转动角度的成像区域内，利用CBCT围绕标定板进行旋转，并以不同的扫描角度拍摄一定数量的二维图像，其中，不同二维图像之间对应的CBCT的射线源的位置不同，标定板中按照等间距设置多个标定物，每个标定物对应一个标定点。
在一些实施例中，图4的建立模块402基于标定点在世界坐标系下的坐标，建立世界坐标系与相机坐标系之间的第二转换关系；利用投影法将标定点在相机坐标系中的坐标转变为图像坐标系中的坐标，以建立相机坐标系与图像坐标系之间的第三转换关系；基于图像坐标系的原点在像素坐标系下的坐标以及像素的尺寸，建立图像坐标系与像素坐标系之间的第四转换关系；将第二转换关系、第三转换关系和第四转换关系进行整合，得到标定点在世界坐标系下的坐标与像素坐标系下的坐标之间的第一转换关系。
在一些实施例中，图4的建立模块402利用第三转换关系和第四转换关系生成CBCT的拍摄相机对应的内参矩阵，并获取第二转换关系中CBCT的拍摄相机对应的外参矩阵；基于CBCT的拍摄相机的内参矩阵、外参矩阵以及标定点分别在世界坐标系和像素坐标系下的坐标进行整合得到第一转换关系。
在一些实施例中，图4的第一拟合模块403将基于第一转换关系建立的第一损失函数表示为：
其中，M表示二维图像的个数，N表示标定点的个数，mmn表示标定点在像素坐标系中的坐标，A表示内参矩阵，Rm、Tm构成外参矩阵，Mn表示标定点在世界坐标系中的坐标。
在一些实施例中，图4的第二拟合模块405将基于焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立的第二损失函数表示为：
ei(x0，y0，z0，r)＝(xi-x0)2+(yi-y0)2+(zi-z0)2-r2
其中，(x0,y0,z0)表示CBCT旋转中心的坐标，r表示CBCT旋转半径，(xi,yi,zi)表示不同扫描角度下焦点位置在世界坐标系的坐标。
在一些实施例中，图4的第二拟合模块405将第二损失函数进行最小化，利用最小二乘拟合的方式，在参数空间内寻找使第二损失函数的值最小化时的CBCT旋转中心的坐标以及CBCT旋转半径，并将第二损失函数最小化时对应的CBCT旋转半径作为CBCT旋转半径的测量结果。
应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请实施例的实施过程构成任何限定。
图5是本申请实施例提供的电子设备5的结构示意图。如图5所示，该实施例的电子设备5包括：处理器501、存储器502以及存储在该存储器502中并且可以在处理器501上运行的计算机程序503。处理器501执行计算机程序503时实现上述各个方法实施例中的步骤。或者，处理器501执行计算机程序503时实现上述各装置实施例中各模块/单元的功能。
示例性地，计算机程序503可以被分割成一个或多个模块/单元，一个或多个模块/单元被存储在存储器502中，并由处理器501执行，以完成本申请。一个或多个模块/单元可以是能够完成特定功能的一系列计算机程序指令段，该指令段用于描述计算机程序503在电子设备5中的执行过程。
电子设备5可以是桌上型计算机、笔记本、掌上电脑及云端服务器等电子设备。电子设备5可以包括但不仅限于处理器501和存储器502。本领域技术人员可以理解，图5仅仅是电子设备5的示例，并不构成对电子设备5的限定，可以包括比图示更多或更少的部件，或者组合某些部件，或者不同的部件，例如，电子设备还可以包括输入输出设备、网络接入设备、总线等。
处理器501可以是中央处理单元(Central Processing Unit，CPU)，也可以是其它通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其它可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
存储器502可以是电子设备5的内部存储单元，例如，电子设备5的硬盘或内存。存储器502也可以是电子设备5的外部存储设备，例如，电子设备5上配备的插接式硬盘，智能存储卡(Smart Media Card，SMC)，安全数字(Secure Digital，SD)卡，闪存卡(Flash Card)等。进一步地，存储器502还可以既包括电子设备5的内部存储单元也包括外部存储设备。存储器502用于存储计算机程序以及电子设备所需的其它程序和数据。存储器502还可以用于暂时地存储已经输出或者将要输出的数据。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述或记载的部分，可以参见其它实施例的相关描述。
本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
在本申请所提供的实施例中，应该理解到，所揭露的装置/计算机设备和方法，可以通过其它的方式实现。例如，以上所描述的装置/计算机设备实施例仅仅是示意性的，例如，模块或单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通讯连接可以是通过一些接口，装置或单元的间接耦合或通讯连接，可以是电性，机械或其它的形式。
作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本申请各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读存储介质中。基于这样的理解，本申请实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，计算机程序可以存储在计算机可读存储介质中，该计算机程序在被处理器执行时，可以实现上述各个方法实施例的步骤。计算机程序可以包括计算机程序代码，计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。计算机可读介质可以包括：能够携带计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(Read-Only Memory，ROM)、随机存取存储器(Random Access Memory，RAM)、电载波信号、电信信号以及软件分发介质等。需要说明的是，计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如，在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括电载波信号和电信信号。
以上实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的精神和范围，均应包含在本申请的保护范围之内。

Claims (10)

1.一种CBCT旋转半径的测量方法，其特征在于，包括：
利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，并获取每个所述二维图像中的标定点分别在世界坐标系和像素坐标系下的坐标；
建立每个所述二维图像中的标定点在所述世界坐标系下的坐标与所述像素坐标系下的坐标之间的第一转换关系；
基于每个所述二维图像中的标定点对应的所述第一转换关系建立第一损失函数，基于所述第一损失函数得到内参矩阵和外参矩阵；
利用所述外参矩阵对应的逆矩阵，将每个所述二维图像对应的焦点位置转换到世界坐标系下，得到所述焦点位置在世界坐标系下的坐标；
基于所述焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于所述第二损失函数得到CBCT旋转半径。
2.根据权利要求1所述的方法，其特征在于，所述利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，包括：
将所述标定板放置在CBCT的有效转动角度的成像区域内，利用CBCT围绕所述标定板进行旋转，并以不同的扫描角度拍摄一定数量的二维图像，其中，不同所述二维图像之间对应的CBCT的射线源的位置不同，所述标定板中按照等间距设置多个标定物，每个所述标定物对应一个标定点。
3.根据权利要求1所述的方法，其特征在于，所述建立每个所述二维图像中的标定点在所述世界坐标系下的坐标与在所述像素坐标系下的坐标之间的第一转换关系，包括：
基于所述标定点在世界坐标系下的坐标，建立所述世界坐标系与所述相机坐标系之间的第二转换关系；
利用投影法将所述标定点在所述相机坐标系中的坐标转变为图像坐标系中的坐标，以建立所述相机坐标系与所述图像坐标系之间的第三转换关系；
基于所述图像坐标系的原点在像素坐标系下的坐标以及像素的尺寸，建立所述图像坐标系与所述像素坐标系之间的第四转换关系；
将所述第二转换关系、所述第三转换关系和所述第四转换关系进行整合，得到所述标定点在所述世界坐标系下的坐标与所述像素坐标系下的坐标之间的第一转换关系。
4.根据权利要求3所述的方法，其特征在于，将所述第二转换关系、所述第三转换关系和所述第四转换关系进行整合得到所述第一转换关系，包括：
利用所述第三转换关系和所述第四转换关系生成所述CBCT的拍摄相机对应的内参矩阵，并获取所述第二转换关系中所述CBCT的拍摄相机对应的外参矩阵；
基于所述CBCT的拍摄相机的内参矩阵、外参矩阵以及所述标定点分别在所述世界坐标系和所述像素坐标系下的坐标进行整合得到所述第一转换关系。
5.根据权利要求1所述的方法，其特征在于，将基于所述第一转换关系建立的所述第一损失函数表示为：
其中，M表示二维图像的个数，N表示标定点的个数，mmn表示标定点在像素坐标系中的坐标，A表示内参矩阵，Rm、Tm构成外参矩阵，Mn表示标定点在世界坐标系中的坐标。
6.根据权利要求1所述的方法，其特征在于，将基于所述焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立的所述第二损失函数表示为：
ei(x0，y0，z0，r)＝(xi-x0)2+(yi-y0)2+(zi-z0)2-r2
其中，(x0,y0,z0)表示CBCT旋转中心的坐标，r表示CBCT旋转半径，(xi,yi,zi)表示不同扫描角度下焦点位置在世界坐标系的坐标。
7.根据权利要求1所述的方法，其特征在于，基于所述焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于所述第二损失函数得到CBCT旋转半径，包括：
将所述第二损失函数进行最小化，利用最小二乘拟合的方式，在参数空间内寻找使所述第二损失函数的值最小化时的CBCT旋转中心的坐标以及CBCT旋转半径，并将所述第二损失函数最小化时对应的CBCT旋转半径作为CBCT旋转半径的测量结果。
8.一种CBCT旋转半径的测量装置，其特征在于，包括：
拍摄模块，被配置为利用CBCT对放置在成像区域内的标定板进行旋转拍摄，得到不同扫描角度下拍摄的多个二维图像，并获取每个所述二维图像中的标定点分别在世界坐标系和像素坐标系下的坐标；
建立模块，被配置为建立每个所述二维图像中的标定点在所述世界坐标系下的坐标与所述像素坐标系下的坐标之间的第一转换关系；
第一拟合模块，被配置为基于每个所述二维图像中的标定点对应的所述第一转换关系建立第一损失函数，基于所述第一损失函数得到内参矩阵和外参矩阵；
转换模块，被配置为利用所述第一损失函数最小化时的外参矩阵对应的逆矩阵，将每个所述二维图像对应的焦点位置转换到世界坐标系下，得到所述焦点位置在世界坐标系下的坐标；
第二拟合模块，被配置为基于所述焦点位置在世界坐标系下的坐标、CBCT旋转中心以及CBCT旋转半径建立第二损失函数，基于所述第二损失函数得到CBCT旋转半径。
9.一种电子设备，包括存储器，处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现权利要求1至7中任一项所述的方法。
10.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1至7中任一项所述的方法。
CN202210770183.XA 2022-06-30 2022-06-30 Cbct旋转半径的测量方法、装置、设备及存储介质 Pending CN117392230A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210770183.XA CN117392230A (zh) 2022-06-30 2022-06-30 Cbct旋转半径的测量方法、装置、设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210770183.XA CN117392230A (zh) 2022-06-30 2022-06-30 Cbct旋转半径的测量方法、装置、设备及存储介质

Publications (1)

Publication Number Publication Date
CN117392230A true CN117392230A (zh) 2024-01-12

Family

ID=89470756

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210770183.XA Pending CN117392230A (zh) 2022-06-30 2022-06-30 Cbct旋转半径的测量方法、装置、设备及存储介质

Country Status (1)

Country Link
CN (1) CN117392230A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117934885A (zh) * 2024-03-21 2024-04-26 沈阳市口腔医院 基于cbct数据与三维面扫数据的坐标自动匹配方法及系统

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117934885A (zh) * 2024-03-21 2024-04-26 沈阳市口腔医院 基于cbct数据与三维面扫数据的坐标自动匹配方法及系统

Similar Documents

Publication Publication Date Title
US9561387B2 (en) Ambiguity-free optical tracking system
US7620144B2 (en) Parallel stereovision geometry in image-guided radiosurgery
Dey et al. Automatic fusion of freehand endoscopic brain images to three-dimensional surfaces: creating stereoscopic panoramas
US8121380B2 (en) Computerized imaging method for a three-dimensional reconstruction from two-dimensional radiological images; implementation device
US9715739B2 (en) Bone fragment tracking
US11257241B2 (en) System and method for component positioning by registering a 3D patient model to an intra-operative image
CN111132730B (zh) 与放射治疗设备一起使用的患者监测系统的校准方法
CN108245788B (zh) 一种双目测距装置及方法、包括该装置的加速器放疗系统
CN112258593B (zh) 单目相机下ct或pet-ct智能定位扫描方法
US20180345040A1 (en) A target surface
Lin et al. A novel approach of surface texture mapping for cone-beam computed tomography in image-guided surgical navigation
CN112401919B (zh) 一种基于摆位模型的辅助摆位方法及系统
Meng et al. An automatic markerless registration method for neurosurgical robotics based on an optical camera
Xu et al. Design and validation of a spinal surgical navigation system based on spatial augmented reality
Guéziec et al. Providing visual information to validate 2-D to 3-D registration
Selby et al. Patient positioning with X-ray detector self-calibration for image guided therapy
CN117392230A (zh) Cbct旋转半径的测量方法、装置、设备及存储介质
Jenkins et al. Using a handheld stereo depth camera to overcome limited field‐of‐view in simulation imaging for radiation therapy treatment planning
COMLEKCILER et al. Three-dimensional repositioning of jaw in the orthognathic surgery using the binocular stereo vision
Chen et al. An augmented reality platform for planning of minimally invasive cardiac surgeries
JP2000020696A (ja) 医用画像合成装置
CN112472293A (zh) 一种术前三维影像与术中透视图像的配准方法
Corbi et al. Joint calibration of RGB and X-ray cameras
Zhao et al. Augmented Reality Calibration with Stereo Image Registration for Surgical Navigation
Bertelsen et al. Proof‐of‐concept of a robotic‐driven photogrammetric scanner for intra‐operative knee cartilage repair

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