一种由二维X光图像重建三维体图像的方法及其应用与流程

本发明涉及计算机视觉和图像处理技术领域，尤其涉及一种由二维x光片重建三维体图像的方法，可应用于对颅骨形态进行重建，生成三维颅面体图像。

背景技术：

从二维x光图像重建三维结构的技术可以分为基于立体视觉的方法、2d-3d配准方法以及基于统计形状的方法。其中，基于立体视觉的方法通常使用一组从不同视点采集的x光图像，利用shape-from-x技术估计三维几何信息。但是生长数据库中在同一时刻常常只有一张x光图像，并不足以利用立体视觉方法重建三维几何结构。此外，只有图像轮廓或者具有丰富纹理部分的几何结构细节可以得到有效重建。基于图像灰度的2d-3d配准方法可以从一张x光图像中重建三维体图像，将常用的3d-3d图像配准测度调整后可用于2d-3d图像配准问题，例如随机秩相关，稀疏直方图，以及互信息的变分近似等。在基于图像灰度的2d-3d配准过程中，二维以及三维图像的概率密度估计，以及利用数字重建x光图像(drr)技术从三维变形体图像中估计对应的二维投影，都通常对应着较为复杂的计算以及较大的时间代价。此外多数基于灰度的2d-3d配准方法仅处理刚性变换，这并不适合从存档的x光图像中重建对应的三维颅面体图像。基于统计的方法试图最小化三维可变形模型的投影与二维x光图像之间差异。其中可变形模型通常由对一组已知的三维模型进行统计分析得到。考虑到复杂的颅面结构，难以调节统计表面模型使颅面结构细节与目标x光图像一致。而基于灰度的统计模型中包含稠密的的几何形态细节，并可用于三维重建。但基于灰度的统计模型在重建过程中仍然需要费时的drr技术，在每步迭代中从变形的统计灰度模型中估计对应的二维投影。基于部分最小二乘回归技术(partialleastsquaresregression(plsr))可以在三维体图像估计过程中避免drr的计算，但是在学习映射函数时，需要额外的三维代理(surrogate)网格模型，从而会增加系统对于数据的需求，同时加大回归模型的训练代价。

在实际应用中，三维几何形态是度量颅面结构生长的重要依据。在医学临床方面，在锥束ct(cbct)等三维成像技术投入临床使用之前，用于记录颅面生长的数据库仅包含二维的x光图像，并没有三维颅面形态的体图像数据，需要对传统的存档文件进行三维重建得到三维颅面形态。上述现有的三维重建方法存在对数据的需求较高、计算较复杂和时间代价较大等不足。

技术实现要素：

为了克服上述现有技术的不足，本发明提供一种由二维x光片重建三维体图像的方法及其应用，使用锥束ct(cbct)图像作为训练数据训练回归模型，本发明方法可应用于对颅骨形态进行重建，生成三维颅面体图像。

本发明的原理是：本发明基于回归森林(随机森林)从x光片重建三维体图像，首先利用成对的三维体图像与其对应的二维x光片训练回归模型。该回归模型被用于描述稠密的三维体图像与二维投影之间的对应关系。在训练过程中，通过最大化节点分裂中体图像分布的信息增益，以获取最优的节点分裂策略，其中二值的分裂测试函数基于最优选择的x光图像的特征通道(channel)。在线重建过程中，新的x光图像被输入回归森林以预测对应的三维体图像。在高斯分布假设下，可以基于回归森林利用混和高斯模型预测三维体图像，该模型对所有独立决策树中得到的体图像分布进行加权累积。引入迭代修正过程对三维体图像进行调整以最大化体图像相对于输入x光片的似然概率。除去离线过程中生成训练数据之外，该方法在重建中不需要费时的数字重建x光图像(drr)计算以及3d表面代理模型。基于该方法从x光片中重建的三维体图像，可以覆盖颅面结构中的形态细节，而传统的基于立体视觉以及统计表面模型仅能获取稀疏几何结构。

本发明提供的技术方案是：

一种由二维x光图像重建三维体图像的方法，包括训练回归模型阶段和在线重建三维体图像阶段；所述训练回归模型阶段利用成对的三维体图像与其对应的二维x光图像训练得到回归森林模型；所述在线重建三维体图像阶段输入一张二维x光图像，采用训练回归模型阶段学习获得的回归森林模型进行预测，得到对应的三维体图像；包括如下步骤：

一、在训练回归模型阶段，将体图像和对应二维x光图像对(v,c)作为训练数据集，构造回归森林模型；具体包括：

11)针对训练数据集中的体图像v，建立体图像的统计灰度模型；针对训练数据集中与体图像v对应的二维x光图像，通过计算得到相应的图像特征；

12)通过样本自举和最大化节点分裂中的信息增益，学习得到回归森林f中的决策树结构，构造回归森林模型；

二、在线重建阶段，输入一张二维x光图像，采用训练回归模型阶段学习获得的回归森林f进行预测，得到对应的三维体图像；包括如下步骤：

21)从输入的二维x光图像中提取得到形状和纹理特征向量；

22)将步骤21)得到的形状和纹理特征向量投入步骤12)所述回归森林f，经过所述回归森林f中的决策树分枝点中的二值函数判定，所述二维x光图像到达回归森林f的叶子节点；

23)使用叶子节点的中心μi作为当前决策树对重建体图像参数给出的预测；

24)由回归森林f所有决策树对应叶子节点中心的加权组合，得到体图像参数的初始预测；

25)采用混合高斯函数表示由回归森林f所预测体图像的似然概率；通过最大化似然概率定义体图像参数的能量函数；

26)通过将能量函数的一阶偏导设置为0，利用迭代割线方法预测体图像参数，当迭代更新小于预先设置的阈值时迭代终止；

27)由步骤11)所述统计灰度模型公式计算得到对应的三维体图像，实现由二维x光图像重建得到三维体图像。

针对上述由二维x光图像重建三维体图像的方法，进一步地，步骤11)具体通过主元分析方法建立所述体图像的统计灰度模型，所述主元分析方法包括：计算训练集中的体图像与参照体图像之间的稠密对应；计算训练集中的体图像相对于参照体图像的偏移场；对所述偏移场进行主元分析，获取体图像偏移场对应的子空间；所述体图像的统计灰度模型表示为式1：

式1中，v表示偏移场子空间中基元b的系数；nr表示参照体图像的个数；nc表示子空间中占优主元的个数；为训练集中的体图像相对于参照体图像所计算的偏移场的均值；体图像参数v由一个nc×nr维的子空间系数向量表示，如式2：

其中，为体图像v的向量表示，即体图像参数；由此得到与式1对应的通过子空间投影合成的体图像。

针对上述由二维x光图像重建三维体图像的方法，进一步地，步骤11)针对训练数据集 中与体图像v对应的二维x光图像，通过计算得到相应的图像特征；所述计算包括如下步骤：

1121)采用金字塔结构在不同的空间设置中描述所述二维x光图像的形状和纹理信息；

1122)利用hog(histogramoforientedgradients)与lbp(localbinarypatterns)特征描述所述二维x光图像的局部灰度梯度的分布和图像块的纹理；

1123)在金字塔每层均匀划分的网格单元中计算所述二维x光图像的特征，将得到的特征向量记作式3：

式3中，h表示hog与lbp的组合特征，np表示金字塔的层数，表示第i层金字塔中网格个数；c代表x光片的图像特征。

针对上述由二维x光图像重建三维体图像的方法，进一步地，步骤12)所述信息增益采用矩阵的迹计算得到，用式4表示：

式4中，对应父节点vj中的样本个数；为结点对应的协方差矩阵的值；为结点对应的协方差矩阵的迹；对应子节点中的样本个数，节点分裂的二值测试函数定义为其中[·]是0-1指示函数；cq为x光片的第q维特征。

针对上述由二维x光图像重建三维体图像的方法，进一步地，步骤24)所述所有树对应叶子节点中心的加权组合表示为式5：

其中，nt为回归森林中决策树的个数；采用平均移位方法获取叶子节点中的占优模态，叶子节点中的占优模态体图像分布使用高斯函数来描述，其中μi和∑i表示占优模态的均值和方差；

通过式5得到体图像参数的初始预测。

针对上述由二维x光图像重建三维体图像的方法，进一步地，步骤25)所述由回归森林f所预测体图像的似然概率由如下混合高斯函数表示：

式6中，f为回归森林，c为输入x光图像特征，nt为回归森林中决策树的个数,μi和∑i表示占优模态的均值和方差；

体图像参数的能量函数通过最大化似然概率定义如下：

式7中，第一项为随机森林所得到的似然概率的负对数，第二项是正则化项，用于避免体图像参数即体图像偏移场子空间坐标的退化；将正则化系数λ设置为0.01，通过将能量函数的一阶偏导设置为0，得到体图像参数，表示如下：

其中，μj为第j棵决策树所给出的预测；系数κ定义为：

上述由二维x光片重建三维体图像的方法应用于对颅骨形态进行重建，生成三维颅面体图像；具体地，使用锥束ct(cbct)图像作为训练数据训练回归模型，即所述步骤11)训练数据集中的体图像v使用临床采集的锥束ct图像；训练数据集中与体图像v对应的二维x光图像由数字重建x光图像(drr)技术构造。在所述步骤12)之前，对所述训练数据集中的体图像v施加随机的刚性变换的扰动，用于仿真在二维x光图像采集过程中人头姿态的轻微变化。

由二维x光图像重建三维体图像的方法在重建三维颅面体图像的应用中，通过重建得到的体图像中三维骨结构表面与真实锥束ct图像中提取的对应结构的表面之间的距离误差小于0.4mm。所述结构包括前额、上颌骨和下颌骨。

与现有技术相比，本发明的有益效果是：

本发明提供一种由二维x光片重建三维体图像的方法，包括训练回归模型阶段和在线重建三维体图像阶段；所述训练回归模型阶段利用成对的三维体图像与其对应的二维x光片训练得到回归森林模型；所述在线重建三维体图像阶段输入一张二维x光片，采用训练回归模型阶段学习获得的回归森林模型进行预测，得到对应的三维体图像。除去离线过程中生成训练数据之外，本发明方法在重建中不需要费时的数字重建x光图像(drr)计算以及3d表面代理模型。基于该方法从x光片中重建的三维体图像，可以覆盖颅面结构中的形态细节，而传统的基于立体视觉以及统计表面模型仅能获取稀疏几何结构。

利用本发明提供的方法，可以从一张x光片重建三维体图像，其中基于回归森林的形态预测仅仅包括分枝节点中的二值比较，大大降低了回归预测的计算代价，同时利用迭代修正可以有效改善重建体图像精度。重建得到的三维体图像中包含颅面形态细节，可用于口腔临床对于颅面生长规律的分析。

附图说明

图1是本发明实施例中提供的由二维x光片重建三维颅面体图像的方法的流程框图。

图2是本发明实施例在训练阶段作为训练数据集采用的二维x光图像和相应的体图像；

其中，(a)为侧位片(二维x光图像)；(b)为相应的体图像。

图3是本发明实施例在在线重建阶段采用的新的二维x光片图像，通过采用金字塔结构在不同的空间设置描述x光片图像的形状以及纹理信息。

图4是本发明实施例经过在线重建阶段采用新的二维x光片图像经过三维重建得到的相应的体图像；

其中，(a)与(b)为重建得到的体图像在不同视点下的显示。

具体实施方式

下面结合附图，通过实施例进一步描述本发明，但不以任何方式限制本发明的范围。

本发明提供的方法基于输入的二维x光片图像进行三维图像重建，得到对应的三维体图像，该方法利用回归森林描述三维体图像与对应二维投影之间的关系，引入迭代修正过程以调整重建的三维体图像使之与输入二维x光图像一致。该方法可用于从口腔临床采集的生长数据库重建得到三维形态信息，并用于对颅面生长规律的分析。

以下实施例利用口腔临床采集的生长数据库，对颅面体进行三维重建，恢复三维颅面体 图像，从而得到颅面形态信息，并用于对颅面生长规律的分析。

图1是本发明实施例中提供的由二维x光片重建三维颅面体图像的方法的流程框图，包括训练阶段和在线重建阶段；

一、训练阶段：将体图像和对应二维x光图像对(v,c)作为训练数据集，构造对应的回归森林；图2是本发明实施例在训练阶段作为训练数据集采用的二维x光图像和相应的体图像；其中，(a)为侧位片(二维x光图像)；(b)为相应的体图像。

训练阶段包括如下步骤：

11)训练数据集中的体图像v使用临床采集的锥束ct图像，由步骤1)计算得到统计灰度模型；

训练数据集中对应的二维x光图像由数字重建x光图像(drr)技术构造，由步骤2)计算得到对应的图像特征；

12)对训练数据集中的体图像v施加随机的刚性变换的扰动，用于仿真在二维x光片(图像)采集过程中人头姿态的轻微变化；

13)通过样本自举和最大化节点分裂中的信息增益，通过步骤3)的方法学习得到回归森林f中的决策树结构；

二、在线重建阶段，针对一张二维x光片，采用阶段一学习获得的回归森林f进行预测，得到对应的三维体图像；本实施例中，图3是本发明实施例在在线重建阶段采用的新的二维x光片图像。图4是本发明实施例经过在线重建阶段采用新的二维x光片图像经过三维重建得到的相应的体图像；

在线重建阶段，通过步骤4)进行颅面体图像结构预测，得到三维重建体图像；具体包括如下步骤：

21)从输入的二维x光片中提取得到形状与纹理特征向量；

22)将步骤21)得到的形状与纹理特征向量投入回归森林(随机森林)，经过一系列决策树分枝点中的二值函数判定，输入的二维x光片最终到达叶子节点；

23)使用叶子节点的中心μi作为当前决策树对重建体图像参数的预测；

24)由回归森林f所有树对应叶子节点中心的加权组合，得到体图像参数的初始预测；

25)采用混合高斯函数表示由回归森林f所预测体图像的似然概率；通过最大化似然概率定义体图像参数的能量函数；

26)通过将能量函数的一阶偏导设置为0，利用迭代割线方法预测体图像参数，当迭代 更新小于预先设置的阈值则迭代终止；

27)由统计灰度模型公式计算得到对应的三维体图像。

本发明方法中包括如下几个关键步骤：

步骤1：计算体图像的统计灰度模型；

考虑到体图像中通常以百万为量级的体素，本发明利用主元分析方法建立体图像的统计灰度模型，该方法从训练体图像集中随机选取参照体图像，包括计算训练集中的体图像与参照体图像之间的稠密对应；计算训练集中的体图像相对于参照体图像的偏移场；对上述偏移场进行主元分析，获取体图像偏移场对应的子空间。

首先估计训练集中的体图像到参照体图像的稠密偏移场。使用demons配准算法估计训练集中的体图像与参照体图像vr之间的非刚性配准，其中稠密的偏移场记做δvr。为了避免参照体图像的选择所带来的对偏移场的有偏估计，系统使用一组随机选择的参照体图像，对应的可变形的体图像表示为式1所示的体图像的统计灰度模型：

式1中，v表示偏移场子空间中基元b的系数；nr表示参照体图像的个数；nc表示子空间中占优主元的个数；为训练集中的体图像相对于参照体图像所计算的偏移场的均值；体图像v由一个nc×nr维的子空间系数向量表示为：

由此得到与公式1对应通过子空间投影合成的体图像。

步骤2：计算得到x光片的图像特征；

系统采用金字塔结构在不同的空间设置中描述x光片图像的形状以及纹理信息，其中利用hog(histogramoforientedgradients，方向梯度直方图)与lbp(localbinarypatterns，局部二值模式)特征描述x光片的局部灰度梯度的分布以及图像块的纹理。在金字塔每层均匀划分的网格单元中计算x光片的图像特征，将得到的特征向量记作式3：

式3中，h表示hog与lbp的组合特征，np表示金字塔的层数，表示第i层金字塔中网格个数；c代表x光片的图像特征。

图3是本发明实施例在在线重建阶段采用的新的二维x光片图像，如图3所示，本发明通过采用金字塔结构在不同的空间设置描述x光片图像的形状以及纹理信息，并计算得到x光片的图像特征。

步骤3：训练得到回归森林；

回归森林具有高效的在线测试以及泛化能力。系统利用回归森林从二维x光片估计对应的三维体图像。给定体图像以及对应二维x光片图像对(v,c)作为训练数据，构造对应的回归森林。训练集中的体图像v使用临床采集的锥束ct图像，并由步骤1)计算得到统计灰度模型；对应的二维x光图像则由数字重建x光图像(drr)技术构造，并由步骤2)计算得到对应的图像特征。对训练集中的体图像施加随机的刚性变换的扰动，用于仿真在二维x光片采集过程中人头姿态的轻微变化。通过样本自举以及最大化节点分裂中的信息增益学习回归森林中的决策树结构。假设节点中样本满足高斯分布，分裂节点所获取的信息增益由子节点中样本的协方差矩阵σ的值决定，考虑到高维特征空间中协方差矩阵的低秩问题，采用矩阵的迹(trace)计算如下信息增益：

式4中，ij为分裂结点vj的信息增益；对应父节点vj中的样本个数；为结点对应的协方差矩阵的值；为结点对应的协方差矩阵的迹；对应子节点中的样本个数，节点分裂的二值测试函数定义为其中cq为x光片的第q维特征；τ为待优化确定的阈值；[·]是0-1指示函数。考虑到训练数据中对于三维体图像的刚性扰动，通过最大化信息增益所获取的最优分裂策略所对应的特征通道更倾向于对应姿态无关的二维x光片图像特征。

步骤4：颅面体图像结构预测，得到三维重建体图像；

本实施例中，在线重建阶段对颅面体图像结构进行预测，具体地，给定输入的二维x光片，通过训练好的回归森林预测对应的三维体图像。回归森林由一组决策树构成，其中叶子结点中的样本可以看作是对原始数据集的聚类。由每个叶子节点中的样本可以学习该叶子结点对应聚类的数据分布。例如在数据高斯分布假设下，可以利用高斯分布函数统计该类数据的均值与方差并得到对应概率分布函数。

从输入的x光片中提取的形状与纹理特征向量被投入随机森林，经过一系列决策树分枝 点中二值函数判定，输入样本最终到达叶子节点。

直观地可以使用叶子节点的中心μi作为当前决策树对重建体数据参数所给出的预测。回归森林所给出的体图像参数的初始预测可由所有树对应叶子节点中心的加权组合得到。

其中，由所有树对应叶子节点中心的加权组合得到的体图像参数的初始预测；nt为回归森林中决策树的个数。这里采用平均移位方法(meanshift)获取叶子节点中的占优模态。叶子节点中的占优模态体图像分布使用高斯函数来描述，其中μi和∑i表示占优模态的均值和方差。由回归森林f所预测体图像的似然概率由如下混合高斯函数表示：

式6中，为随机森林所得到的似然概率函数；f为回归森林；c为输入x光图像特征；nt为回归森林中决策树的个数；μi和∑i表示占优模态的均值和方差；为使用高斯函数描述的占优模态体图像分布。

体图像参数的能量函数通过最大化似然概率定义如下：

其中，为体图像参数的能量函数；第一项为随机森林所得到的似然概率的负对数；f为回归森林；c为输入x光图像特征；第二项是正则化项，用于避免体图像参数即体图像偏移场子空间坐标的退化。在实验中常系数(正则化系数)λ设置为0.01。通过将能量函数的一阶偏导设置为0，可以得到体图像参数的表示如下：

其中，为待求解的体图像参数；nt为回归森林中决策树的个数；μj为第j棵决策树所给出的预测；系数κ定义为：

其中常系数(正则化系数)λ设置为0.01；μ和∑表示占优模态的均值和方差；为使用高斯函数描述的占优模态体图像分布。由于系数κ依赖于体图像参数如上的体图像参数没有解析解。如果将正则化系数λ设置为0，并假设不同决策树对体图像的分布预测一致，则式5中的初始预测可以看作对式7能量函数最优解(式8)的特殊形式。本方法中利用迭代割线方法预测体图像参数，当迭代更新小于预先设置的阈值则迭代终止。在系统中大约需要20步迭代就可以收敛。给定由回归森林估计的体图像参数(体图像偏移场子空间坐标向量)对应的三维体图像可由统计灰度模型(式1)计算得到，由此获得重建的三维颅面体图像。

为了验证上述基于回归森林重建三维颅面体图像方法的精度，本实施例通过计算重建图像中三维骨结构表面与真实cbct图像中提取的对应结构的表面之间的距离，实验测试了前额、上颌骨与下颌骨的外表面，结果表明，距离误差小于0.4mm，可满足临床口腔从颅面生长数据库中分析颅面生长的精度要求。

利用本发明的方法，可由一张x光片中重建三维颅面体图像，该方法基于统计灰度模型以及回归森林构造三维体图像与二维x光片之间的关系。在数据的高斯分布假设下，回归森林可以得到混合高斯模型形式的体图像分布的似然概率。系统结合正则化约束以避免统计灰度模型中的子空间系数的退化。系统中三维体图像的重建以迭代方式进行求解，用于保证重建体图像投影与输入x光片一致。该系统有效克服了传统方法中需要费时的drr图像投影计算以及三维表面代理模型估计映射函数等问题，同时重建的体图像可以覆盖颅面结构中的形态细节。

需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。