CN108205805B - 锥束ct图像间体素稠密对应的自动建立方法 - Google Patents

锥束ct图像间体素稠密对应的自动建立方法 Download PDF

Info

Publication number
CN108205805B
CN108205805B CN201611184654.XA CN201611184654A CN108205805B CN 108205805 B CN108205805 B CN 108205805B CN 201611184654 A CN201611184654 A CN 201611184654A CN 108205805 B CN108205805 B CN 108205805B
Authority
CN
China
Prior art keywords
voxel
voxels
image
similarity
random forest
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
CN201611184654.XA
Other languages
English (en)
Other versions
CN108205805A (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.)
Peking University
Original Assignee
Peking University
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 Peking University filed Critical Peking University
Priority to CN201611184654.XA priority Critical patent/CN108205805B/zh
Publication of CN108205805A publication Critical patent/CN108205805A/zh
Application granted granted Critical
Publication of CN108205805B publication Critical patent/CN108205805B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/18—Image warping, e.g. rearranging pixels individually
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • 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/20—Special algorithmic details
    • G06T2207/20076—Probabilistic image processing
    • 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/20—Special algorithmic details
    • G06T2207/20081—Training; Learning

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Multimedia (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)

Abstract

本发明公布了一种锥束CT图像间体素稠密对应的自动建立方法，分为对级联测地随机森林进行训练过程和在线测试过程；首先根据锥束CT图像定义图像子集，基于表观随机森林获取图像子集间体素的相似度，然后基于级联测地随机森林更新图像子集间体素的相似度，再利用正则化机制获取原始图像间稠密体素对应，由此实现锥束CT图像间体素稠密对应的快速自动建立。本发明建立对应的方法有效克服了现有方法对标注数据的大量需求以及稠密体素对应效率低的问题，基于本发明的体素对应可得到锥束CT图像之间的非刚性变形以及配准，用于估计不同锥束CT图像之间的差异及临床正畸治疗评价。

Description

锥束CT图像间体素稠密对应的自动建立方法
技术领域
本发明涉及口腔临床医学和计算机视觉领域，具体涉及一种自动建立锥束CT图像之间体素稠密对应关系(相似度)的方法。
背景技术
在正畸学临床中使用锥束CT图像记录并度量治疗前后颅面结构的变化，其中图像之间体素的稠密对应以及图像配准是量化图像结构变化的关键。考虑到锥束CT图像中较低的信噪比、图像中颅面结构中由于治疗与生长造成的形态差异、数据采集中微小的姿态变化，以及三维锥束CT极大的数据规模，在锥束CT图像之间建立体素的稠密对应是一项有挑战性的任务。
现有技术中，基于常见的互信息、归一相关、以及拉普拉斯嵌入流形距离等测度来计算图像的非刚性配准以及体素的稠密对应，一般费时并易于出现局部最小。最近出现的基于缩减子集的方法用于对锥束CT图像进行刚性的重叠，但是，由于一般临床正畸治疗会持续较长的时间并包含非刚性的形态变化，因而刚性的重叠方法并不适用于对锥束CT图像建立体素稠密对应。随机森林技术可用于体素分类以及对应，但是，采用有监督的分类随机森林算法则需要大量的标注数据或者基于超体素分解获取的标签。基于超体素分解的方法可以自动获取对应体素的标签，但是，基于超体素分解的方法仅仅可以使用来自一个锥束CT图像的体素训练随机森林，难以从有限的样本中获取泛化的分类以及对应关系(相似度)。
发明内容
为了克服上述现有技术的不足，本发明提供一种自动建立锥束CT图像之间体素稠密对应的方法，基于非监督级联随机森林技术快速建立锥束CT图像之间体素的稠密对应，进而可以得到锥束CT图像之间的非刚性变形以及配准，可用于估计不同锥束CT图像之间的差异及临床正畸治疗评价。
本发明的原理是：针对获取锥束CT图像之间体素的稠密对应问题，本发明提出基于非监督级联随机森林技术建立体素的稠密对应，其中以三维SIFT(Scale-invariantfeature transform，尺度不变特征变换)技术抽取卷积平滑差异图像的局部极值构造图像子集，并在该图像子集上建立表观随机森林。以非监督聚类随机森林算法生成表观随机森林，获取体素之间的相似度，进而获取图像子集之间体素映射函数。除了常见的基于上下文体素灰度的特征描述之外，本发明利用测地坐标以克服近邻结构的对应混淆。在图像子集中基于体素相似度矩阵构造加权的无向图，其中在空间近邻的体素之间建立边连接。体素的测地坐标定义为该无向图中的最短路径长度，并在体素的测地坐标上建立测地随机森林。由于测地坐标可以描述体素关于结构与背景之间的差异，由测地随机森林所得到的体素相似度反应了体素是否位于相同解剖结构，并用于修正基于表观随机森林得到的体素相似度。由于图像子集中的加权无向图会随着体素的相似度矩阵更新而发生变化，因而可以在迭代更新的测地坐标上构造级联测地随机森林。级联测地随机森林中的每一层都可用于修正当前图像子集体素的相似度。利用多数投票机制建立图像子集间的体素相似度。最后利用正则化机制将图像子集中的体素对应扩散到整个锥束CT图像，获取锥束CT图像之间体素的稠密对应。
本发明提供的技术方案是：
一种锥束CT图像间体素稠密对应的自动建立方法，分为对级联测地随机森林进行训练过程和在线测试过程；首先根据锥束CT图像定义图像子集，基于表观随机森林获取图像子集间体素的相似度，然后基于级联测地随机森林更新图像子集间体素的相似度，再利用正则化机制获取原始图像间稠密体素对应，由此实现锥束CT图像间体素稠密对应的自动建立；具体包括如下步骤：
(一)对级联测地随机森林进行训练的过程，包括：
1)基于锥束CT图像子集的表观信息，以非监督聚类随机森林算法训练表观随机森林，从表观随机森林获取图像子集间体素的相似度，构建体素相似度矩阵A；连接相邻体素得到无向图，根据体素相似度矩阵A设置所得到的无向图中的边的权值；包括：
11)随机从图像子集中选取体素构造表观随机森林中的决策树，再预测体素之间的相似度；
12)通过式1计算得到每棵决策树生成体素对(vk,vl)为相似体素的概率：
pa(φ(vk)＝vl)＝Pkl/νkl (式1)
其中，Pkl为体素vk与vl从根节点到叶节点共有遍历路径的长度；νkl＝max(Pk,Pl)为体素vk与vl从根节点到叶节点的遍历长度Pk与Pl的最大值；φ为待求解的图像子集之间的映射函数；
13)定义体素相似度矩阵A，根据体素相似度矩阵A，设置连接相邻体素所得到的无向图中的边的权值；所述相似度矩阵A中的元素akl定义为体素对(vk,vl)从nt棵独立的决策树所生成的相似概率的均值，
Figure BDA0001186108570000031
为图像子集间体素的相似度；所述无向图中连接相邻体素的边权值由相似度矩阵A中与该相邻体素对应的元素确定；
2)基于级联测地森林更新图像子集间体素的相似度：
21)定义测地坐标g为从边界背景到当前体素在图像子集的无向图中的最短路径长度；根据步骤13)得到的无向图中连接相邻体素的边的权值，估计得到图像子集体素的测地坐标，并在图像子集体素的测地坐标上以非监督聚类森林算法建立测地随机森林；
22)由所述测地随机森林获取体素相似度，更新步骤13)从表观随机森林中获取的体素相似度矩阵A；从测地随机森林中返回的体素相似度反映体素之间在测地距离意义上的相似；更新后的体素相似度
Figure BDA0001186108570000032
定义为式2：
Figure BDA0001186108570000033
其中，
Figure BDA0001186108570000034
是从测地随机森林各个决策树中得到的体素属于相同结构概率的平均值；akl表示从表观随机森林中获取的体素相似度；φ表示在图像子集体素v之间的映射函数；
23)根据相似度矩阵更新测地坐标，依据更新后的测地坐标建立一个新的测地随机森林；对级联测地随机森林进行训练时，进行nk次的测地坐标更新得到nk层的级联测地随机森林，在第i次迭代中，相似度
Figure BDA0001186108570000035
其中
Figure BDA0001186108570000036
分别对应第i步迭代中体素对(vk,vl)的相似度，第i-1步迭代中体素对(vk,vl)的相似度，以及第i-1步迭代中测地随机森林得到的体素相似度；
(二)在线测试过程，包括：
1)从表观随机森林中获取图像子集Sr中的体素与图像子集St中的体素为对应体素的概率；
2)根据级联测地随机森林返回的体素属于同一个结构的概率，对表观随机森林获取的体素相似度进行多次迭代更新；在第i次迭代中，由当前的相似度矩阵确定无向图中邻接边的权重进而估计测地坐标，将该测地坐标输入到对应的测地随机森林得到体素位于相同结构的概率，并以此更新图像子集之间的体素相似度；每次迭代更新的图像子集的相似度都用于估计图像子集之间体素的映射函数
Figure BDA0001186108570000037
多数投票机制被用于从多次迭代得到的映射之中计算得到最后的映射函数φ；
3)根据图像子集的映射函数φ，利用正则化机制获取原始图像间稠密体素对应；
假定参照图像与目标图像之间的位移场记做C，其元素ci对应从参照图像Vr中的体素
Figure BDA0001186108570000041
到该体素在目标图像Vt中的对应体素之间的位移；根据图像子集的映射函数φ，预先确定一组体素对应，其中
Figure BDA0001186108570000042
通过式3估计整个图像对应的稠密位移场：
Figure BDA0001186108570000043
式3中，E为能量函数；C为图像对应的稠密位移场；μ1,μ2为常数系数；式3包括三项：
第一项是施加了位移场的变形后的参照图像与目标图像之间的图像灰度差异，其中
Figure BDA0001186108570000044
分别对应施加了位移场ci的参照图像的体素
Figure BDA0001186108570000045
的灰度与在目标图像中的对应体素
Figure BDA0001186108570000046
的灰度；Φ表示待求解的原始锥束CT图像Vr与Vt之间的稠密体素映射函数；当Vr中的体素在Vt中没有对应时，系数γi设置为0；否则，γi设置为1；
第二项中，
Figure BDA0001186108570000047
是由表观森林与级联测地森林获取的图像子集对应后，计算得到体素
Figure BDA0001186108570000048
与其在目标图像子集中的对应体素之间的位移；ck表示在待求解稠密位移场中体素
Figure BDA0001186108570000049
的位移；通过第二项约束待求解的位移场与之前由表观森林与级联测地森林获取的图像子集对应一致；
第三项是平滑项，通过约束相邻的体素之间保持相似的位移以获取平滑的位移场，其中
Figure BDA00011861085700000410
是位移场的梯度；
求解式3，得到最终原始锥束CT图像Vr与Vt之间的稠密体素映射函数Φ为：
Figure BDA00011861085700000411
其中，
Figure BDA00011861085700000412
分别为参照图像与目标图像Vt中的体素；ci为体素
Figure BDA00011861085700000413
的位移；
由此实现快速获取原始锥束CT图像之间体素的稠密对应。
针对上述锥束CT图像间体素稠密对应的自动建立方法，进一步地，所述根据锥束CT图像定义图像子集，具体采用三维尺度不变特征变换SIFT方法抽取卷积平滑差异图像的局部极值构造图像子集。
针对上述锥束CT图像间体素稠密对应的自动建立方法，进一步地，具体采用非监督的表观随机森林构造体素对应。
针对上述锥束CT图像间体素稠密对应的自动建立方法，进一步地，所述训练过程中，步骤11)随机从图像子集中选取体素构造表观随机森林中的决策树具体通过最大化信息增益得到节点的最优分裂策略，并基于该策略将节点的体素分到左子节点以及右子节点；在非监督的表观随机森林训练中，定义节点分裂中的信息增益为节点分裂后的左右子节点中体素的协方差矩阵的迹和的负数；由此建立决策树；根据建立的决策树估计并预测体素之间的相似度，具体是：将体素投入决策树的根节点后，依据节点中保存的最优分裂策略，该体素最终落入叶子节点中；根据两个体素共有遍历路径的长度定义该体素对的相似度(如式1)。
针对上述锥束CT图像间体素稠密对应的自动建立方法，进一步地，所述训练过程中，步骤21)建立测地随机森林时，在规则采样的锥束CT图像中，背景节点B定义在锥束CT图像的六个边界平面上。
针对上述锥束CT图像间体素稠密对应的自动建立方法，进一步地，所述在线测试过程中，步骤2)具体利用多数投票机制从多次迭代得到的映射中计算最后的映射函数φ。
针对上述锥束CT图像间体素稠密对应的自动建立方法，进一步地，所述在线测试过程中，步骤3)利用Levenberg-Marquardt算法求解式3能量最小化问题。
与现有技术相比，本发明的有益效果是：
本发明提供一种建立锥束CT图像之间体素稠密对应的方法，基于非监督随机森林技术快速建立锥束CT图像之间体素的稠密对应。利用本发明提供的方法，可以快速建立锥束CT图像之间体素的稠密对应。基于该体素对应可以得到锥束CT图像之间的非刚性变形以及配准，用于估计不同锥束CT图像之间的差异及临床正畸治疗评价。
本发明方法基于非监督的表观随机森林，在训练阶段不需要预先标注的大规模数据集；利用测地坐标并构造测地随机森林，迭代更新体素属于相同结构的概率，并以此修正基于体素上下文表观信息获取的体素相似度。由于使用缩减的图像子集训练表观随机森林，表观随机森林的训练集中可以包含大量的锥束CT图像数据。本发明方案在图像子集之间可以更快地获取对应，并据此利用正则化机制得到原始大规模锥束CT图像之间体素的稠密对应。因此，本发明方法有效克服了传统对应算法对标注数据的大量需求以及稠密体素对应效率低的问题。
附图说明
图1是本发明方法的流程框图。
图2本发明实施例采用图像子集建立体素对应的实施流程示意图。
具体实施方式
下面结合附图，通过实施例进一步描述本发明，但不以任何方式限制本发明的范围。
本发明提供一种建立锥束CT图像之间体素稠密对应的方法，基于非监督随机森林技术快速建立锥束CT图像之间体素的稠密对应。利用本发明提供的方法，可以快速建立锥束CT图像之间体素的稠密对应。基于该体素对应可以得到锥束CT图像之间的非刚性变形以及配准，用于估计不同锥束CT图像之间的差异及临床正畸治疗评价。图1所示是本发明方法的流程，包括训练阶段和在线测试阶段，首先根据锥束CT图像定义图像子集，再建立锥束CT图像之间体素的稠密对应；体素稠密对应的自动建立主要包括如下步骤：
1)基于表观随机森林获取图像子集间体素的相似度；
本发明利用非监督的表观随机森林构造体素对应，以避免费力的锥束CT图像数据标注。
2)基于级联测地随机森林更新图像子集间体素的相似度；
3)利用正则化机制获取原始图像间稠密体素对应。
图2是本发明实施例采用图像子集建立体素对应的实施流程示意图。采用本发明方法构造两个锥束CT图像(参照图像与目标图像)之间体素的稠密对应。实施例中的训练集包含100个锥束CT图像，利用三维SIFT技术从锥束CT图像中提取大约5万特征点，并将图像中提取的特征点定义为图像子集。在所有图像的特征点上构造表观随机森林与级联测地随机森林。每个随机森林由10棵左右的决策树组成。
1)基于表观随机森林获取图像子集间体素的相似度
随机森林具有高效的在线测试以及泛化能力，广泛用于医学图像处理中的结构检测、分割以及对应。通常随机森林以有监督的方式进行训练，并需要一个预先标注的图像集，或者基于超体素分割获取的受限体素标签。为了避免费力的锥束CT图像数据标注，本发明利用非监督的随机森林技术构造体素对应。
11)随机从图像子集中选取体素构造表观随机森林中的决策树；估计并预测体素之间的相似度；
其中，通过最大化信息增益得到节点的最优分裂策略，并基于该策略将节点的体素分到左子节点以及右子节点。在非监督的表观随机森林训练中，节点分裂中的信息增益定义为节点分裂后的左右子节点中体素的协方差矩阵的迹和的负数。在建立决策树后，将体素投入决策树的根节点后，依据节点中保存的最优分裂策略，该体素会最终落入叶子节点中。根据两个体素共有遍历路径的长度定义该体素对的相似度(式1)，不同图像之间的相似体素即对应体素。
12)每棵决策树生成体素对(vk,vl)为对应体素(相似体素)的概率；
每棵决策树可生成体素对(vk,vl)为相似体素的概率，表示为式1：
pa(φ(vk)＝vl)＝Pkl/νkl (式1)
其中，Pkl为体素vk与vl从根节点到叶节点共有遍历路径的长度；νkl＝max(Pk,Pl)为体素vk与vl从根节点到叶节点的遍历长度Pk与Pl的最大值；φ为待求解的图像子集之间的映射函数。
13)定义体素相似度矩阵A，根据体素相似度矩阵A，设置连接相邻体素所得到的无向图中的边的权值；
相似度矩阵A中的元素akl定义为体素对(vk,vl)从nt棵独立的决策树所生成的相似概率的均值，即
Figure BDA0001186108570000071
其中φ:Sr→St，定义为图像子集Sr与St之间的映射函数。在图像子集中通过连接近邻体素(例如k近邻)定义无向图结构，其中连接相邻体素的边权值由相似度矩阵A中的对应该相邻体素之间的相似度确定。
2)基于级联测地森林更新图像子集间体素的相似度
21)估计图像子集体素的测地坐标，并依据图像子集的测地坐标构造测地随机森林；
由于空间近邻的体素可能属于不同的结构，使用欧式坐标作为体素的空间特征有明显的缺陷，本系统引入测地坐标估计体素的结构先验。测地坐标g定义为从边界背景到当前体素在图像子集的无向图中的最短路径长度，即g(v)＝mind(v,B)，其中d(v,B)为体素v与边界节点B之间候选路径距离。具有相同测地坐标的体素被假设位于相同的解剖结构上。在规则采样的锥束CT图像中，背景节点B定义在锥束CT图像的六个边界平面上。
一旦从表观随机森林得到图像子集体素的相似度矩阵，就可以定义图像子集的无向图中边连接的权值，进而估计体素的测地坐标。在图像子集体素的测地坐标上建立测地随机森林。测地随机森林返回的体素相似度反映体素之间在测地距离意义上的相似。例如测地随机森林返回的相似体素更有可能位于相同结构上。
22)由测地随机森林获取体素相似度，更新步骤13)从表观随机森林中获取的体素相似度矩阵A；
测地随机森林获取体素相似度反映体素之间在测地距离意义上的相似。由测地随机森林获取的体素相似度用于修正从表观随机森林中获取的体素相似度矩阵A。该更新被认为是对相同解剖结构上的体素相似度进行增强，反之抑制在不同解剖结构上体素的相似度。更新后的体素相似度
Figure BDA0001186108570000081
定义为式2：
Figure BDA0001186108570000082
其中，
Figure BDA0001186108570000083
是从测地随机森林各个决策树中得到的体素v属于相同结构概率的平均值；akl表示从表观随机森林中获取的体素相似度；φ表示在图像子集体素之间的映射函数。
23)对级联测地随机森林进行训练时，更新测地坐标，依据更新后的测地坐标建立一个新的测地随机森林，进行nk次的测地坐标更新得到nk层的级联测地随机森林；当进行在线测试时，用训练好的级联测地随机森林返回的体素属于同一个结构的概率，对表观随机森林获取的体素相似度进行多次更新；
一旦相似度矩阵发生变化，图像子集中的无向图的权值也随之改变，进而从无向图中估计的测地坐标也进行更新，依据更新后的测地坐标建立一个新的测地随机森林。在级联测地随机森林的训练过程中，进行nk次的测地坐标更新得到nk层的级联测地随机森林，在第i次迭代中，相似度
Figure BDA0001186108570000084
其中
Figure BDA0001186108570000085
分别对应第i步迭代中体素对(vk,vl)的相似度，第i-1步迭代中体素对(vk,vl)的相似度，以及第i-1步迭代中测地随机森林得到的体素相似度。
在线测试过程中，从表观森林中可以获取图像子集Sr中的体素与图像子集St中的体素为对应体素的概率，随后，级联测地随机森林返回的体素属于同一个结构的概率对表观随机森林获取的体素相似度进行多次更新。在第i次迭代中，由当前的相似度矩阵确定无向图中邻接边的权重进而估计测地坐标，将该测地坐标输入到对应的测地随机森林得到体素位于相同结构的概率，并以此更新图像子集之间的体素相似度。每次迭代更新的图像子集的相似度都可以用于估计图像子集之间体素的映射函数
Figure BDA0001186108570000086
其中对应体素
Figure BDA0001186108570000087
与
Figure BDA0001186108570000088
之间的相似度
Figure BDA0001186108570000089
最大。多数投票机制被用于从多次迭代得到的映射之中计算最后的映射函数φ。
3)利用正则化机制获取原始图像间稠密体素对应
从表观随机森林以及级联测地随机森林可以获取图像子集之间的体素对应，正则化技术被用于获取原始图像中的稠密体素对应。假定参照图像与目标图像之间的位移场记做C，其元素ci对应从参照图像Vr中的体素
Figure BDA0001186108570000091
到该体素在目标图像Vt中的对应体素之间的位移。当给定图像子集的映射函数φ,可以预先知道参照图像子集Sr与目标图像子集St中一组体素对应，其中
Figure BDA0001186108570000092
估计整个图像对应的稠密位移场则被定义为如式3的正则化问题，其中能量函数E定义为：
Figure BDA0001186108570000093
式3中，C为图像对应的稠密位移场；μ1,μ2为常数系数；
第一项是施加了位移场的变形后的参照图像与目标图像之间的图像灰度差异，其中
Figure BDA0001186108570000094
分别对应施加了位移场ci的参照图像的体素
Figure BDA0001186108570000095
的灰度与在目标图像中的对应体素
Figure BDA0001186108570000096
的灰度；Φ表示待求解的原始锥束CT图像Vr与Vt之间的稠密体素映射函数；由于Vr与Vt之间的形态变化，当Vr中的体素在Vt中没有对应时，此时系数γi设置为0；否则γi设置为1；
第二项对应从表观随机森林与级联测地随机森林获取的图像子集中的对应体素之间的位移
Figure BDA0001186108570000097
与原始图像中的对应体素的位移要求一致；第二项中，
Figure BDA0001186108570000098
是由表观森林与级联测地森林获取的图像子集对应后，计算得到体素
Figure BDA0001186108570000099
与其在目标图像子集中的对应体素之间的位移，ck表示在待求解稠密位移场中体素
Figure BDA00011861085700000910
的位移；通过第二项约束待求解的位移场与之前由表观森林与级联测地森林获取的图像子集对应一致；
第三项是平滑项，通过约束相邻的体素之间保持相似的位移以获取平滑的位移场，其中
Figure BDA00011861085700000911
是位移场的梯度。
利用Levenberg-Marquardt(LM)算法求解该能量(式3)最小化问题。最终原始锥束CT图像Vr与Vt之间的稠密体素映射函数定义为：
Figure BDA00011861085700000912
其中
Figure BDA00011861085700000913
分别为参照图像与目标图像Vt中的体素，ci为体素
Figure BDA00011861085700000914
的位移。
利用本发明的方法，可以快速建立锥束CT图像之间体素的稠密对应，其中表观随机森林与级联的测地随机森林返回图像子集之间的映射函数，在该映射的指导下，通过正则化机制快速获取原始锥束CT图像之间体素的稠密对应。该算法基于非监督的随机森林技术，在训练阶段不需要预先标注的大规模数据集。该算法利用测地坐标并构造测地随机森林，迭代更新体素属于相同结构的概率，并以此修正基于体素上下文表观信息获取的体素相似度。由于使用缩减的图像子集训练表观随机森林与级联测地随机森林，训练集中可以包含大量的锥束CT图像数据。在图像子集之间可以更快地获取对应，并据此利用正则化机制得到原始大规模锥束CT图像之间体素的稠密对应。该方法有效克服了传统对应算法对标注数据的大量需求以及稠密体素对应效率低的问题。
需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。

Claims (6)

1.一种锥束CT图像间体素稠密对应的自动建立方法，分为对级联测地随机森林进行训练过程和在线测试过程；首先根据锥束CT图像定义图像子集，基于表观随机森林获取图像子集间体素的相似度，然后基于级联测地随机森林更新图像子集间体素的相似度，再利用正则化机制获取原始图像间稠密体素对应，由此实现锥束CT图像间体素稠密对应的快速自动建立；具体包括如下步骤：
(一)对级联测地随机森林进行训练的过程，包括：
1)基于锥束CT图像子集的表观信息，以非监督聚类随机森林算法训练表观随机森林，从表观随机森林获取图像子集间体素的相似度，构建体素相似度矩阵A；连接相邻体素得到无向图，根据体素相似度矩阵A设置无向图中的边的权值；包括：
11)随机从图像子集中选取体素构造表观随机森林中的决策树，再预测体素之间的相似度；
12)通过式1计算得到每棵决策树生成体素对(vk,vl)为相似体素的概率：
pa(φ(vk)＝vl)＝Pkl/νkl (式1)
其中，Pkl为体素vk与vl从根节点到叶节点共有遍历路径的长度；νkl＝max(Pk,Pl)为体素vk与vl从根节点到叶节点的遍历长度Pk与Pl的最大值；φ为待求解的图像子集之间的映射函数；
13)定义体素相似度矩阵A，根据体素相似度矩阵A，设置连接相邻体素所得到的无向图中的边的权值；所述相似度矩阵A中的元素akl定义为体素对(vk,vl)从nt棵独立的决策树所生成的相似概率的均值；所述无向图中连接相邻体素的边的权值由相似度矩阵A中与该相邻体素对应的元素确定；
2)基于级联测地森林更新图像子集间体素的相似度：
21)定义测地坐标g为从边界背景到当前体素在图像子集的无向图中的最短路径长度；根据步骤13)得到的无向图中连接相邻体素的边的权值，该权值对应相似度矩阵A中与该相邻体素对应的元素；并由该无向图估计得到图像子集体素的测地坐标，并在图像子集体素的测地坐标上以非监督聚类森林算法建立测地随机森林；
22)由所述测地随机森林获取体素相似度，更新步骤13)从表观随机森林中获取的体素相似度矩阵A；从测地随机森林中返回的体素相似度反映体素之间在测地距离意义上的相似；更新后的体素相似度
Figure FDA0002376212030000011
定义为式2：
Figure FDA0002376212030000021
其中，
Figure FDA0002376212030000022
是从测地随机森林各个决策树中得到的体素属于相同结构概率的平均值；φ表示在图像子集体素v之间的映射函数；
23)根据相似度矩阵更新测地坐标，依据更新后的测地坐标建立一个新的测地随机森林；对级联测地随机森林进行训练时，进行nk次的测地坐标更新得到nk层的级联测地随机森林，在第i次迭代中，相似度
Figure FDA0002376212030000023
其中
Figure FDA0002376212030000024
分别对应第i步迭代中体素对(vk,vl)的相似度，第i-1步迭代中体素对(vk,vl)的相似度，以及第i-1步迭代中测地随机森林得到的体素相似度；
(二)在线测试过程，包括：
1)从表观随机森林中获取图像子集Sr中的体素与图像子集St中的体素为对应体素的概率；并根据表观随机森林定义图像子集体素之间的相似度矩阵，初始定义连接图像子集体素的加权无向图，定义图像子集中体素的测地坐标；
2)根据级联测地随机森林返回的体素属于同一个结构的概率，对表观随机森林获取的体素相似度进行多次迭代更新；在第i次迭代中，由当前的相似度矩阵确定无向图中邻接边的权重进而估计测地坐标，将该测地坐标输入到对应的测地随机森林得到体素位于相同结构的概率，并以此更新图像子集之间的体素相似度；每次迭代更新的图像子集的相似度都用于估计图像子集之间体素的映射函数
Figure FDA0002376212030000025
其中对应体素
Figure FDA0002376212030000026
与
Figure FDA0002376212030000027
之间的相似度
Figure FDA0002376212030000028
最大；多数投票机制被用于从多次迭代得到的映射之中计算得到最后的映射函数φ；
3)根据图像子集的映射函数φ，利用正则化机制获取原始图像间稠密体素对应；
假定参照图像与目标图像之间的位移场记做C，其元素ci对应从参照图像Vr中的体素
Figure FDA0002376212030000029
到该体素在目标图像Vt中的对应体素之间的位移；根据给定图像子集的映射函数φ，预先确定一组体素对应，其中
Figure FDA00023762120300000210
通过式3估计整个图像对应的稠密位移场：
Figure FDA00023762120300000211
式3中，E为能量函数；C为图像对应的稠密位移场；μ1,μ2为常数系数；式3包括三项：
第一项是施加了位移场的变形后的参照图像与目标图像之间的图像灰度差异，其中
Figure FDA0002376212030000031
分别对应施加了位移场ci的参照图像的体素
Figure FDA0002376212030000032
的灰度与在目标图像中的对应体素
Figure FDA0002376212030000033
的灰度；Φ表示待求解的原始锥束CT图像Vr与Vt之间的稠密体素映射函数；当Vr中的体素在Vt中没有对应时，系数γi设置为0；否则，γi设置为1；
第二项中，
Figure FDA0002376212030000034
是由表观森林与级联测地森林获取的图像子集对应后，计算得到的体素
Figure FDA0002376212030000035
与其在目标图像子集中的对应体素之间的位移；ck表示在待求解稠密位移场中体素
Figure FDA0002376212030000036
的位移；通过第二项约束待求解的位移场与之前由表观森林与级联测地森林获取的图像子集对应一致；
第三项是平滑项，通过约束相邻的体素之间保持相似的位移以获取平滑的位移场，其中▽C是位移场的梯度；
求解式3，得到锥束CT图像Vr与Vt之间的稠密体素映射函数为：
Figure FDA0002376212030000037
其中，
Figure FDA0002376212030000038
分别为参照图像与目标图像Vt中的体素；ci为体素
Figure FDA0002376212030000039
的位移；
由此实现快速获取原始锥束CT图像之间体素的稠密对应。
2.如权利要求1所述自动建立方法，其特征是，所述根据锥束CT图像定义图像子集，具体采用三维尺度不变特征变换SIFT方法抽取卷积平滑差异图像的局部极值构造图像子集。
3.如权利要求1所述自动建立方法，其特征是，所述训练过程中，步骤11)随机从图像子集中选取体素构造表观随机森林中的决策树具体通过最大化信息增益得到节点的最优分裂策略，并基于该策略将节点的体素分到左子节点以及右子节点；在非监督的表观随机森林训练中，定义节点分裂中的信息增益为节点分裂后的左右子节点中体素的协方差矩阵的迹和的负数；由此建立决策树；根据建立的决策树估计并预测体素之间的相似度，具体是：将体素投入决策树的根节点后，依据节点中保存的最优分裂策略，该体素最终落入叶子节点中；再根据两个体素共有遍历路径的长度定义该体素对的相似度。
4.如权利要求1所述自动建立方法，其特征是，所述训练过程中，步骤21)建立测地随机森林时，在规则采样的锥束CT图像中，背景节点B定义在锥束CT图像的六个边界平面上。
5.如权利要求1所述自动建立方法，其特征是，所述在线测试过程中，步骤2)具体利用多数投票机制从多次迭代得到的映射中计算最后的映射函数φ。
6.如权利要求1所述自动建立方法，其特征是，所述在线测试过程中，步骤3)利用Levenberg-Marquardt算法求解式3的能量最小化问题。
CN201611184654.XA 2016-12-20 2016-12-20 锥束ct图像间体素稠密对应的自动建立方法 Active CN108205805B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201611184654.XA CN108205805B (zh) 2016-12-20 2016-12-20 锥束ct图像间体素稠密对应的自动建立方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201611184654.XA CN108205805B (zh) 2016-12-20 2016-12-20 锥束ct图像间体素稠密对应的自动建立方法

Publications (2)

Publication Number Publication Date
CN108205805A CN108205805A (zh) 2018-06-26
CN108205805B true CN108205805B (zh) 2020-06-02

Family

ID=62603267

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201611184654.XA Active CN108205805B (zh) 2016-12-20 2016-12-20 锥束ct图像间体素稠密对应的自动建立方法

Country Status (1)

Country Link
CN (1) CN108205805B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110222778B (zh) * 2019-06-11 2021-04-13 中国科学院自动化研究所 基于深度森林的在线多视角分类方法、系统、装置
CN114997278B (zh) * 2022-05-09 2023-04-07 浙江大学 基于计算机算法模型的工程数字化信息分析方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6421454B1 (en) * 1999-05-27 2002-07-16 Litton Systems, Inc. Optical correlator assisted detection of calcifications for breast biopsy
CN102945554A (zh) * 2012-10-25 2013-02-27 西安电子科技大学 基于学习和加速鲁棒surf特征的目标跟踪方法
CN104392250A (zh) * 2014-11-21 2015-03-04 浪潮电子信息产业股份有限公司 一种基于MapReduce的图像分类方法
CN104956397A (zh) * 2012-12-06 2015-09-30 西门子产品生命周期管理软件公司 基于自动空间上下文的3d图像中的多对象分段
CN105229699A (zh) * 2013-03-28 2016-01-06 外密景专家公司 基于医学图像评估血管网络的计算机实施的方法及其用途
CN105427325A (zh) * 2015-12-07 2016-03-23 苏州大学 基于随机森林和单调下降函数的肺肿瘤的自动分割方法
CN105528595A (zh) * 2016-02-01 2016-04-27 成都通甲优博科技有限责任公司 在无人机航拍图像中对输电线路绝缘子的识别定位方法
CN106204514A (zh) * 2015-04-30 2016-12-07 中国科学院深圳先进技术研究院 一种基于三维ct图像的肝脏定位方法及装置

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6421454B1 (en) * 1999-05-27 2002-07-16 Litton Systems, Inc. Optical correlator assisted detection of calcifications for breast biopsy
CN102945554A (zh) * 2012-10-25 2013-02-27 西安电子科技大学 基于学习和加速鲁棒surf特征的目标跟踪方法
CN104956397A (zh) * 2012-12-06 2015-09-30 西门子产品生命周期管理软件公司 基于自动空间上下文的3d图像中的多对象分段
CN105229699A (zh) * 2013-03-28 2016-01-06 外密景专家公司 基于医学图像评估血管网络的计算机实施的方法及其用途
CN104392250A (zh) * 2014-11-21 2015-03-04 浪潮电子信息产业股份有限公司 一种基于MapReduce的图像分类方法
CN106204514A (zh) * 2015-04-30 2016-12-07 中国科学院深圳先进技术研究院 一种基于三维ct图像的肝脏定位方法及装置
CN105427325A (zh) * 2015-12-07 2016-03-23 苏州大学 基于随机森林和单调下降函数的肺肿瘤的自动分割方法
CN105528595A (zh) * 2016-02-01 2016-04-27 成都通甲优博科技有限责任公司 在无人机航拍图像中对输电线路绝缘子的识别定位方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
The basics of maxillofacial cone beam computer tomography;Farman AG等;《Semin Orthod》;20091231;第15卷(第1期);第2-13页 *
应用锥束计算机体层摄影术评价颅颌面的不对称;张晓芸等;《北京大学学报(医学版)》;20130109;第45卷(第1期);第156-161页 *

Also Published As

Publication number Publication date
CN108205805A (zh) 2018-06-26

Similar Documents

Publication Publication Date Title
Harrison et al. Progressive and multi-path holistically nested neural networks for pathological lung segmentation from CT images
US20190294970A1 (en) Systems and methods for polygon object annotation and a method of training an object annotation system
CN106203432B (zh) 一种基于卷积神经网显著性图谱的感兴趣区域的定位系统
CN107633522B (zh) 基于局部相似性活动轮廓模型的脑部图像分割方法和系统
WO2019167883A1 (ja) 機械学習装置および方法
JP6325322B2 (ja) 医用画像処理装置、医用画像処理方法および医用画像処理プログラム
CN113826143A (zh) 特征点检测
CN110246580B (zh) 基于神经网络和随机森林的颅侧面影像分析方法和系统
CN111369525A (zh) 图像分析方法、设备和存储介质
CN110796691B (zh) 一种基于形状上下文和hog特征的异源图像配准方法
KR20210010920A (ko) 심장 관류 자기 공명 이미징을 사용하여 허혈성 심장 질환을 검출하기 위한 완전 정량적 픽셀 단위 심근 혈류 및 심근 관류 예비능 맵을 자동으로 생성하고 분석하는 방법 및 시스템
CN113298855B (zh) 基于自动勾画的图像配准方法
EP3012781A1 (en) Method and apparatus for extracting feature correspondences from multiple images
JP6431404B2 (ja) 姿勢推定モデル生成装置及び姿勢推定装置
CN112364881B (zh) 一种进阶采样一致性图像匹配方法
CN112102294A (zh) 生成对抗网络的训练方法及装置、图像配准方法及装置
Song et al. Multi-scale feature based land cover change detection in mountainous terrain using multi-temporal and multi-sensor remote sensing images
CN108205805B (zh) 锥束ct图像间体素稠密对应的自动建立方法
Tang et al. Retinal image registration based on robust non-rigid point matching method
CN108597589B (zh) 模型生成方法、目标检测方法及医学成像系统
CN111209946B (zh) 三维图像处理方法、图像处理模型训练方法及介质
CN115880358A (zh) 定位模型的构建方法、影像标志点的定位方法及电子设备
CN112581513B (zh) 锥束计算机断层扫描图像特征提取与对应方法
Saif et al. Computer Vision-based Efficient Segmentation Method for Left Ventricular Epicardium and Endocardium using Deep Learning
Chen et al. Fully-automatic landmark detection in skull X-ray 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