CN114842154B - 一种基于二维x射线图像重建三维影像的方法和系统 - Google Patents

一种基于二维x射线图像重建三维影像的方法和系统 Download PDF

Info

Publication number
CN114842154B
CN114842154B CN202210777073.6A CN202210777073A CN114842154B CN 114842154 B CN114842154 B CN 114842154B CN 202210777073 A CN202210777073 A CN 202210777073A CN 114842154 B CN114842154 B CN 114842154B
Authority
CN
China
Prior art keywords
dimensional
image
transformation
ray
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
CN202210777073.6A
Other languages
English (en)
Other versions
CN114842154A (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.)
Jiangsu Jicui Sukesi Technology Co ltd
Original Assignee
Jiangsu Jicui Sukesi Technology 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 Jiangsu Jicui Sukesi Technology Co ltd filed Critical Jiangsu Jicui Sukesi Technology Co ltd
Priority to CN202210777073.6A priority Critical patent/CN114842154B/zh
Publication of CN114842154A publication Critical patent/CN114842154A/zh
Application granted granted Critical
Publication of CN114842154B publication Critical patent/CN114842154B/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
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00—Manipulating 3D models or images for computer graphics
    • G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Architecture (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及一种基于二维X射线图像重建三维影像的方法和系统。主要是获取包含目标结构的三维模板影像和至少两幅不同方向拍摄的待重建目标结构的X射线图，将各个X射线图校准对齐到同一参考空间坐标系，基于各个方向的X射线图和三维模板影像，将三维模板影像初始对齐到参考空间坐标系中，基于各个方向的X射线图，对三维模板影像通过迭代的方式进行若干次二维/三维相似变换配准，获得初级三维重建影像，基于各个方向的X射线图，对初级三维重建影像通过迭代的方式进行若干次包含畸变校正的二维/三维非线性变换配准，获得最终三维重建影像。该方法重建后的三维影像外轮廓形状精确度较高，且能够较好地重建内部的重要结构。

Description

一种基于二维X射线图像重建三维影像的方法和系统
技术领域
本发明涉及图像处理领域，特别是涉及基于二维X射线图像重建三维影像的方法和系统。
背景技术
临床上，在制定手术计划前，通常会需要对相关部位拍摄二维X射线图像或三维计算机断层扫描影像（CT影像）。在拍摄X射线图像时，人体所吸收的辐射剂量较小，然而由于X射线图像只有平面投影图像数据（丢失了透射方向上的信息），医生在制定方案时需要通过间接估计来做出决策。相较而言， CT影像可构建出精确的三维结构信息，但辐射剂量带来的身体损伤也更大。
因此，临床上渴望一种介于两者之间的技术手段，即在不增加人体辐射伤害的情况下能够满足临床术前计划需求，但现有的技术都是通过二维X射线图像重建三维影像的外部轮廓结构，缺乏内部结构信息。
发明内容
基于此，提供一种基于二维X射线图像重建三维影像的方法。该方法重建后的三维影像外轮廓形状精确度较高，且能够较好地重建内部的轮廓结构。
一种基于二维X射线图像重建三维影像的方法，包括：
获取包含目标结构的三维模板影像和至少两幅不同方向拍摄的待重建目标结构的X射线图，将各个X射线图校准对齐到同一参考空间坐标系，
基于各个方向的X射线图和三维模板影像，将三维模板影像初始对齐到参考空间坐标系中，
基于各个方向的X射线图，对三维模板影像通过迭代的方式进行若干次二维/三维相似变换（带缩放因子的刚性变换）配准，获得初级三维重建影像，
基于各个方向的X射线图，对初级三维重建影像通过迭代的方式进行若干次包含畸变校正的二维/三维非线性变换配准，获得最终三维重建影像。
本申请由于使用的三维模板影像不是由患者直接拍摄，而是利用以往的已经采集拍摄的CT影像来制作获得，患者只需拍摄至少2张X射线图，这样可在不增加人体辐射伤害的情况下满足临床术前计划需求。而且，该方法重建后的三维影像外轮廓形状精确度较高。且本申请由于使用的是三维模板影像，其包含了内部结构信息，所以重建后获得的最终三维重建影像的能够较好地重建内部重要结构。
在其中一个实施例中，所述三维模板影像是基于一组CT影像获得。
在其中一个实施例中，所述三维模板影像是基于一组CT影像获得，具体包括：
a.对一组CT影像进行三维影像分割预处理，获得只包含目标结构的三维体素影像数据集
Figure 229240DEST_PATH_IMAGE001
，
b.对三维体素影像数据集
Figure 100002_DEST_PATH_IMAGE002
进行影像配准操作，具体包括：
(1)从三维体素影像数据集
Figure 150929DEST_PATH_IMAGE002
中任选一个作为参考影像
Figure 759765DEST_PATH_IMAGE003
，利用三维仿射变换配准方法把三维体素影像数据集中的其它各个影像与参考影像
Figure 120208DEST_PATH_IMAGE003
对齐，得到相对于参考影像
Figure 100002_DEST_PATH_IMAGE004
消除了位置和线性几何差异的影像数据集
Figure DEST_PATH_IMAGE005
，
(2)基于三维B样条自由形变配准方法，把每个影像数据
Figure 100002_DEST_PATH_IMAGE006
配准到参考影像
Figure 255523DEST_PATH_IMAGE003
上，得到三维B样条变换参数集合
Figure 161162DEST_PATH_IMAGE007
和一组自由形变后形貌对准的目标结构的影像数据集合
Figure 100002_DEST_PATH_IMAGE008
，
(3)对影像数据集合
Figure 444245DEST_PATH_IMAGE008
与三维B样条变换参数集合
Figure 827953DEST_PATH_IMAGE007
分别求平均，得到均值影像
Figure 145670DEST_PATH_IMAGE009
和三维B样条均值变换参数
Figure 100002_DEST_PATH_IMAGE010
，
(4)计算三维B样条均值变换参数
Figure 222211DEST_PATH_IMAGE010
的逆变换
Figure 992590DEST_PATH_IMAGE011
来重采样均值影像
Figure 242305DEST_PATH_IMAGE009
，以产生三维体素影像数据
Figure 100002_DEST_PATH_IMAGE012
，
Figure 182291DEST_PATH_IMAGE013
，
c.对三维体素影像数据
Figure 100002_DEST_PATH_IMAGE014
进行处理来生成三维模板影像，具体包括：
(1)从三维体素影像数据
Figure 429732DEST_PATH_IMAGE014
获取目标结构的外轮廓二值掩模，将目标结构的外轮廓二值掩模乘以三维体素影像数据
Figure 687407DEST_PATH_IMAGE014
得到独立的目标结构的影像数据，
(2)从三维体素影像数据
Figure 678497DEST_PATH_IMAGE014
获取目标结构的内部结构二值掩模，基于外轮廓二值掩模和内部结构二值掩模产生三维表面形状网格模型，在三维表面形状网格模型上确定三维解剖结构标记点集，
(3)对三维表面形状网格模型进行稀疏化处理，将稀疏化处理后的目标结构的外表面和内表面的网格点合并，得到目标结构的内外表面点集，
所述三维模板影像包括：目标结构的三维影像数据、目标结构的三维解剖结构标记点集、目标结构的内外表面点集。
在其中一个实施例中，所述基于各个X射线图和三维模板影像，将三维模板影像初始对齐到参考空间坐标系中，具体包括：
在X射线图像上获取目标结构的二维解剖结构标记点集，将二维解剖结构标记点集做反投影，计算得到参考空间下X射线图像上解剖结构标记点集的三维位置，将反投影计算得到的三维解剖结构标记点集和三维模板影像中定义的三维解剖结构标记点集做配对点云刚性配准，获得三维刚性变换，用该三维刚性变换将三维模板影像初始对齐到参考空间坐标系中。
在其中一个实施例中，基于各个X射线图，对三维模板影像通过迭代的方式进行若干次二维/三维相似变换（带缩放因子的刚性变换）配准，获得初级三维重建影像，具体包括：
a.基于拍摄各X射线图的外参投影（矩阵）参数，从三维模板影像投影生成各个对应投影方向上的数字重建射线（DRR）图像，并将三维模板影像中目标结构的内外表面点集投影到各个DRR图像所在平面上，获得对应的表面投影点集；
b.计算每个投影方向的DRR图像与对应X射线图的二维相似变换图像配准，得到二维相似变换，用二维相似变换更新对应的表面投影点集的位置，获得各个投影方向的表面投影点集的新位置；
c.基于各个投影方向上表面投影点集的新位置进行反投影，获得目标结构的新的内外表面点集，目标结构的新的内外表面点集与旧的内外表面点集做配对点云配准，得到三维带缩放因子的刚性变换；
d.用该三维带缩放因子的刚性变换去重采样三维模板影像，重采样后的三维模板影像作为下次迭代的输入，
e.经过一次或多次迭代后，获得初级三维重建影像。
在其中一个实施例中，所述基于各个X射线图，对初级三维重建影像通过迭代的方式进行若干次包含畸变校正的二维/三维非线性变换配准，获得最终三维重建影像，具体包括：
a.基于拍摄各X射线图的外参投影（矩阵）参数，从初级三维重建影像投影生成各个对应投影方向上的DRR图像，并将目标结构的内外表面点集投影到各个DRR图像所在平面上，获得对应的表面投影点集；
b.计算每个投影方向的DRR图像与对应X射线图的二维B样条自由形变图像配准，得到二维B样条变换，用二维B样条变换更新对应的表面投影点集的位置，获得各个表面投影点集的新位置；
c.基于各个表面投影点集的新位置进行反投影，获得目标结构的新的内外表面点集，目标结构的新的内外表面点集与旧的内外表面点集之间计算薄板样条变换；
d.基于薄板样条变换，计算三维形变场，对该三维形变场做畸变校正，得到校正后的三维形变场，利用校正后的三维形变场去重采样初级三维重建影像，重采样后的初级三维重建影像作为下次迭代的输入，
e.经过一次或多次迭代后，获得最终三维重建影像。
一种图像处理系统，包括图像获取模块和图像处理模块，所述图像获取模块用于获取包含目标结构的三维模板影像和至少两幅不同方向拍摄的待重建目标结构的X射线图，
所述图像处理模块用于按照所述的基于二维X射线图像重建三维影像的方法对三维模板影像进行处理。
在其中一个实施例中，所述目标结构为股骨和胫骨。
一种计算机装置，包括：处理器、存储器、通信接口和通信总线，所述处理器、存储器和通信接口通过所述通信总线完成相互间的通信，所述存储器用于存放至少一个可执行指令，所述可执行指令使所述处理器执行所述的基于二维X射线图像重建三维影像的方法对应的操作。
一种存储介质，所述计算机存储介质中存储有至少一个可执行指令，所述可执行指令使处理器执行所述的基于二维X射线图像重建三维影像的方法对应的操作。
附图说明
图1为本申请的实施例的基于二维X射线图像重建三维影像的方法的流程图。
图2为本申请的实施例的股骨和胫骨结构的三维模板影像创建流程图。
图3为本申请的实施例的三维到二维投影与二维到三维反投影的示意图。
图4为本申请的实施例的基于三角测量法求两幅X射线图的二维解剖结构标记点的反投影交点的示意图。
图5为本申请的实施例的基于表面点集进行投影/反投影的配对点相似变换（带缩放因子的刚性变换）配准的示意图。
图6为本申请的实施例的断层截面展示图。按照从左向右的顺序，第一幅图是薄板样条变换生成的三维自由形变场的示意图。第二幅图是从该三维自由形变场建立B样条网格，B样条网格控制点的形变示意图。第三幅图是经过正则化约束后的B样条网格控制点形变的示意图。第四幅图是从校正后的B样条变换产生的三维自由形变场的示意图。
图7为本申请的实施例的基于二维X射线图像重建三维影像的方法在不同阶段获得的效果图。第一排是不同阶段根据本申请的方法得到的股骨和胫骨的DRR图像，其中，第一幅图是将三维模板影像初始对齐到参考空间坐标系中投影得到的DRR图像，第二幅图是基于各X射线图，对三维模板影像通过迭代的方式进行三次二维/三维相似变换配准，获得的初级三维重建影像投影得到的DRR图像，第三幅图是基于各个X射线图，对初级三维重建影像通过迭代的方式进行若干次（这里只迭代了一次）包含畸变校正的二维/三维非线性变换配准，获得的最终三维重建影像投影得到的DRR图像。第二排是将第一排的图像抽取得到轮廓图，并将轮廓图叠加到患者下肢X射线图后得到的图像。
具体实施方式
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似改进，因此本发明不受下面公开的具体实施例的限制。
需要说明的是，当元件被称为“固定于”另一个元件，它可以直接在另一个元件上或者也可以存在居中的元件。当一个元件被认为是“连接”另一个元件，它可以是直接连接到另一个元件或者可能同时存在居中元件。
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
如图1所示，本申请的实施例提供了一种基于二维X射线图像重建三维影像的方法，该方法具体包括：
获取包含目标结构的三维模板影像和至少两幅不同方向拍摄的待重建目标结构的X射线图。上述不同方向拍摄的目标结构的X射线图也就是从不同拍摄角度获得的同一目标结构的X射线图。
接着，将各个X射线图校准对齐到同一参考空间坐标系。为了校准对齐到同一参考空间坐标系，可进行校准和空间变换链计算，这些都可应用现有的方法实现。
基于各个方向的X射线图和三维模板影像，将三维模板影像初始对齐到参考空间坐标系中。
基于各个方向的X射线图，对三维模板影像通过迭代的方式进行若干次二维/三维相似变换（即带缩放因子的刚性变换）配准，获得初级三维重建影像。
基于各个方向的X射线图，对初级三维重建影像通过迭代的方式进行若干次包含畸变校正的二维/三维非线性变换配准，获得最终三维重建影像。
具体的，上述X射线图可为患者直接拍摄获得，待重建的目标结构为患者的特定人体结构。例如，该目标结构为患者腿部的股骨和胫骨，患者拍摄两幅或多幅下肢X射线图，每幅X射线图中都包含股骨和胫骨结构。可以理解，上述目标结构也可以是人体的其它骨头结构。
具体的，上述三维模板影像可基于以往已经采集拍摄的CT影像获得。也就是可使用其他人已经拍摄好的三维影像来获得上述三维模板影像。上述三维模板影像的目标结构为与待重建的目标结构对应的特定人体结构。例如，该目标结构为腿部的股骨和胫骨结构。
可以理解，所述三维模板影像可以是基于一组CT影像获得，也可以是基于其它包含三维数据信息的影像数据获得。
以下通过股骨和胫骨结构为例，如图2所示,详细介绍如何通过一组CT影像获得三维模板影像。
a.对一组K个单侧下肢CT影像（包含股骨和胫骨）进行三维影像分割预处理，获得只包含目标结构的三维体素影像数据集
Figure 705228DEST_PATH_IMAGE001
，具体包括：
（1）将其中不同侧（左侧/右侧）的CT影像做镜像翻转，使各个影像都处于（解剖几何）同侧。
（2）对每个CT影像实施图像分割操作，提取股骨和胫骨结构，得到股骨和胫骨区域二值掩模（1-目标体素/0-背景体素）。以上操作可通过在CT断层影像上由人工勾画实现，也可以采用基于机器学习或深度学习的三维医学影像分割方法来实现。
（3）将提取出的股骨和胫骨区域二值掩模乘以原始的下肢CT影像得到只包含股骨和胫骨结构的干净三维体素影像数据集
Figure 123571DEST_PATH_IMAGE001
。
b.对三维体素影像数据集
Figure 416012DEST_PATH_IMAGE002
进行影像配准操作，以此来实现空间和形貌的对齐，具体包括：
(1)从三维体素影像数据集
Figure 460060DEST_PATH_IMAGE002
中任选一个作为参考影像
Figure 154347DEST_PATH_IMAGE003
，利用三维仿射变换配准方法把三维体素影像数据集中的其它各个影像与参考影像
Figure 743591DEST_PATH_IMAGE003
对齐，通过重采样，得到相对于参考影像
Figure 975858DEST_PATH_IMAGE004
消除了位置和线性几何差异的影像数据集
Figure 308750DEST_PATH_IMAGE005
。
(2)基于三维B样条自由形变配准方法，把每个影像数据
Figure 44494DEST_PATH_IMAGE006
配准到参考影像
Figure 132536DEST_PATH_IMAGE003
上，得到三维B样条变换参数集合
Figure 337252DEST_PATH_IMAGE007
和一组自由形变后形貌对准的股骨和胫骨的影像数据集合
Figure 270573DEST_PATH_IMAGE008
。这里
Figure 126403DEST_PATH_IMAGE015
，
Figure 100002_DEST_PATH_IMAGE016
。
(3)对影像数据集合
Figure 509979DEST_PATH_IMAGE008
与三维B样条变换参数集合
Figure 998730DEST_PATH_IMAGE007
分别求平均，得到均值影像
Figure 939004DEST_PATH_IMAGE009
和三维B样条均值变换参数
Figure 655199DEST_PATH_IMAGE010
。
(4)由于任意选择的参考影像
Figure 819464DEST_PATH_IMAGE003
带有偏向性，计算均值变换参数
Figure 998773DEST_PATH_IMAGE010
的逆变换
Figure 726426DEST_PATH_IMAGE011
来重采样均值影像
Figure 42001DEST_PATH_IMAGE009
，以产生相对该组无偏向的三维体素影像数据
Figure 642747DEST_PATH_IMAGE012
，
Figure 293040DEST_PATH_IMAGE013
。这里采用现有的求三维B样条变换逆变换的方法来近似计算
Figure 575116DEST_PATH_IMAGE011
，使得以下公式1的代价函数值最小。
Figure 73094DEST_PATH_IMAGE017
（1） 。
c.对三维体素影像数据
Figure 703795DEST_PATH_IMAGE014
进行处理来生成三维模板影像，具体包括：
(1)从三维体素影像数据
Figure 106964DEST_PATH_IMAGE014
获取股骨和胫骨的外轮廓二值掩模
Figure 100002_DEST_PATH_IMAGE018
和
Figure 927152DEST_PATH_IMAGE019
，将股骨和胫骨的外轮廓二值掩模分别乘以三维体素影像数据
Figure 466587DEST_PATH_IMAGE014
得到独立的股骨和胫骨的影像数据
Figure 100002_DEST_PATH_IMAGE020
和
Figure 268190DEST_PATH_IMAGE021
。
(2)从三维体素影像数据
Figure 643807DEST_PATH_IMAGE014
获取股骨和胫骨的内部结构二值掩模
Figure 100002_DEST_PATH_IMAGE022
和
Figure 516954DEST_PATH_IMAGE023
。其中的内部结构也就是骨髓腔结构。具体的，由人工通过工具一层层地从
Figure 723945DEST_PATH_IMAGE014
的断层图像上分别勾画股骨和胫骨内的骨壁结构，得到股骨和胫骨的骨髓腔结构的二值掩模
Figure 775077DEST_PATH_IMAGE022
和
Figure 881399DEST_PATH_IMAGE023
。
接着，基于外轮廓二值掩模
Figure 43390DEST_PATH_IMAGE018
和
Figure 557417DEST_PATH_IMAGE019
以及内部结构二值掩模
Figure 841768DEST_PATH_IMAGE022
和
Figure 926399DEST_PATH_IMAGE023
产生三维表面形状网格模型，该步骤可由现有的三维表面网格生成工具实现。在三维表面形状网格模型上确定三维解剖结构标记点集
Figure 100002_DEST_PATH_IMAGE024
和
Figure 141348DEST_PATH_IMAGE025
，该步骤由人工选择。
(3)对三维表面形状网格模型进行稀疏化处理，将稀疏化处理后的同一骨结构的外表面和骨髓腔表面的网格点合并，得到股骨的内外表面点集
Figure 100002_DEST_PATH_IMAGE026
和胫骨的内外表面点集
Figure 447565DEST_PATH_IMAGE027
。
按照上述步骤生成的三维模板影像包括以下数据：股骨和胫骨影像数据：
Figure 840500DEST_PATH_IMAGE020
和
Figure 740323DEST_PATH_IMAGE021
。股骨模型：
Figure 493384DEST_PATH_IMAGE018
，
Figure 467156DEST_PATH_IMAGE024
，
Figure 280260DEST_PATH_IMAGE026
。胫骨模型：
Figure 401800DEST_PATH_IMAGE019
，
Figure 709285DEST_PATH_IMAGE025
，
Figure 52410DEST_PATH_IMAGE027
。
需要说明的是，本申请的实施例中，为了重建出患者下肢股骨和胫骨的三维形貌体积影像，需要至少两幅经过校准并对齐到同一参考空间坐标系下的患者下肢X射线图。经过校准每一幅X射线图都有一个3
Figure 100002_DEST_PATH_IMAGE028
4大小的投影矩阵
Figure 787148DEST_PATH_IMAGE029
，用于从三维模板影像投影生成二维数字重建射线图像以及将二维成像平面点反投影回到三维空间。投影和反投影过程参见图3，从中可推导出
Figure 582934DEST_PATH_IMAGE029
的表达式以及投影/反投影的计算方式。
具体地，对于参考空间坐标系下一点
Figure DEST_PATH_IMAGE030
，它在以
Figure 621341DEST_PATH_IMAGE031
为原点的三维投影坐标系下记为
Figure DEST_PATH_IMAGE032
。这里，R为从参考空间坐标系到投影坐标系的旋转矩阵，
Figure 569705DEST_PATH_IMAGE032
可通过公式2所示的齐次坐标变换得到。
Figure 724612DEST_PATH_IMAGE033
(2)。
定义二维投影平面坐标系，原点位于平面左上角，记为
Figure DEST_PATH_IMAGE034
，该原点在三维投影坐标系下记为
Figure 945378DEST_PATH_IMAGE035
。
Figure 594665DEST_PATH_IMAGE032
点投影到该成像平面上的点记为
Figure DEST_PATH_IMAGE036
，该点位置在三维投影坐标系下记为
Figure 912383DEST_PATH_IMAGE037
。平面点
Figure 988923DEST_PATH_IMAGE036
的坐标可通过公式3计算得到，其中
Figure DEST_PATH_IMAGE038
和
Figure 431406DEST_PATH_IMAGE039
为像素物理尺寸，而根据比例关系
Figure DEST_PATH_IMAGE040
，公式3可转换成公式4。
Figure 133651DEST_PATH_IMAGE041
(3)。
Figure DEST_PATH_IMAGE042
(4)。
给定
Figure 243559DEST_PATH_IMAGE043
，结合公式4可以推导得到公式5：
Figure DEST_PATH_IMAGE044
(5)。
结合公式2和5，得到如下投影公式6：
Figure 491000DEST_PATH_IMAGE045
(6)。
这里定义
Figure DEST_PATH_IMAGE046
，
Figure 420779DEST_PATH_IMAGE047
，则投影矩阵
Figure DEST_PATH_IMAGE048
。
为了得到参考空间坐标系下任意一点
Figure 678715DEST_PATH_IMAGE030
对应的投影像素坐标，采用如下步骤：
步骤1、首先计算中间变量
Figure 456178DEST_PATH_IMAGE049
。
步骤2、再计算投影平面上的像素位置
Figure DEST_PATH_IMAGE050
。
对于反投影计算，需要确定两点要素：1）投射焦点，即三维投影坐标系的原点在参考空间坐标系下的位置
Figure 123788DEST_PATH_IMAGE031
；2）投影平面每一点像素的反投影方向向量，记为
Figure 619492DEST_PATH_IMAGE051
。这里定义矩阵
Figure DEST_PATH_IMAGE052
，
Figure 397961DEST_PATH_IMAGE053
，则可推导出以下关系式：
Figure DEST_PATH_IMAGE054
，投射焦点
Figure 295510DEST_PATH_IMAGE055
。从
Figure DEST_PATH_IMAGE056
开始推导，
Figure 71705DEST_PATH_IMAGE057
，而由于
Figure 851442DEST_PATH_IMAGE050
，则投影平面上任意一点像素的反投影方向向量为
Figure DEST_PATH_IMAGE058
。
基于以上三维模板影像和投影计算定义，按图1所示流程可对患者下肢股骨和胫骨结构三维形貌体积影像进行重建。具体包括以下内容。
在上述基础上，所述基于各个X射线图和三维模板影像，将三维模板影像初始对齐到参考空间坐标系中，具体包括：在X射线图像上获取目标结构的二维解剖结构标记点集，将二维解剖结构标记点集做反投影，计算得到参考空间下X射线图像上解剖结构标记点集的三维位置，将反投影计算得到的三维解剖结构标记点集和三维模板影像中定义的三维解剖结构标记点集做配对点云刚性配准，获得三维刚性变换，用该三维刚性变换将三维模板影像初始对齐到参考空间坐标系中。
具体的，根据前述创建三维模板影像时定义的股骨和胫骨的解剖结构标记点集
Figure 371285DEST_PATH_IMAGE059
和
Figure DEST_PATH_IMAGE060
，在每幅X射线照射图像上标记出对应标记点的像素坐标，分别记为
Figure 372608DEST_PATH_IMAGE061
和
Figure DEST_PATH_IMAGE062
，这里的
Figure 132753DEST_PATH_IMAGE063
指的是第
Figure 586737DEST_PATH_IMAGE063
个方向的图。
通过前述推导的反投影计算公式以及图4所示的三角测量法，可得到参考空间坐标系下的患者股骨和胫骨解剖结构标记点集
Figure DEST_PATH_IMAGE064
和
Figure 723320DEST_PATH_IMAGE065
。
三角测量法用于求两条射线的交点，在本申请中，由于存在图像噪声和校准误差，两条反射线通常不会完全相交，可首先计算两条射线间距离最短时对应的空间直线，该直线与两条射线交点的中点位置即为反投影交点。对于两幅以上X射线图的情景，可以计算两两相邻图像间的反投影交点，再取均值。
基于两组配对点集（股骨）
Figure 579150DEST_PATH_IMAGE024
和
Figure 572514DEST_PATH_IMAGE064
以及（胫骨）
Figure 264526DEST_PATH_IMAGE025
和
Figure 459927DEST_PATH_IMAGE065
，采用现有的配对点配准方法分别计算带缩放的三维刚性变换
Figure DEST_PATH_IMAGE066
和
Figure 920996DEST_PATH_IMAGE067
，该方法是一种闭合解形式，可直接计算得到配对点集间的旋转、平移和缩放因子。
需要说明的是，在配准时，以
Figure 85261DEST_PATH_IMAGE024
和
Figure 513837DEST_PATH_IMAGE025
为目标点集，将
Figure 992223DEST_PATH_IMAGE064
和
Figure 370115DEST_PATH_IMAGE065
配准对齐到目标点集上，得到的变换可以用来重采样影像数据
Figure 157811DEST_PATH_IMAGE020
和
Figure 558836DEST_PATH_IMAGE021
。此外，分别计算
Figure 903230DEST_PATH_IMAGE066
和
Figure 322579DEST_PATH_IMAGE067
的逆变换
Figure DEST_PATH_IMAGE068
和
Figure 31909DEST_PATH_IMAGE069
，并对表面点集
Figure 169498DEST_PATH_IMAGE026
和
Figure 989686DEST_PATH_IMAGE027
进行相应的变换。这样得到对齐到参考空间坐标系的三维模板影像。
具体的，
Figure DEST_PATH_IMAGE070
、
Figure 529121DEST_PATH_IMAGE071
，而
Figure DEST_PATH_IMAGE072
、
Figure 409352DEST_PATH_IMAGE073
。
在上述基础上，基于各个X射线图，对三维模板影像通过迭代的方式进行若干次二维/三维相似变换配准，获得初级三维重建影像，具体包括：
a.基于拍摄各X射线图的外参投影（矩阵）参数，从三维模板影像投影生成各个对应投影方向上的DRR图像，并将三维模板影像中目标结构的内外表面点集投影到各个DRR图像所在平面上，获得对应的表面投影点集；
b.计算每个投影方向的DRR图像与对应X射线图的二维相似变换图像配准，得到二维相似变换，用二维相似变换更新对应的表面投影点集的位置，获得各个投影方向的表面投影点集的新位置；
c.基于各个投影方向上表面投影点集的新位置进行反投影，获得目标结构的新的内外表面点集，目标结构的新的内外表面点集与旧的内外表面点集做配对点云配准，得到三维带缩放因子的刚性变换；
d.用该三维带缩放因子的刚性变换去重采样三维模板影像，重采样后的三维模板影像作为下次迭代的输入，
e.经过一次或多次迭代后，获得初级三维重建影像。
以下详细说明上述获得初级三维重建影像的迭代过程。
具体的，基于上述对齐到参考空间坐标系的三维模板影像，进一步做空间和姿势的对齐。首先从三维影像
Figure DEST_PATH_IMAGE074
和
Figure DEST_PATH_IMAGE075
投射产生数字重建射线图像，也就是DRR图像，同时将
Figure DEST_PATH_IMAGE076
和
Figure DEST_PATH_IMAGE077
投射到这些数字重建射线图像上。数字重建射线图像用于对CT影像进行投影来模拟X射线照射图像，本质上是一种遵循Beer-Lambert定律的计算机图形学光线跟踪技术。
具体的，定义第
Figure 96555DEST_PATH_IMAGE063
个方向的X射线图为
Figure DEST_PATH_IMAGE078
，第
Figure 963842DEST_PATH_IMAGE063
个方向的股骨和胫骨数字重建射线图像为
Figure DEST_PATH_IMAGE079
和
Figure DEST_PATH_IMAGE080
，执行如下步骤：
步骤1、计算参考图像
Figure DEST_PATH_IMAGE081
和待配准图像
Figure 233149DEST_PATH_IMAGE078
间的二维/二维带缩放刚性变换配准，由于分为股骨和胫骨结构，所以将
Figure 346599DEST_PATH_IMAGE078
分别乘以股骨和胫骨结构的二值掩模图得到
Figure DEST_PATH_IMAGE082
和
Figure 396463DEST_PATH_IMAGE083
，在第
Figure 807722DEST_PATH_IMAGE063
个方向配准得到变换
Figure DEST_PATH_IMAGE084
和
Figure 72481DEST_PATH_IMAGE085
。
步骤2、利用变换
Figure 278204DEST_PATH_IMAGE084
和
Figure 628413DEST_PATH_IMAGE085
将
Figure 843363DEST_PATH_IMAGE076
和
Figure 24946DEST_PATH_IMAGE077
在该方向的投影点变换到新的位置，如图5所示。对不同方向上都变换到新位置的表面投影点集进行反投影重建，得到两组新的点集
Figure DEST_PATH_IMAGE086
和
Figure 604832DEST_PATH_IMAGE087
。
步骤3、在
Figure 239075DEST_PATH_IMAGE076
和
Figure 742869DEST_PATH_IMAGE086
之间与
Figure 231488DEST_PATH_IMAGE077
和
Figure 795324DEST_PATH_IMAGE087
之间分别进行配对点云配准，得到三维带缩放刚性变换
Figure DEST_PATH_IMAGE088
和
Figure 109674DEST_PATH_IMAGE089
。需要说明的是，这里是以
Figure 479476DEST_PATH_IMAGE076
和
Figure 573334DEST_PATH_IMAGE077
为目标点集，将
Figure 557339DEST_PATH_IMAGE086
和
Figure 166175DEST_PATH_IMAGE087
配准对齐到目标点集上；此外，计算
Figure 11771DEST_PATH_IMAGE088
和
Figure 209403DEST_PATH_IMAGE089
的逆变换
Figure DEST_PATH_IMAGE090
和
Figure 115042DEST_PATH_IMAGE091
。
步骤4、计算复合的三维带缩放刚性变换：
Figure DEST_PATH_IMAGE092
、
Figure 398125DEST_PATH_IMAGE093
、
Figure DEST_PATH_IMAGE094
以及
Figure 47412DEST_PATH_IMAGE095
。
步骤5、对
Figure 99551DEST_PATH_IMAGE020
和
Figure 176091DEST_PATH_IMAGE021
分别用变换
Figure DEST_PATH_IMAGE096
和
Figure 946470DEST_PATH_IMAGE097
进行重采样，利用变换
Figure DEST_PATH_IMAGE098
和
Figure 399448DEST_PATH_IMAGE099
将表面点集
Figure 306093DEST_PATH_IMAGE026
和
Figure 553535DEST_PATH_IMAGE027
变换到参考空间坐标中，得到一组消除了位置和线性几何差异的模板影像：
Figure DEST_PATH_IMAGE100
，
Figure 811209DEST_PATH_IMAGE101
；
Figure DEST_PATH_IMAGE102
，
Figure 802299DEST_PATH_IMAGE103
。
重复步骤1到5，进行多次迭代优化。需要说明的是，每次迭代的输入三维模板影像都是上一次迭代得到的三维模板影像。步骤4的复合变换在第二次迭代时就是基于上一次迭代的变换结果了，以股骨结构为例，第k次迭代得到：
Figure DEST_PATH_IMAGE104
、
Figure 846608DEST_PATH_IMAGE105
。
在上述基础上，所述基于各个X射线图，对初级三维重建影像通过迭代的方式进行若干次包含畸变校正的二维/三维非线性变换配准，获得最终三维重建影像，具体包括：
a. 基于拍摄各X射线图的外参投影（矩阵）参数，从初级三维重建影像投影生成各个对应投影方向上的DRR图像，并将目标结构的内外表面点集投影到各个DRR图像所在平面上，获得对应的表面投影点集；
b. 计算每个投影方向的DRR图像与对应X射线图的二维B样条自由形变图像配准，得到二维B样条变换，用二维B样条变换更新对应的表面投影点集的位置，获得各个表面投影点集的新位置；
c. 基于各个表面投影点集的新位置进行反投影，获得目标结构的新的内外表面点集，目标结构的新的内外表面点集与旧的内外表面点集之间计算薄板样条变换；
d. 基于薄板样条变换，计算三维形变场，对该三维形变场做畸变校正，得到校正后的三维形变场，利用校正后的三维形变场去重采样初级三维重建影像，重采样后的初级三维重建影像作为下次迭代的输入，
e. 经过一次或多次迭代后，获得最终三维重建影像。
以下详细说明上述获得最终三维重建影像的迭代过程。
具体的，上述获得最终三维重建影像的迭代过程与获得初级三维重建影像的迭代过程相似。但有几点不同，在获得最终三维重建影像的迭代过程中：1）在X射线图与数字重建射线图像之间做的是二维/二维非线性图像配准，对应图1所示的子流程“基于表面点投影/反投影的配对点非线性变换配准”；2）由点集
Figure DEST_PATH_IMAGE106
和
Figure 451902DEST_PATH_IMAGE107
反投影重建得到
Figure DEST_PATH_IMAGE108
和
Figure DEST_PATH_IMAGE109
后，计算的是薄板样条变换；3）计算三维形变场，对该三维形变场做畸变校正，也就是图1所示的子流程“基于正则化约束的三维自由形变场畸变校正”。
具体的，三维薄板样条变换的数学形式如公式7所示，由前四项仿射变换与后面的非线性变换项两部分组合而成。
Figure DEST_PATH_IMAGE110
(7)。
这里，空间中任意一点记为
Figure DEST_PATH_IMAGE111
，
Figure DEST_PATH_IMAGE112
为薄板样条控制点，而
Figure DEST_PATH_IMAGE113
。令
Figure DEST_PATH_IMAGE114
，
Figure DEST_PATH_IMAGE115
，
Figure DEST_PATH_IMAGE116
，
Figure 383823DEST_PATH_IMAGE117
，
Figure DEST_PATH_IMAGE118
，
Figure 913025DEST_PATH_IMAGE119
，三维薄板样条变换就由这组参数
Figure DEST_PATH_IMAGE120
，
Figure 731945DEST_PATH_IMAGE121
，
Figure DEST_PATH_IMAGE122
，
Figure 836036DEST_PATH_IMAGE123
，
Figure DEST_PATH_IMAGE124
，
Figure 553457DEST_PATH_IMAGE125
决定，而通过配对点集合
Figure DEST_PATH_IMAGE126
和
Figure 73300DEST_PATH_IMAGE127
可求出这些参数。具体地，把目标点集
Figure 887672DEST_PATH_IMAGE127
带入公式7中求解变换系数，建立如下所示的线性系统：
Figure DEST_PATH_IMAGE128
(8)。
Figure 840628DEST_PATH_IMAGE129
(9)。
Figure DEST_PATH_IMAGE130
是N
Figure 294612DEST_PATH_IMAGE028
4大小的矩阵，任意行为
Figure 431195DEST_PATH_IMAGE131
，
Figure DEST_PATH_IMAGE132
是4
Figure 224707DEST_PATH_IMAGE028
4的零矩阵，
Figure 218071DEST_PATH_IMAGE133
为全0列向量，而
Figure DEST_PATH_IMAGE134
，
Figure 97034DEST_PATH_IMAGE135
，
Figure DEST_PATH_IMAGE136
。公式8所述的线性最小二乘问题可通过诸如QR分解方法计算薄板样条变换的系数。
在本申请中，
Figure 958680DEST_PATH_IMAGE108
和
Figure 482065DEST_PATH_IMAGE109
为源点集，
Figure 98860DEST_PATH_IMAGE106
和
Figure 278169DEST_PATH_IMAGE107
为目标点集。考虑到股骨和胫骨的关节联动性，这里没有独立地计算股骨和胫骨的三维薄板样条变换，而是将源点集和目标点集分别合并得到
Figure 553292DEST_PATH_IMAGE137
和
Figure DEST_PATH_IMAGE138
，计算整体的变换
Figure 321397DEST_PATH_IMAGE139
，可用于从
Figure DEST_PATH_IMAGE140
和
Figure 843514DEST_PATH_IMAGE141
重采样插值重建患者的下肢股骨和胫骨三维形貌体积影像。
进一步，基于正则化约束的三维自由形变场畸变校正的具体方法如下：
实际在使用三维薄板样条变换插值重建时，会存在局部拓扑形状扭曲的情况，主要是由X射线图噪声导致的二维/二维自由形变配准误差引起。所以本申请提出了一种对三维自由形变场进行正则化约束的方法，其步骤如下：
步骤1、在参考空间坐标系中确定影像数据
Figure 510119DEST_PATH_IMAGE140
和
Figure 588933DEST_PATH_IMAGE141
的空间范围，得到包围立方体空间：原点
Figure DEST_PATH_IMAGE142
，长宽高
Figure 205685DEST_PATH_IMAGE143
。定义各向同性的体素尺寸
Figure DEST_PATH_IMAGE144
，如0.5mm，则得到大小为
Figure 164283DEST_PATH_IMAGE145
的离散体素空间
Figure DEST_PATH_IMAGE146
。
步骤2、利用三维薄板样条变换
Figure 52604DEST_PATH_IMAGE139
计算每个体素点
Figure 122060DEST_PATH_IMAGE147
变换后的新位置
Figure DEST_PATH_IMAGE148
，这里，
Figure 412227DEST_PATH_IMAGE149
、
Figure DEST_PATH_IMAGE150
、
Figure 541726DEST_PATH_IMAGE151
为
Figure DEST_PATH_IMAGE152
的体素下标，得到三维自由形变场
Figure 104295DEST_PATH_IMAGE153
，
Figure DEST_PATH_IMAGE154
。
步骤3、从三维自由形变场
Figure 462595DEST_PATH_IMAGE153
建立一个样条控制节点都位于离散体素点上的三维B样条网格，其变换记为
Figure 122115DEST_PATH_IMAGE155
，该变换形式如公式10和公式11，其中公式11为B样条基函数。
Figure DEST_PATH_IMAGE156
(10)。
Figure 173248DEST_PATH_IMAGE157
(11)。
这里的B样条系数
Figure DEST_PATH_IMAGE158
是未知的，需要求解，而在求解前先采用正则化方法对三维自由形变场做自由形变约束/畸变校正。由于样条变换本身具有一定的平滑作用，且直接对整个三维自由形变场做正则约束计算代价较大，所以本申请只对变换
Figure 285429DEST_PATH_IMAGE159
的控制节点（体素下标为
Figure DEST_PATH_IMAGE160
）的位置形变（记为
Figure 447420DEST_PATH_IMAGE161
）做约束，本申请采用现有的正则化方法得到校正的节点位置形变，记为
Figure DEST_PATH_IMAGE162
，计算过程如公式12和13所示。
Figure 633551DEST_PATH_IMAGE163
(12)。
Figure DEST_PATH_IMAGE164
(13)。
这里，将所有节点位置的形变按方向分离构成三组向量
Figure 110712DEST_PATH_IMAGE165
，
Figure DEST_PATH_IMAGE166
和
Figure 460922DEST_PATH_IMAGE167
，并计算它们的离散余弦变换，在公式12的基础上计算离散余弦逆变换对每个节点分别进行校正。
Figure DEST_PATH_IMAGE168
为常量系数，而
Figure 613554DEST_PATH_IMAGE169
表示多维拉普拉斯算子的特征值，
Figure DEST_PATH_IMAGE170
，
Figure 982088DEST_PATH_IMAGE171
为三维B样条网格的尺寸，而
Figure DEST_PATH_IMAGE172
如公式14所示。
Figure 109444DEST_PATH_IMAGE173
(14)。
步骤4、正则化后的节点位置形变量
Figure DEST_PATH_IMAGE174
可帮助计算系数
Figure 133900DEST_PATH_IMAGE175
，本申请采用现有的方法得到变换
Figure 700011DEST_PATH_IMAGE155
的B样条系数，利用得到的
Figure 188630DEST_PATH_IMAGE155
计算
Figure 486887DEST_PATH_IMAGE152
中每个体素的位置形变，再用于对股骨影像
Figure 60957DEST_PATH_IMAGE140
和
Figure 430758DEST_PATH_IMAGE141
进行重采样重建得到患者股骨和胫骨三维形貌体积影像
Figure DEST_PATH_IMAGE176
和
Figure 445988DEST_PATH_IMAGE177
。图6展示了约束前和经过正则化约束后的三维自由形变场及重建结果。
步骤5、根据需要可以重复迭代优化，从第二次迭代开始，可利用现有方法求解
Figure 180725DEST_PATH_IMAGE155
的逆变换
Figure DEST_PATH_IMAGE178
，然后对
Figure 976512DEST_PATH_IMAGE106
和
Figure 87687DEST_PATH_IMAGE107
实施变换得到:
Figure 302898DEST_PATH_IMAGE179
；
Figure DEST_PATH_IMAGE180
。
基于新的模板影像
Figure 395487DEST_PATH_IMAGE176
、
Figure 491619DEST_PATH_IMAGE177
、
Figure 140907DEST_PATH_IMAGE181
和
Figure DEST_PATH_IMAGE182
重复执行图1所示的两步子流程：1）基于表面点投影/反投影的配对点非线性变换配准；2）基于正则化约束的三维自由形变场畸变校正。也就是不断进行迭代，直到达到理想的效果。
通过图7可以看出，通过本申请获得的最终的三维重建影像在投影方向上精确度很高。
本申请的实施例还提供了一种图像处理系统，包括图像获取模块和图像处理模块，所述图像获取模块用于获取包含目标结构的三维模板影像和至少两幅不同方向拍摄的待重建目标结构的X射线图，所述图像处理模块用于按照所述的基于二维X射线图像重建三维影像的方法对三维模板影像进行处理。
本申请的实施例还提供了一种计算机装置，包括：处理器、存储器、通信接口和通信总线，所述处理器、存储器和通信接口通过所述通信总线完成相互间的通信，所述存储器用于存放至少一个可执行指令，所述可执行指令使所述处理器执行所述的基于二维X射线图像重建三维影像的方法对应的操作。
本申请的实施例还提供了一种存储介质，所述计算机存储介质中存储有至少一个可执行指令，所述可执行指令使处理器执行所述的基于二维X射线图像重建三维影像的方法对应的操作。
通过以上分析可知，本申请具有以下有益效果：
1.提出了一种从患者下肢X射线图重建三维形貌体积影像的方法，在不增加辐射剂量的前提下提供了类似CT影像的术前计划输入；
2.提出了一种不依赖统计先验模型的重建方法，统计先验模型需要大量数据才能保证较好的泛化和特异性，而本申请的方案提出的自由形变策略则不受限制；
3.提出了一种新颖的患者下肢三维形貌影像重建方案，巧妙合理地将薄板样条变换、B样条插值平滑、自适应正则化畸变校正串联在一起，并且允许重复迭代优化重建结果；
4.提出了一种同时重建股骨和胫骨结构的策略，考虑了膝关节联动性，在关节处的重建误差会降低，避免比如股骨和胫骨结构相连的错误。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (8)

1.一种基于二维X射线图像重建三维影像的方法，其特征在于，包括：
获取包含目标结构的三维模板影像和至少两幅不同方向拍摄的待重建目标结构的X射线图，将各个X射线图校准对齐到同一参考空间坐标系，
基于各个方向的X射线图和三维模板影像，将三维模板影像初始对齐到参考空间坐标系中，
基于各个方向的X射线图，对三维模板影像通过迭代的方式进行若干次二维/三维相似变换配准，获得初级三维重建影像，
基于各个方向的X射线图，对初级三维重建影像通过迭代的方式进行若干次包含畸变校正的二维/三维非线性变换配准，获得最终三维重建影像；
所述三维模板影像是基于一组CT影像获得；
所述三维模板影像是基于一组CT影像获得，具体包括：
对一组CT影像进行三维影像分割预处理，获得只包含目标结构的三维体素影像数据集
Figure DEST_PATH_IMAGE002
，
对三维体素影像数据集
Figure DEST_PATH_IMAGE004
进行影像配准操作，具体包括：
从三维体素影像数据集
Figure DEST_PATH_IMAGE004A
中任选一个作为参考影像
Figure DEST_PATH_IMAGE006
，利用三维仿射变换配准方法把三维体素影像数据集中的其它各个影像与参考影像
Figure 579891DEST_PATH_IMAGE006
对齐，得到相对于参考影像
Figure DEST_PATH_IMAGE008
消除了位置和线性几何差异的影像数据集
Figure DEST_PATH_IMAGE010
，
基于三维B样条自由形变配准方法，把每个影像数据
Figure DEST_PATH_IMAGE012
配准到参考影像
Figure 811546DEST_PATH_IMAGE006
上，得到三维B样条变换参数集合
Figure DEST_PATH_IMAGE014
和一组自由形变后形貌对准的目标结构的影像数据集合
Figure DEST_PATH_IMAGE016
，
对影像数据集合
Figure 528966DEST_PATH_IMAGE016
与三维B样条变换参数集合
Figure 861859DEST_PATH_IMAGE014
分别求平均，得到均值影像
Figure DEST_PATH_IMAGE018
和三维B样条均值变换参数
Figure DEST_PATH_IMAGE020
，
计算三维B样条均值变换参数
Figure 613914DEST_PATH_IMAGE020
的逆变换
Figure DEST_PATH_IMAGE022
来重采样均值影像
Figure 872595DEST_PATH_IMAGE018
，以产生三维体素影像数据
Figure DEST_PATH_IMAGE024
，
Figure DEST_PATH_IMAGE026
，
对三维体素影像数据
Figure DEST_PATH_IMAGE028
进行处理来生成三维模板影像，具体包括：
从三维体素影像数据
Figure 952677DEST_PATH_IMAGE028
获取目标结构的外轮廓二值掩模，将目标结构的外轮廓二值掩模乘以三维体素影像数据
Figure 151577DEST_PATH_IMAGE028
得到独立的目标结构的影像数据，
从三维体素影像数据
Figure 271323DEST_PATH_IMAGE028
获取目标结构的内部结构二值掩模，基于外轮廓二值掩模和内部结构二值掩模产生三维表面形状网格模型，在三维表面形状网格模型上确定三维解剖结构标记点集，
对三维表面形状网格模型进行稀疏化处理，将稀疏化处理后的目标结构的外表面和内表面的网格点合并，得到目标结构的内外表面点集，
所述三维模板影像包括：目标结构的三维影像数据、目标结构的三维解剖结构标记点集、目标结构的内外表面点集。
2.根据权利要求1所述的基于二维X射线图像重建三维影像的方法，其特征在于，所述基于各个X射线图和三维模板影像，将三维模板影像初始对齐到参考空间坐标系中，具体包括：
在X射线图像上获取目标结构的二维解剖结构标记点集，将二维解剖结构标记点集做反投影，计算得到参考空间下X射线图像上解剖结构标记点集的三维位置，将反投影计算得到的三维解剖结构标记点集和三维模板影像中定义的三维解剖结构标记点集做配对点云刚性配准，获得三维刚性变换，用该三维刚性变换将三维模板影像初始对齐到参考空间坐标系中。
3.根据权利要求2所述的基于二维X射线图像重建三维影像的方法，其特征在于，基于各个X射线图，对三维模板影像通过迭代的方式进行若干次二维/三维相似变换配准，获得初级三维重建影像，具体包括：
a.基于拍摄各X射线图的外参投影参数，从三维模板影像投影生成各个对应投影方向上的DRR图像，并将三维模板影像中目标结构的内外表面点集投影到各个DRR图像所在平面上，获得对应的表面投影点集；
b.计算每个投影方向的DRR图像与对应X射线图的二维相似变换图像配准，得到二维相似变换，用二维相似变换更新对应的表面投影点集的位置，获得各个投影方向的表面投影点集的新位置；
c.基于各个投影方向上表面投影点集的新位置进行反投影，获得目标结构的新的内外表面点集，目标结构的新的内外表面点集与旧的内外表面点集做配对点云配准，得到三维带缩放因子的刚性变换；
d.用该三维带缩放因子的刚性变换去重采样三维模板影像，重采样后的三维模板影像作为下次迭代的输入，
e.经过一次或多次迭代后，获得初级三维重建影像。
4.根据权利要求3所述的基于二维X射线图像重建三维影像的方法，其特征在于，所述基于各个X射线图，对初级三维重建影像通过迭代的方式进行若干次包含畸变校正的二维/三维非线性变换配准，获得最终三维重建影像，具体包括：
a. 基于拍摄各X射线图的外参投影参数，从初级三维重建影像投影生成各个对应投影方向上的DRR图像，并将目标结构的内外表面点集投影到各个DRR图像所在平面上，获得对应的表面投影点集；
b. 计算每个投影方向的DRR图像与对应X射线图的二维B样条自由形变图像配准，得到二维B样条变换，用二维B样条变换更新对应的表面投影点集的位置，获得各个表面投影点集的新位置；
c. 基于各个表面投影点集的新位置进行反投影，获得目标结构的新的内外表面点集，目标结构的新的内外表面点集与旧的内外表面点集之间计算薄板样条变换；
d. 基于薄板样条变换，计算三维形变场，对该三维形变场做畸变校正，得到校正后的三维形变场，利用校正后的三维形变场去重采样初级三维重建影像，重采样后的初级三维重建影像作为下次迭代的输入，
e. 经过一次或多次迭代后，获得最终三维重建影像。
5.根据权利要求1所述的基于二维X射线图像重建三维影像的方法，其特征在于，所述目标结构为股骨和胫骨。
6.一种图像处理系统，其特征在于，包括图像获取模块和图像处理模块，所述图像获取模块用于获取包含目标结构的三维模板影像和至少两幅不同方向拍摄的待重建目标结构的X射线图，
所述图像处理模块用于按照权利要求1至5中任意一项所述的基于二维X射线图像重建三维影像的方法对三维模板影像进行处理。
7.一种计算机装置，其特征在于，包括：处理器、存储器、通信接口和通信总线，所述处理器、存储器和通信接口通过所述通信总线完成相互间的通信，所述存储器用于存放至少一个可执行指令，所述可执行指令使所述处理器执行如权利要求1至5中任意一项所述的基于二维X射线图像重建三维影像的方法对应的操作。
8.一种存储介质，其特征在于，所述存储介质中存储有至少一个可执行指令，所述可执行指令使处理器执行如权利要求1至5中任意一项所述的基于二维X射线图像重建三维影像的方法对应的操作。
CN202210777073.6A 2022-07-04 2022-07-04 一种基于二维x射线图像重建三维影像的方法和系统 Active CN114842154B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210777073.6A CN114842154B (zh) 2022-07-04 2022-07-04 一种基于二维x射线图像重建三维影像的方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210777073.6A CN114842154B (zh) 2022-07-04 2022-07-04 一种基于二维x射线图像重建三维影像的方法和系统

Publications (2)

Publication Number Publication Date
CN114842154A CN114842154A (zh) 2022-08-02
CN114842154B true CN114842154B (zh) 2022-11-15

Family

ID=82574074

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210777073.6A Active CN114842154B (zh) 2022-07-04 2022-07-04 一种基于二维x射线图像重建三维影像的方法和系统

Country Status (1)

Country Link
CN (1) CN114842154B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116128936A (zh) * 2023-02-15 2023-05-16 北京纳通医用机器人科技有限公司 配准方法、装置、设备及存储介质
CN117671221B (zh) * 2024-02-01 2024-05-03 江苏一影医疗设备有限公司 基于膝关节有限角图像的数据修正方法、装置及存储介质

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3175790B1 (en) * 2013-11-04 2021-09-08 Ecential Robotics Method for reconstructing a 3d image from 2d x-ray images
CN108765474A (zh) * 2018-04-17 2018-11-06 天津工业大学 一种针对ct与光学扫描牙齿模型的高效配准方法

Also Published As

Publication number Publication date
CN114842154A (zh) 2022-08-02

Similar Documents

Publication Publication Date Title
US10217217B2 (en) Systems and methods for obtaining 3-D images from X-ray information
CN114842154B (zh) 一种基于二维x射线图像重建三维影像的方法和系统
US20210012492A1 (en) Systems and methods for obtaining 3-d images from x-ray information for deformed elongate bones
US7394946B2 (en) Method for automatically mapping of geometric objects in digital medical images
EP1598778B1 (en) Method for automatically mapping of geometric objects in digital medical images
Karade et al. 3D femur model reconstruction from biplane X-ray images: a novel method based on Laplacian surface deformation
US20210007806A1 (en) A method for obtaining 3-d deformity correction for bones
EP1868157A1 (en) Shape reconstruction using X-ray images
CN107154038B (zh) 一种肋骨可视化的肋骨骨折辅助诊断方法
JP2005521502A (ja) 胸部および腹部の画像モダリティの重ね合わせ
CN111260748B (zh) 基于神经网络的数字合成x射线层析成像方法
Kolditz et al. Comparison of extended field-of-view reconstructions in C-arm flat-detector CT using patient size, shape or attenuation information
CN112562082A (zh) 一种三维人脸重建方法及系统
Tong et al. X-ray2Shape: reconstruction of 3D liver shape from a single 2D projection image
WO2019180746A1 (en) A method for obtaining 3-d deformity correction for bones
Papademetris et al. Articulated rigid registration for serial lower-limb mouse imaging
Balestra et al. Articulated statistical shape model-based 2D-3D reconstruction of a hip joint
Zheng 3D volumetric intensity reconsturction from 2D x-ray images using partial least squares regression
WO2019180747A1 (en) Systems and methods for obtaining patient specific instrument designs
Li et al. Combination of intensity-based image registration with 3D simulation in radiation therapy
Yu et al. 2d-3d regularized deformable b-spline registration: Application to the proximal femur
Aghajani et al. Intensity based image registration by minimizing the complexity of weighted subtraction under illumination changes
CN110473241A (zh) 图像配准方法、存储介质及计算机设备
CN116485850A (zh) 基于深度学习的手术导航影像的实时非刚性配准方法及系统
WO2022229816A1 (en) 3d reconstruction of anatomical images

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
CP03 Change of name, title or address
CP03 Change of name, title or address

Address after: Room 601, Building D, Zhonghe (Suzhou) Science and Technology Innovation Port, No. 588 Xiangrong Road, High Speed Rail New City, Xiangcheng District, Suzhou City, Jiangsu Province, 215000 (6th and 7th floors)

Patentee after: Jiangsu Jicui sukesi Technology Co.,Ltd.

Country or region after: Zhong Guo

Address before: 215000 18th floor, Ziguang building (Qidi building), No. 99, nantiancheng Road, Xiangcheng District, Suzhou City, Jiangsu Province

Patentee before: Jiangsu Jicui sukesi Technology Co.,Ltd.

Country or region before: Zhong Guo