CN113781652B - 一种基于能量密度区域收缩的多层次概率重建方法 - Google Patents
一种基于能量密度区域收缩的多层次概率重建方法 Download PDFInfo
- Publication number
- CN113781652B CN113781652B CN202110935293.2A CN202110935293A CN113781652B CN 113781652 B CN113781652 B CN 113781652B CN 202110935293 A CN202110935293 A CN 202110935293A CN 113781652 B CN113781652 B CN 113781652B
- Authority
- CN
- China
- Prior art keywords
- index
- grid
- feasible region
- region
- matrix
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 239000011159 matrix material Substances 0.000 claims abstract description 38
- 230000008602 contraction Effects 0.000 claims abstract description 6
- 230000001172 regenerating effect Effects 0.000 claims description 7
- 230000001174 ascending effect Effects 0.000 claims description 3
- 238000005315 distribution function Methods 0.000 claims description 3
- 238000012217 deletion Methods 0.000 abstract description 3
- 230000037430 deletion Effects 0.000 abstract description 3
- 230000000877 morphologic effect Effects 0.000 abstract description 3
- 230000009467 reduction Effects 0.000 abstract description 2
- 210000001519 tissue Anatomy 0.000 description 35
- 238000004088 simulation Methods 0.000 description 30
- 230000008569 process Effects 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 8
- 230000005540 biological transmission Effects 0.000 description 7
- 230000003287 optical effect Effects 0.000 description 7
- 238000003325 tomography Methods 0.000 description 7
- 238000005415 bioluminescence Methods 0.000 description 6
- 230000029918 bioluminescence Effects 0.000 description 6
- 238000003384 imaging method Methods 0.000 description 6
- 238000010521 absorption reaction Methods 0.000 description 5
- 238000001914 filtration Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 238000012986 modification Methods 0.000 description 4
- 230000004048 modification Effects 0.000 description 4
- 210000000056 organ Anatomy 0.000 description 4
- 241001529936 Murinae Species 0.000 description 3
- 241000699666 Mus <mouse, genus> Species 0.000 description 3
- 238000009825 accumulation Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 230000003902 lesion Effects 0.000 description 3
- 241000699670 Mus sp. Species 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 210000003484 anatomy Anatomy 0.000 description 2
- 238000013473 artificial intelligence Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000005284 excitation Effects 0.000 description 2
- 210000002216 heart Anatomy 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 210000003734 kidney Anatomy 0.000 description 2
- 210000004185 liver Anatomy 0.000 description 2
- 210000004072 lung Anatomy 0.000 description 2
- 210000003205 muscle Anatomy 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 239000000523 sample Substances 0.000 description 2
- 210000002784 stomach Anatomy 0.000 description 2
- 241001465754 Metazoa Species 0.000 description 1
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- WYTGDNHDOZPMIW-RCBQFDQVSA-N alstonine Natural products C1=CC2=C3C=CC=CC3=NC2=C2N1C[C@H]1[C@H](C)OC=C(C(=O)OC)[C@H]1C2 WYTGDNHDOZPMIW-RCBQFDQVSA-N 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000009543 diffuse optical tomography Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000002224 dissection Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000003064 k means clustering Methods 0.000 description 1
- 238000004020 luminiscence type Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000010172 mouse model Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
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
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
- G06F18/232—Non-hierarchical techniques
- G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
- G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Software Systems (AREA)
- Bioinformatics & Computational Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Computer Graphics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Geometry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Image Analysis (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
本发明公开了一种基于能量密度区域收缩的多层次概率重建方法，获取生物组织表面光分布信息、网格节点坐标矩阵、网格内部四面体矩阵、网格内部四面体索引矩阵和系统矩阵作为输入进行重建；引入动态协方差可行域空间联合硬阈值收缩方法，克服了传统迭代算法只会缩小可行域区域的现象，避免单纯的空间缩小和节点删减，更加拟合光源整体空间，实现了光源在形态学上的高精度重建，针对多光源情形可以根据聚类结果各自进行可行域的运算。
Description
技术领域
本发明涉及医学成像技术领域，具体涉及的是一种基于能量密度区域收缩的多层次概率重建方法。
背景技术
光学分子影像(Optical Molecular Imaging，OMI)能够使医师以非侵入的方式在分子水平上观察疾病及其发展，相对于其它影像技术来说成本较为低廉，所以近年来被广泛应用于相关领域，逐渐成为一种完整的从诊断、监测到治疗的成像方式。
生物组织对其中传播的光子会有包括吸收、反射、折射、散射、发光、光化学、光声诸多复杂的效应，传输过程有漫反射、漫透射、弹道传输和蛇形传输，最终在表面采集到的光学信号相比于光源放射出的原始信号，无论在信号强度上面还是在信号出射方向上都会产生显著的变化。生物发光其主要在可见光和近可见光波段，波长较长，穿透能力相对较差，所以会有明显的散射、吸收现象。
为了更好的确定病灶的位置，获取光源的深度信息，依靠计算机技术的辅助，需要外部激发光的荧光分子断层成像(Fluorescence Molecular Tomography，FMT)技术与无需外部激发光的生物发光断层成像(Bioluminescence Tomography，BLT)先后被提出。这两种技术具有类似的原理，都是以生物组织内光传播的物理模型为基础，通过采集生物组织表面出射的光子，探测得到光子的强度分布信息，随后结合预先知道的生物组织内部的解剖结构以及相应各个组织的光学参数(如吸收系数、散射系数等)，通过使用不同的重建算法，反向计算内部的光源目标的空间位置以及信号强度分布，进而定量的展示生物组织内部的早期病灶变化情况。
可以发现，基于断层成像技术的三维重建方法可以有效的提升当前分子影像技术的适用范围和功能性。三维重建方法都是基于生物组织表面所出射的光学分子强度分布来针对组织内部的目标进行重建的。这个过程是一个明显的逆问题，有着很大的病态性，深层原因就是信息不对等。当组织内部的光源释放光子时，从光源表面直到探测器表面之间充满了光源所释放的光子，也就是能量分布信息，但重建过程要从结果去逆向的求解光源分布时，只能探测到生物组织表面的光子分布，而无法去收集生物组织中的传输中的光子分布数据，这导致实质上获得的只占所有能量分布信息的一少部分。如果将该问题类比于线性代数中的方程组问题的话，相当于整个方程组会拥有无数多个解，这就是三维重建逆问题病态性的实质。
为了解决这类问题，目前有着各种优化方法，其中迭代重建和人工智能方法是效果较为显著的两种方法。但是人工智能方法只能在同一个网格或相同结构的网格上进行重建，这一局限性严重制约了其应用范围；迭代重建则是根据一定的规则，单纯的缩减整个可行域空间，往往破坏了光源的几何特征和空间结构信息，导致结果出现较大的偏差，且往往采用最优的某一次迭代结果，甚至是仅仅采用最终次结果，当中间某次结果出现偏差时，错误便会累计到后面次，污染结果。
发明内容
针对上述技术中存在的不足，本发明提出了一种基于能量密度区域收缩的多层次概率重建方法，解决了现有方法的错误累积、删减节点和多个光源的问题。
为了实现上述目的，本发明采用的技术方案如下：
一种基于能量密度区域收缩的多层次概率重建方法，获取生物组织表面光分布信息、网格节点坐标矩阵、网格内部四面体矩阵、网格内部四面体索引矩阵和系统矩阵作为输入进行重建，包括：
步骤A：初始化网格节点全局索引得到包含所有网格节点索引的初始索引列Index_Sp；对初始索引列Index_Sp对应的网格节点形成的初始区域Index_Sp进行重建得到Xi，内部光源的能量强度值；
步骤B：结合Xi、网格节点坐标矩阵、网格内部四面体矩阵、网格内部四面体索引矩阵和系统矩阵得到当前次迭代的2范数错误率EL2和余弦相似度Ecos，并得到当前次迭代的全局概率权重值Xip；
其中，i表示迭代次数，取值为自然数；
步骤C：根据Xip得到当前次迭代的概率均值坐标XP、YP和ZP，网格节点坐标X的方差、Y的方差和Z的方差，XY的协方差、YZ的协方差和ZX的协方差，根据如下形式组装该次的协方差矩阵MCov：
根据协方差矩阵MCov，进一步得到MCov特征值和MＣov特征向量；
步骤D：利用MCov特征向量确定可行域的偏转角度和网格节点坐标的偏转角度，利用MCov特征值确定可行域的边长，可行域中包含的初始区域Index_Sp中的节点即为该次迭代更新得到的可行域Index_Sp_ROI；
步骤E：计算可行域Index_Sp_ROI中的网格节点坐标与概率均值坐标之间的距离误差，按照距离误差从小到大的升序排列可行域Index_Sp_ROI中的网格节点；
步骤F：衰减因子
更新可行域Index_Sp_ROI，以可行域Index_Sp_ROI的网格节点个数除以β2的商为更新标准，如果该商大于2，则下一次可行域的半边长度是2倍的MCov特征值；
如果该商小于等于2但是大于1，则下一次可行域的半边长度为1倍的MCov特征值；
如果该商小于等于1，则下一次可行域的半边长度为0.5倍的MCov特征值；
步骤G：按照初始区域Index_Sp中的网格节点的能量强度降序排列得到可行域Index_Sp_Descend，将当前可行域Index_Sp_Descend的网格节点个数作为初始节点数，根据公式更新β值和可行域Index_Sp_Descend的网格节点个数；
步骤H：新索引列Index_Sp为Index_Sp_ROI索引列拼接Index_Sp_Descend索引列，新索引列Index_Sp对应的网格节点形成新索引区域Index_Sp；新索引区域Index_Sp的网格节点个数小于等于1，或者迭代次数已经达到预设的最大迭代次数Lmax，根据每一次保存的XiP和全局概率权重值Pi，加权得到输出的X；否则返回步骤A；
可选的，计算得到全局概率权重值Pi的均值和标准差，形成高斯分布函数，仅保留一倍标准差的结果；
重新生成Pi，进而加权重新生成X。
本发明的基于能量密度区域收缩的多层次概率重建方法，根据所有次的迭代结果全局分析得到的，不会导致某些正确的节点在某些次被删去，也不会使得错误结果被接下来的迭代一直累积。在迭代过程中，本发明引入坐标协方差可行域和硬收缩联合方式作为每一次可行域划分的标准，根据每一次节点数量的动态变化速率，灵活调整可行域大小，并非一直向下递减收缩节点。本发明基于聚类方法，根据聚类出的初始能量团，各自计算内部的可行域，互不干扰，在考虑空间分布信息的情况下进行迭代。
附图说明
附图是用来提供对本公开的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本公开，但并不构成对本公开的限制。在附图中：
图1为南加州大学数字鼠仿真模型，所有的器官和组织均已标记；
图2为通过Amira软件剖分出来的不含光源的数字鼠仿真模型的.grid.am文件；
图3为通过COMSOL Multiphysics软件生成并剖分的半径为1mm的光源的展示图；
图4为通过Mose软件生物组织表面光分布信息的前向仿真结果，光源大小和形状为图3点光源，位置为(20，8，15)mm；
图5为通过Mose软件生物组织表面光分布信息的前向仿真结果，光源为图3的点光源，两个光源的空间位置分别为(22，8，15)mm和(15，8，15)mm；
图6为本发明单光源重建的整体迭代流程图；
图7为本发明双光源重建的整体迭代流程图；
图8为单光源仿真实验的实验结果，其中真实光源位置用箭头标记；
图9为图8沿着XOY平面在Z＝15mm处的切面图，圆圈代表真实光源位置，用箭头标记；
图10为图8对应的可行域收缩过程图；
图11为双光源仿真实验的实验结果，真实光源位置用箭头标记；
图12为图11沿着XOY平面在Z＝15mm处的切面图，黑色圆圈代表真实光源位置，用箭头标记；
图13为图11对应的可行域收缩过程图。
具体实施实例
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的重建方法，处理对象基于生物发光断层成像技术，以生物组织内光传播的物理模型为基础，通过采集生物组织表面出射的光子，探测得到光子的强度分布信息，随后结合预先知道的生物组织内部的解剖结构以及相应各个组织的光学参数(如吸收系数、散射系数等)，通过使用不同的重建算法，反向计算内部的光源目标的空间位置以及信号强度分布，进而定量的展示生物组织内部的早期病灶变化情况。具体涉及的是生物发光断层成像技术领域中的一种基于多层次概率模型的自适应重建迭代框架，该系统可用于快速寻找病变区域的准确位置，拥有定位精度高，不依赖于特定参数和特定算法等优点，可以显著提高目前重建算法的自身缺陷。
对于本发明中出现的名词解释如下：
生物组织表面光分布信息：指的是生物组织的表面能量分布向量。
网格指的是在模型重建过程中划分的多个立方体结构，立方体结构的8个顶点分别为网格节点，相邻的节点通过立方体结构的边长连接，在一个分析模型中，会出现多个网格，并含有多个网格节点。
网格节点坐标矩阵：指的是组成每一个四面体的四个顶点(节点)的坐标就是网格节点的坐标，其XYZ坐标可以组成网格节点坐标矩阵，
网格内部四面体矩阵：指的是每一个顶点坐标进行编号，可以得到组成每一个四面体的四个顶点的编号，这就是四面体矩阵。
网格内部四面体索引：指的是在三维重建中，使用四面体网格作为基本单元，四面体网格组成了整个生物组织的空间区域，对这些四面体进行编号，就可以得到四面体索引，
系统矩阵：指的是描述光子在组织中传播的数学关系与逻辑的数学因子组成的参数矩阵。
全局索引或索引：全局索引指的是初始第一次重建使用所有网格节点进行重建，这些网格节点的编号集合就是全局索引，每个网格节点编号就是一个索引；
索引区域：某一次的重建所使用的索引编号对应的网格节点组成的区域称为索引区域；可行域的结构等同于索引区域，通常认为可行域为索引区域可行域的边长，即为相邻节点在同一坐标方向之间的连线的长度，半边长即为边长的一半。
本发明的主要研究成果如下：(1)提出基于多级概率迭代的自适应重建框架，克服了传统迭代重建算法依赖于之前结果，会产生的错误累计现象，不依赖于某一次结果，本发明将每一次的结果都赋予其概率权重，最后统计各自的权重加权得到最终结果，创新性的保留每一次迭代的结果，根据其内部错误率灵活调整权重，避免了某些次结果丢失的问题；(2)引入动态协方差可行域空间联合硬阈值收缩方法，克服了传统迭代算法只会缩小可行域区域的现象，避免单纯的空间缩小和节点删减，更加拟合光源整体空间，实现了光源在形态学上的高精度重建，针对多光源情形可以根据聚类结果各自进行可行域的运算；(3)针对最后的重建结果，引入高斯分布模型，根据迭代过程中所有的错误率，将其按照高斯分布排布，自动删去不符合高斯分布结果的某些次重建信息，避免人为因素干扰。本发明在医学分子影像、重建方法等领域有重要的应用价值。本发明基于概率权重结果引入高斯分布过滤方法。根据高斯分布在均值±1倍标准差的区间内可信度为0.6826，保留一倍标准差的结果，滤除所有区间外次数的结果。滤除之后重新生成全局概率权重值，进而加权重建生成内部光源的能量强度值即为最终的结果。整个过程无需人为参数的设定。
该模型能够克服传统迭代重建方法依赖于上一次重建结果的限制，能够记录下所有迭代次的结果，根据其内部错误率自行分析权重，为每一次的结果赋予对应的概率值，在最后进行加权得到最终的重建结果。该方法根据每一次重建的结果，判断其空间内均值和协方差，同步使用硬阈值收缩方法判断当前可行域区域节点数，根据网格节点数量判断协方差可行域空间的边长，实现迭代过程中网格节点并非单纯减少，而是自动根据情况有增有减。针对多光源的情形，根据聚类结果划分各自的空间位置，并行进行收缩。该方法能够克服基于传统迭代重建算法单纯收缩可行域空间，导致光源形状信息损失的问题，在对目标形状的重建上，具有更高的精确度和泛化性。
实施例一：本发明的重建方法包括以下步骤：
在南加州大学数字鼠仿真模型上进行仿真实验数据集的构造。南该数字鼠仿真模型是由真实试验动物通过CT技术获取的成像切片数据中所提取的组织信息，经过断层成像重建手段获得去数字鼠模型，模型尺寸为长3.8cm，宽2.08cm，高3.5cm。为了降低计算复杂度，节约系统资源，去除小鼠的头部和尾部区域，只保留躯干部分。根据其内部组织结构，可划分为六个器官，及其肌体剩余组织，从上到下分别为心脏、肺、肝、胃、肾和肌肉。如图1所示。
下面详细说明构造一个前向能量仿真数据的操作过程。步骤如下：
步骤一：数据获取；
I、获得含有光源的.raw文件。
南加州大学数字鼠仿真模型是尺寸为38.0mm×20.8mm×35.0mm的.raw文件，其内部不含有光源。
II、获得整个仿真模型的.grid.am文件。
Amira软件是用于网格剖分的工具，对I中获得的.raw文件进行剖分，获得整个仿真模型的网格结构，文件格式为.grid.am文件。该步骤获得的示例整个仿真模型的.grid.am文件如图2所示。
III、获得仿真模型中所要使用的光源.mphtxt文件。
COMSOL Multiphysics软件是用来生成标准几何物理网格的工具，可以根据需要设定特定的几何形状，文件格式为.mphtxt文件。该步骤获得的光源网格文件如图3所示。
IV、获得生物组织表面光分布信息的前向仿真结果。
生物组织内部的特异性发光探针自发或受激辐射发出光子，光子在由内到外的传输过程中经历了反射、吸收、散射等多种物理光学过程，传输到生物体表面形成光分布信息。生物表面光分布信息与光源信息之间具有一对一的非线性映射关系。将II中获得的整个仿真模型的.grid.am文件通过前向仿真获得生物组织表面光分布信息。Mose软件是用来读入II中获得的网格文件，可以将.grid.am文件根据其内部的组织信息转化为.off文件，每种器官对应一个.off文件。随后Mose软件可以读入自身生成的.off文件，基于蒙特卡洛方法进行光子传输的模拟，从而进行前向能量仿真。前向仿真获得的文件为.CW格式，单光源和双光源生物组织表面光分布信息的前向仿真结果示例分别如图4和图5所示。
V、获得生物组织表面光分布信息，生成预运行数据。
使用Matlab编程语言从IV中获得的格式为.CW的前向仿真结果文件中提取出生物组织表面光分布信息，并相应的读入II中所生成的.grid.am网格文件。根据网格文件和前向结果文件，生成表面能量分布向量、网格节点坐标(X、Y和Z)矩阵、网格内部四面体矩阵、网格内部四面体索引矩阵和系统矩阵，将这些结果保存成Matlab对应的.mat文件，包含b_analytic.mat、Nodes.mat、InsideElement.mat、InsideElementIndex.mat、G.mat、triMat.mat共计6个mat文件。
步骤二、模型重建；将步骤一得到的.mat文件作为输入进行重建，重建过程主要包含如下几个步骤：
A.初始化全局索引，根据Nodes.mat文件，生成当前初始区域的所有节点的初始索引列Index_Sp；
B.针对初始索引列Index_Sp,使用A.N.Tikhonov的Solution of IncorrectlyFormulated Problems and the Regularization Method文章的Tikhonov正则化重建算法对当前的初始索引列Index_Sp对应的网格节点形成的初始区域Index_Sp进行重建得到Xi，内部光源的能量强度值；
C.根据当前得到的Xi，结合读入的G.mat文件和前向仿真能量.CW文件，可以计算出当前次的2范数错误率EL2和余弦相似度Ecos，根据这两个值结合前次的Xi可以计算出当前的全局概率权重值Pi，Pi＝0.5*(PL2+cos)，其中 进一步得到Xi对应的XiP＝i*i；i表示当前进行到的第i次迭代，i的取值为自然数；
当前次迭代的全局概率权重值Xip；
D.根据Xip得到当前次迭代的概率均值坐标XP、YP和ZP，网格节点坐标X的方差、Y的方差和Z的方差，XY的协方差、YZ的协方差和ZX的协方差；
E.根据如下形式组装该次的协方差矩阵MCov：
| X的方差 | XY的协方差 | XZ的协方差 |
| XY的协方差 | Y的方差 | YZ的协方差 |
| XZ的协方差 | YZ的协方差 | Z的方差 |
进一步得到MCov的特征值Val和特征向量Vec；
F.根据E中获得的Val和Vec，利用Vec确定可行域的偏转角度和网格节点坐标的偏转角度，利用Val确定可行域区间的每个方向上的边长，可以确定可行域区域的八个顶点，剖分八个顶点围城的区域，得到当前次的可行域区域，可行域中包含的初始区域Index_Sp中的节点即为该次迭代更新得到的可行域Index_Sp_ROI；
G.记录当前的可行域Index_Sp_ROI中的网格节点个数，保存为Cut_Num，采用Mohamed A.Naser的Improved bioluminescence and fluorescence reconstructionalgorithms using diffuse optical tomography,normalized data,and optimizedselection of the permissible source region文章，根据公式
计算出当前的衰减因子β,计算可行域Index_Sp_ROI中的网格节点坐标与概率均值坐标之间的距离误差，按照距离误差从小到大的升序排列可行域Index_Sp_ROI中的网格节点；
H.根据Cut_Num除以β2得到的商大小关系，如果该商大于2，则下一次的可行域空间的半边长是2倍的Val长度，如果该商小于等于2但是大于1，则下一次可行域空间半边长为1倍的Val长度，若该商小于等于1，则下一次可行域空间半边长为0.5倍的Val长度；
I.按照初始区域Index_Sp中的网格节点的能量强度降序排列得到可行域Index_Sp_Descend，将当前的网格节点个数作为初始节点数，根据公式更新β值和Cut_Num；
J.新索引列Index_Sp为索引列Index_Sp_ROI拼接索引列Index_Sp_Descend；
K.新索引列Index_Sp对应的网格节点形成新索引区域Index_Sp；新索引区域Index_Sp的网格节点个数小于等于1，或者迭代次数已经达到预设的最大迭代次数Lmax，根据每一次保存的XiP和全局概率权重值Pi，加权得到输出的X，即得到预测的内部光源的能量强度值；否则返回步骤A；
步骤三、实验结果的后处理阶段；
由于在步骤二中保留了每一次的结果，即使自适应概率权重值Pi会自适应的降低内部错误率较高的次数，也就是降低某些重建次影响力，仍然会有部分错误的结果被累计到最终结果X中。为此，根据所有结果的Pi可以得到其均值和标准差，进而可以画出高斯分布函数；根据高斯分布在均值±1倍标准差的区间内可信度为0.6826，保留一倍标准差的结果，滤除所有区间外次数的结果；滤除之后重新生成Pi，进而加权重新生成X即为最终的结果。
整体的流程图如图6和图7所示，图6为单光源流程，图7为双光源流程。针对双光源乃至多光源情形，只需要在步骤二的D之前，将光源结果进行K-Means聚类得到对应的空间，后续的步骤均按照不同光源各自区域分别处理即可。
实施例二：
为了验证所提框架的良好性能，本发明设置了单光源仿真实验。本发明使用南加州大学数字鼠仿真模型，模型尺寸为长3.8cm，宽2.08cm，高3.5cm。为了降低计算复杂度，节约系统资源，去除小鼠的头部和尾部区域，只保留躯干部分。根据其内部组织结构，可划分为六个器官，及其肌体剩余组织，从上到下分别为心脏、肺、肝、胃、肾和肌肉。该实验中测试样本的光源中心坐标设置在(20,8,15)mm处，光源半径设置为1mm。然后，通过前向仿真获得生物组织表面光分布信息。实验结果表明本发明的框架具有良好的目标定位和形态学重建性能。测试样本的重建光源中心坐标为(19.92，8.38mm，14.68)mm,定位误差LE为0.50mm。重建光源和真实光源之间的Dice系数为0.78，真是光源体积比上重建光源比率为0.67。单光源仿真实验的实验结果如图8和图9所示。可行域变化如图10所示。
实施例三：
为了验证多个光源下本发明的性能，本发明设置了双光源仿真实验。实验设置同实施实例1。该实验中测试样本的光源中心坐标设置在(22,8,15)mm和(15，8，15)mm处，光源半径设置为1mm。测试样本的重建光源中心坐标为(14.39,8.30,14.91)mm和(22.01,8.76,14.95)mm,定位误差LE分别为0.66mm和0.54mm。重建光源和真实光源之间的Dice系数分别为为0.69和0.63，真是光源体积比上重建光源比率分别为0.62和0.92。单光源仿真实验的实验结果如图11和图12所示。可行域变化如图13所示。
通过上述两个实施案例可以表面，本发明提出的方法具有良好的目标定位和光源几何结构恢复性能，且多光源的实验也表面了该框架具有识别多光源并很好的复原的能力。
对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。
以上结合附图详细描述了本公开的优选实施方式，但是，本公开并不限于上述实施方式中的具体细节，在本公开的技术构思范围内，可以对本公开的技术方案进行多种简单变型，这些简单变型均属于本公开的保护范围。
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合，为了避免不必要的重复，本公开对各种可能的组合方式不再另行说明。
此外，本公开的各种不同的实施方式之间也可以进行任意组合，只要其不违背本公开的思想，其同样应当视为本公开所公开的内容。
Claims (2)
1.一种基于能量密度区域收缩的多层次概率重建方法，其特征在于，获取生物组织表面光分布信息、网格节点坐标矩阵、网格内部四面体矩阵、网格内部四面体索引矩阵和系统矩阵作为输入进行重建，包括：
步骤A：初始化网格节点全局索引得到包含所有网格节点索引的初始索引列Index_Sp；对初始索引列Index_Sp对应的网格节点形成的初始区域Index_Sp进行重建得到Xi，内部光源的能量强度值；
步骤B：结合Xi、网格节点坐标矩阵、网格内部四面体矩阵、网格内部四面体索引矩阵和系统矩阵得到当前次迭代的2范数错误率EL2和余弦相似度Ecos，并得到当前次迭代的全局概率权重值Xip；
其中，i表示迭代次数，取值为自然数；
步骤C：根据Xip得到当前次迭代的概率均值坐标XP、YP和ZP，网格节点坐标X的方差、Y的方差和Z的方差，XY的协方差、YZ的协方差和ZX的协方差，根据如下形式组装该次的协方差矩阵MCov：
根据协方差矩阵MCov，进一步得到MCov特征值和MCov特征向量；
步骤D：利用MCov特征向量确定可行域的偏转角度和网格节点坐标的偏转角度，利用MCov特征值确定可行域的边长，可行域中包含的初始区域Index_Sp中的节点即为该次迭代更新得到的可行域Index_Sp_ROI；
步骤E：计算可行域Index_Sp_ROI中的网格节点坐标与概率均值坐标之间的距离误差，按照距离误差从小到大的升序排列可行域Index_Sp_ROI中的网格节点；
步骤F：衰减因子
更新可行域Index_Sp_ROI，以可行域Index_Sp_ROI的网格节点个数除以β2的商为更新标准，如果该商大于2，则下一次可行域的半边长度是2倍的MCov特征值；
如果该商小于等于2但是大于1，则下一次可行域的半边长度为1倍的MCov特征值；
如果该商小于等于1，则下一次可行域的半边长度为0.5倍的MCov特征值；
步骤G：按照初始区域Index_Sp中的网格节点的能量强度降序排列得到可行域Index_Sp_Descend，将当前可行域Index_Sp_Descend的网格节点个数作为初始节点数，根据公式更新β值和可行域Index_Sp_Descend的网格节点个数；
步骤H：新索引列Index_Sp为Index_Sp_ROI索引列拼接Index_Sp_Descend索引列，新索引列Index_Sp对应的网格节点形成新索引区域Index_Sp；新索引区域Index_Sp的网格节点个数小于等于1，或者迭代次数已经达到预设的最大迭代次数Lmax，根据每一次保存的XiP和全局概率权重值Pi，加权得到输出的X；否则返回步骤A；
2.根据权利要求1所述的基于能量密度区域收缩的多层次概率重建方法，其特征在于，计算得到全局概率权重值Pi的均值和标准差，形成高斯分布函数，仅保留一倍标准差的结果；
重新生成Pi，进而加权重新生成X。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110935293.2A CN113781652B (zh) | 2021-08-16 | 2021-08-16 | 一种基于能量密度区域收缩的多层次概率重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110935293.2A CN113781652B (zh) | 2021-08-16 | 2021-08-16 | 一种基于能量密度区域收缩的多层次概率重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113781652A CN113781652A (zh) | 2021-12-10 |
| CN113781652B true CN113781652B (zh) | 2024-04-12 |
Family
ID=78837784
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110935293.2A Active CN113781652B (zh) | 2021-08-16 | 2021-08-16 | 一种基于能量密度区域收缩的多层次概率重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113781652B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105326475A (zh) * | 2015-09-16 | 2016-02-17 | 西北大学 | 一种基于多光源分辨的生物发光断层成像重建方法 |
| US10002286B1 (en) * | 2015-04-28 | 2018-06-19 | Carnegie Mellon University | System and method for face recognition robust to multiple degradations |
| CN111771229A (zh) * | 2019-01-30 | 2020-10-13 | 百度时代网络技术（北京）有限公司 | 用于自动驾驶车辆的点云重影效果检测系统 |
| CN112055805A (zh) * | 2019-01-30 | 2020-12-08 | 百度时代网络技术（北京）有限公司 | 用于自动驾驶车辆的点云登记系统 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7978932B2 (en) * | 2007-08-02 | 2011-07-12 | Mauna Kea Technologies | Robust mosaicing method, notably with correction of motion distortions and tissue deformations for in vivo fibered microscopy |
| US9524552B2 (en) * | 2011-08-03 | 2016-12-20 | The Regents Of The University Of California | 2D/3D registration of a digital mouse atlas with X-ray projection images and optical camera photos |
| US10542243B2 (en) * | 2018-04-10 | 2020-01-21 | Intel Corporation | Method and system of light source estimation for image processing |
-
2021
- 2021-08-16 CN CN202110935293.2A patent/CN113781652B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10002286B1 (en) * | 2015-04-28 | 2018-06-19 | Carnegie Mellon University | System and method for face recognition robust to multiple degradations |
| CN105326475A (zh) * | 2015-09-16 | 2016-02-17 | 西北大学 | 一种基于多光源分辨的生物发光断层成像重建方法 |
| CN111771229A (zh) * | 2019-01-30 | 2020-10-13 | 百度时代网络技术（北京）有限公司 | 用于自动驾驶车辆的点云重影效果检测系统 |
| CN112055805A (zh) * | 2019-01-30 | 2020-12-08 | 百度时代网络技术（北京）有限公司 | 用于自动驾驶车辆的点云登记系统 |
Non-Patent Citations (2)
| Title |
|---|
| 一种基于模糊C均值聚类的新型单视图切伦科夫发光断层成像方法;贺小伟;卫潇;曹欣;侯榆青;;光学学报;20180312(07);全文 * |
| 余景景 ; 王海玉 ; 李启越 ; .结合迭代收缩可行域的单视图多光谱生物发光断层成像.光学学报.(12),全文. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113781652A (zh) | 2021-12-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111709953B (zh) | Ct影像的肺叶段分割中的输出方法、装置 | |
| CN109191564B (zh) | 基于深度学习的激发荧光断层成像三维重建方法 | |
| CN108451508B (zh) | 基于多层感知机的生物自发荧光三维成像方法 | |
| CA2481650A1 (en) | Modification of the normalized difference method for real-time optical tomography | |
| Li et al. | Refraction corrected transmission ultrasound computed tomography for application in breast imaging | |
| CN110327018B (zh) | 稀疏度自适应组正交匹配追踪的激发荧光断层重建方法 | |
| CN112767407A (zh) | 一种基于级联门控3DUnet模型的CT图像肾脏肿瘤分割方法 | |
| CN103300829A (zh) | 一种基于迭代重加权的生物自发荧光断层成像方法 | |
| CN104921706A (zh) | 基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法 | |
| CN112734872A (zh) | 一种基于多波长并发重建的荧光分子断层成像方法及系统 | |
| Li et al. | SAR-CGAN: Improved generative adversarial network for EIT reconstruction of lung diseases | |
| CN115272590B (zh) | 光学传输参数空间分布的重建方法、装置、系统及介质 | |
| CN115868923A (zh) | 基于展开式循环神经网络的荧光分子断层成像方法及系统 | |
| CN113781652B (zh) | 一种基于能量密度区域收缩的多层次概率重建方法 | |
| CN115760874A (zh) | 一种基于联合空间域的多尺度U-Net医学图像分割方法 | |
| Sun et al. | An iterative gradient convolutional neural network and its application in endoscopic photoacoustic image formation from incomplete acoustic measurement | |
| CN112465924A (zh) | 一种基于多特征融合的快速医学图像重构方法 | |
| CN111751343A (zh) | 基于瓶颈残差gan的介观荧光分子层析成像方法及系统 | |
| An et al. | Application of machine learning method in optical molecular imaging: a review | |
| Wang et al. | Adaptive machine learning method for photoacoustic computed tomography based on sparse array sensor data | |
| CN116712038A (zh) | 基于螺旋交错稀疏采样的多光谱光声断层成像系统及方法 | |
| CN113288188B (zh) | 基于分组注意力残差网络的锥束x射线发光断层成像方法 | |
| Lyu et al. | Photoacoustic digital brain: numerical modelling and image reconstruction via deep learning | |
| Long et al. | Spatial and channel attention-based conditional Wasserstein GAN for direct and rapid image reconstruction in ultrasound computed tomography | |
| CN107374588B (zh) | 一种基于同步聚类的多光源荧光分子断层成像重建方法 |
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 |