CN116644621A - 一种全光谱后向散射漫反射率的模拟方法 - Google Patents

一种全光谱后向散射漫反射率的模拟方法 Download PDF

Info

Publication number
CN116644621A
CN116644621A CN202310929532.2A CN202310929532A CN116644621A CN 116644621 A CN116644621 A CN 116644621A CN 202310929532 A CN202310929532 A CN 202310929532A CN 116644621 A CN116644621 A CN 116644621A
Authority
CN
China
Prior art keywords
spf
scattering
photon
wavelength
bias
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.)
Granted
Application number
CN202310929532.2A
Other languages
English (en)
Other versions
CN116644621B (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.)
Always Wuxi Medical Technology Co ltd
Original Assignee
Always Wuxi Medical 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 Always Wuxi Medical Technology Co ltd filed Critical Always Wuxi Medical Technology Co ltd
Priority to CN202310929532.2A priority Critical patent/CN116644621B/zh
Publication of CN116644621A publication Critical patent/CN116644621A/zh
Application granted granted Critical
Publication of CN116644621B publication Critical patent/CN116644621B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

本发明涉及全光谱后向散射漫反射率模拟技术领域，公开了一种全光谱后向散射漫反射率的模拟方法，包括S1、根据TT SPF理论推导偏置散射函数TT‑B SPF；S2、求解TT SPF中的波长依赖的参数集合；S3、部署基于TT SPF、TT‑B SPF和IS的模拟方法以计算得到F‑BDR。本发明是基于一种与米氏SPF拟合更好、对后向散射事件模拟更加准确的TT SPF，同时推导TT‑B SPF，并根据精确的米氏散射理论计算波长依赖的散射参数，不仅适用于BDR的快速模拟，也可模拟全光谱的特性，并且本发明所述的图像模拟算法考虑每个波长模拟的单独参数计算，使得模拟出的结果具有波长依赖性，实现了F‑BDR和FD‑OCT图像的快速准确模拟。

Description

一种全光谱后向散射漫反射率的模拟方法
技术领域
本发明涉及全光谱后向散射漫反射率模拟技术领域，具体为一种全光谱后向散射漫反射率的模拟方法。
背景技术
全光谱后向散射漫反射率（Full-spectrum backscattered diffusereflectance, F-BDR）提供了一种对生物组织结构、微观特性以及光谱特性的非侵入性检测指标，已经成为了多种医学诊疗和生物成像技术的理论基础。蒙特卡洛（Monte Carlo,MC）模拟方法长期以来被认为是模拟光与组织相互作用的“金标准”，近些年来也被应用于F-BDR的精确理论建模。然而基于MC的F-BDR的正向模拟方法始终无法被推广使用，有限的模拟F-BDR的方法，更是大多局限于单条BDR（Backscattered diffuse reflectance，后向散射漫反射率）线的模拟，难以扩展到面扫描（B-scan）和体扫描的模拟，这进一步限制了基于F-BDR的成像应用的模拟。究其原因，这是因为精确F-BDR的MC模拟基于米氏散射理论，它依赖于具有极低概率被采集的后向散射光子，同时还需要进行全光谱中每个波长的单独模拟，所以模拟高质量的F-BDR的计算耗时极长，因此现有的基于蒙特卡洛的模拟方法存在受制于不精确的散射相函数模型、极低的后向散射光子探测效率以及全光谱模拟的额外复杂度而导致无法快速精确地模拟全光谱后向散射漫反射率的问题。
针对相关技术中的问题，目前尚未提出有效的解决方案。
发明内容
针对现有技术的不足，本发明提供了一种全光谱后向散射漫反射率的模拟方法，具备适用于BDR的快速模拟，可模拟全光谱的特性，并且本发明所述的图像模拟算法考虑每个波长模拟的单独参数计算，使得模拟出的结果具有波长依赖性等优点，解决了现有的基于蒙特卡洛的模拟方法存在受制于不精确的散射相函数模型、极低的后向散射光子探测效率以及全光谱模拟的额外复杂度而导致无法快速精确地模拟全光谱后向散射漫反射率的问题。
为解决上述现有的基于蒙特卡洛的模拟方法存在受制于不精确的散射相函数模型、极低的后向散射光子探测效率以及全光谱模拟的额外复杂度而导致无法快速精确地模拟全光谱后向散射漫反射率的技术问题，本发明提供如下技术方案：
一种全光谱后向散射漫反射率的模拟方法，包括
S1、根据TT SPF 理论推导偏置散射函数TT-B SPF；
S1.1、定义光子的偏置散射条件；
在光子包根据偏置散射函数TT-B SPF进行偏置散射的假设下，定义偏置散射条件为：偏置散射后的光子行进方向为，探测器方向为/>，则偏置散射夹角/>的余弦值，即/>，偏置散射角/>根据TT-B SPF采样获得，具体采样方式为对TT-B SPF实行拒绝-接受采样方法以得到角度值[Casella G,Robert C P, Wells M T. Generalized accept-reject sampling schemes[J]. LectureNotes-Monograph Series, 2004: 342-347.];
S1.2、根据TT SPF和偏置散射条件推导TT-B SPF的偏置单向归一化系数；
TT SPF由两个分别由前向散射和后向散射主导的单向 SPF线性组合而成，写为，；其中/>是TT SPF,是TT SPF的参数集合向量，/>表示前向散射SPF，/>表示后向散射SPF，/>是散射角度，满足/>，/>和为/>为两个方向的各向异性因子，分别控制了散射向前向和后向偏向的程度，/>和/>为增强因子，分别对/>和/>控制的效应进行增强，/>是权重因子，对前向散射和后向散射进行平衡，并满足/>；对于某一方向散射主导的单向 SPF，/>函数表达式为， ；其中/>分别为各向异性因子和增强因子，x=f代指前向，x=b代指后向，/>是单向归一化系数，表示为，/>根据偏置散射条件将以上归一化系数进行变形，推导出TT-B SPF中偏置单向归一化系数/>应写作，/>；其中/>，/>为偏置各向异性因子，/>为偏置增强因子，两者取值范围同/>与/>, x=f代指前向，x=b代指后向;
S1.3、根据偏置单向归一化系数推导TT-B SPF;将偏置单向归一化系数替换/>定义式中的单向归一化系数，得到偏置单向 SPF表达式为，/> ，，其中/>为偏置单向 SPF的普遍表示，x=f代指前向，x=b代指后向，为偏置散射夹角;随后线性组合前向散射和后向散射的单向 SPF，得到TT-BSPF的表达式为，/> ；其中/>为TT-B SPF，为偏置函数的权重因子，/>是TT-B SPF的参数集合向量。
S2、求解TT SPF中的波长依赖的参数集合；
S3、部署基于TT SPF（Two-term SPF，两项散射相函数）、TT-B SPF（Two-termbiased SPF，偏置两项散射相函数）和IS（Importance sampling，重要性采样） 的模拟方法以计算得到F-BDR；
其中步骤S3中部署模拟方法以得到F-BDR的方法为，
S3.1、加载对应当前波长的光子包并确定波长相关的散射系数；
S3.2、光子移动、吸收与边界判断；
S3.3、基于IS的光子散射过程；
S3.4、判断光子是否存活；
S3.5、判断光子是否被探测；
S3.6、F-BDR计算与FD-OCT（Fourier-domain optical coherence tomography，傅里叶域光学相干断层成像）图像生成。
优选地，所述步骤S2中求解TT SPF中的波长依赖的参数集合的方法为：
S2.1、根据TT SPF中波长依赖的参数集合中的五个参数和当前波长处的米氏SPF对TT SPF进行拟合，以获取优化后的TT SPF中波长依赖的参数集合中的五个参数；
使用非线性最小二乘法建立优化方程，通过优化调整TT SPF的五个参数来拟合当前波长下的米氏SPF和TT SPF，当优化方程达到最小值时，输出此时优化后的TT SPF中波长依赖的参数集合中的五个参数；
具体方法为：在波长为时的米氏SPF定义为/>，数学表达式为：；其中/>是散射角度，/>是粒子尺寸参数，被定义为/>，/>为背景介质的折射率，/>为散射粒子半径，/>是散射效率，/>和/>是散射振幅矩阵的分量。TT SPF在波长为/>时的参数可以通过参数集合向量表示为/>；使用非线性最小二乘法的优化方程可以表示为，/>；其中，/>；是SPF的缩放因子，/>是缩放因子中的权重函数，表示为；其中/>是在散射角/>时的值。通过优化调整参数集合向量/>中的五个参数，使得上述优化方程达到最小值，即此时TT SPF和当前波长的米氏SPF最为接近。使得上述优化方程达到最小值的参数集合向量为最优参数集合向量/>，该向量的五个参数为最优的TTSPF中波长依赖的五个参数。S2.2、使用非线性最小二乘优化得到的优化后的TT SPF中波长依赖的五个参数构造波长依赖的TT SPF组和TT-B SPF组；
通过重复步骤S2.1对所有波长对应的TT SPF和米氏SPF进行拟合，可以得到每个波长下优化后的TT SPF的波长依赖的参数集合，然后将这些参数集合中的五个参数带入TTSPF的表达式构造一组波长依赖的TT SPF；之后将该组TT SPF中的五个参数重用到对应的TT-B SPF中以得到一组波长依赖的TT-B SPF，即该组TT-B SPF中对应于波长的参数集合向量。优选地，关于加载对应当前波长的光子包并确定波长相关的散射系数的方法为，加载对应当前模拟波长/>的光子包/>，光子包以垂直仿体上表面的方向射入仿体；此时该光子包对应的仿体散射系数/>为/>；其中/>是仿体中粒子的数字密度，/>是粒子尺寸参数，/>为散射粒子半径，/>是散射效率。优选地，关于光子移动、吸收与边界判断的具体操作为，光子移动：光子包被加载后，在仿体中单次移动一随机步长，步长/>可表示为/>；其中/>为光子包所在当前仿体层的吸收系数，/>为0到1之间均匀分布的随机数；光子边界判断：移动结束后，判断是否击中下一仿体单元的边界，若击中，则进行反射或透射的判断，若经过反射或透射后被接收，则跳转步骤S3.5；若并未击中边界，则当移动结束后结束，根据/>执行吸收；光子吸收：移动结束后并未击中边界，光子包/>的能量根据/>被仿体吸收，导致能量衰减。优选地，关于基于IS的散射，若光子包并未逸出仿体边界，且并未被探测器接收，则当移动结束后，进行散射事件；在光子移动过程中，三维前进方向定义为/>，探测器方向为，散射后的新前进方向为/>。对于任一波长/>下的MC模拟，所有光子都从探测器垂直射入组织，即/>。优选地，关于基于IS的散射，基于IS确定散射角的流程包括第一次散射事件、分裂、额外偏置散射、权重更新。
优选地，关于存活判断的方法，
若光子散射后逸出仿体边界，则光子寿命终止，加载下一光子；若光子散射后并未逸出仿体边界，且能量大于事先设定的能量阈值，则光子继续存活，并返回步骤S3.2；若光子散射后并未逸出边界，但能量小于事先设定的能量阈值，则进行俄罗斯轮盘赌，即随机生成一个在[0,1]均匀分布的随机数，若该随机数小于等于预先设定好的机会概率，则视为赌赢，否则视为赌输，若光子赌赢，则能量增大十倍并返回步骤S3.2，若赌输，则光子寿命终止，加载下一光子。
优选地，关于探测判断的方法为：若光子在边界处被探测器接收，记录光子包的剩余能量作为对应于波长的BDR />，记录光子包的行走路径长作为光程/>，随后终止该光子包的寿命，加载下一光子包。优选地，关于F-BDR计算与FD-OCT图像生成的方法为，当所有波长的光子包经历步骤S3.1到S3.5，记录下的光子包的剩余能量按照波长和光程进行二维排列，得到全光谱下的F-BDR图：/>, />，/>为光谱最小波长，/>为光谱最大波长；
根据该F-BDR图使用信号计算方法重建FD-OCT图像，并将其作为F-BDR图的另一表现形式。优选地，根据该F-BDR图使用信号计算方法重建FD-OCT图像具体方法为，使用F-BDR值作为强度，光程作为相位，构建干涉条纹并进行傅里叶反变换重建FD-OCT图像。
与现有技术相比，本发明提供了一种全光谱后向散射漫反射率的模拟方法，具备以下有益效果：
本发明是基于一种与米氏SPF拟合更好、对后向散射事件模拟更加准确的TT SPF，同时推导TT-B SPF根据精确的米氏散射理论计算波长依赖的散射参数，不仅适用于BDR的快速模拟，也可模拟全光谱的特性，并且本发明所述的图像模拟算法考虑每个波长模拟的单独参数计算，使得模拟出的结果具有波长依赖性，实现了F-BDR和FD-OCT图像的快速准确模拟。
附图说明
图1为本发明的全光谱后向散射漫反射率的模拟方法流程示意图；
图2为本发明的单条A线的F-BDR模拟结果图；
图3为本发明的面扫描的F-BDR模拟结果图；
图4为本发明的面扫描的FD-OCT图像模拟结果图；
图5为本发明的体扫描的FD-OCT图像模拟结果图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
正如背景技术所介绍的，现有技术中存在的不足，为了解决如上的技术问题，本申请提出了一种全光谱后向散射漫反射率的模拟方法。
请参阅图1-5，一种全光谱后向散射漫反射率的模拟方法，
S1、根据TT SPF 理论推导偏置散射函数TT-B SPF；
S1.1、定义光子的偏置散射条件；在光子包根据偏置散射函数TT-B SPF进行偏置散射的假设下，定义偏置散射条件为：偏置散射后的光子行进方向为，探测器方向为/>，则偏置散射夹角/>的余弦值/>，即/>，偏置散射角/>根据TT-B SPF采样获得，具体采样方式为对TT-B SPF实行拒绝-接受采样方法以得到角度值[Casella G, Robert C P, Wells M T. Generalized accept-rejectsampling schemes[J]. Lecture Notes-Monograph Series, 2004: 342-347.]。S1.2、根据TT SPF和偏置散射条件推导TT-B SPF的偏置单向归一化系数；
TT SPF由两个分别由前向散射和后向散射主导的单向 SPF线性组合而成，可以写为，；其中/>是TT SPF,是TT SPF的参数集合向量，/>表示前向散射SPF，/>表示后向散射SPF，/>是散射角度，满足/>，/>和/>为两个方向的各向异性因子，分别控制了散射向前向和后向偏向的程度，/>和/>为增强因子，分别对/>和/>控制的效应进行增强，/>是权重因子，对前向散射和后向散射进行平衡，并满足/>。对于某一方向散射主导的单向SPF，/>，函数表达式为，/> ；其中/>分别为各向异性因子和增强因子，x=f代指前向，x=b代指后向，/>是单向归一化系数，可以表示为，；根据偏置散射条件将以上归一化系数进行变形，推导出TT-B SPF中偏置单向归一化系数/>应写作，；其中/>，/>为偏置各向异性因子，/>为偏置增强因子，两者取值范围同/>与/>，x=f代指前向，x=b代指后向。S1.3、根据偏置单向归一化系数推导TT-B SPF；将偏置单向归一化系数/>替换定义式中的单向归一化系数，可得偏置单向 SPF表达式/>为，/> ；其中/>为偏置单向 SPF的普遍表示，x=f代指前向，x=b代指后向，/>为偏置散射夹角；随后线性组合前向散射和后向散射的单向 SPF，得到TT-B SPF的表达式为， ；其中/>为TT-B SPF，为偏置函数的权重因子，/>是TT-B SPF的参数集合向量。S2、求解TT SPF中的波长依赖的参数集合；
S2.1、根据TT SPF中波长依赖的参数集合中的五个参数和当前波长处的米氏SPF对TT SPF进行拟合，以获取优化后的TT SPF中波长依赖的参数集合中的五个参数；
使用非线性最小二乘法建立优化方程，通过优化调整TT SPF的五个参数来拟合当前波长下的米氏SPF和TT SPF，当优化方程达到最小值时，输出此时优化后的TT SPF中波长依赖的参数集合中的五个参数；
具体方法为：在波长为时的米氏SPF定义为/>，数学表达式为：；其中/>是散射角度，/>是粒子尺寸参数，被定义为/>，/>为背景介质的折射率，/>为散射粒子半径，/>是散射效率，/>和/>是散射振幅矩阵的分量。TT SPF在波长为/>时的参数可以通过参数集合向量表示为/>；使用非线性最小二乘法的优化方程可以表示为，/>；其中，/>；是SPF的缩放因子，/>是缩放因子中的权重函数，表示为；其中/>是在散射角/>时的值。通过优化调整参数集合向量/>中的五个参数，使得上述优化方程达到最小值，即此时TT SPF和当前波长的米氏SPF最为接近。使得上述优化方程达到最小值的参数集合向量为最优参数集合向量/>，该向量的五个参数为最优的TTSPF中波长依赖的五个参数。S2.2、使用非线性最小二乘优化得到的优化后的TT SPF中波长依赖的五个参数构造波长依赖的TT SPF组和TT-B SPF组；
通过重复步骤S2.1对所有波长对应的TT SPF和米氏SPF进行拟合，可以得到每个波长下优化后的TT SPF的波长依赖的参数集合，然后将这些参数集合中的五个参数带入TTSPF的表达式构造一组波长依赖的TT SPF；之后将该组TT SPF中的五个参数重用到对应的TT-B SPF中以得到一组波长依赖的TT-B SPF，即该组TT-B SPF中对应于波长的参数集合向量；S3、部署基于TT SPF、TT-B SPF和IS 的模拟方法以计算得到F-BDR；
其中步骤S3中部署模拟方法以得到F-BDR的方法为，
S3.1、加载对应当前波长的光子包并确定波长相关的散射系数；加载对应当前模拟波长的光子包/>，光子包以垂直仿体上表面的方向射入仿体；此时该光子包对应的仿体散射系数/>为/>；其中/>是仿体中粒子的数字密度，/>是粒子尺寸参数，/>为散射粒子半径，/>是散射效率。S3.2、光子移动、吸收与边界判断；光子移动：光子包被加载后，在仿体中单次移动一随机步长，步长/>可表示为；其中/>为光子包所在当前仿体层的吸收系数，/>为0到1之间均匀分布的随机数；光子边界判断：移动结束后，判断是否击中下一仿体单元的边界，若击中，则进行反射或透射的判断，若经过反射或透射后被接收，则跳转步骤S3.5；若并未击中边界，则当移动结束后结束，根据/>执行吸收；光子吸收：移动结束后并未击中边界，光子包的能量根据/>被仿体吸收，导致能量衰减。S3.3、基于IS的光子散射过程；
若光子包并未逸出仿体边界，且并未被探测器接收，则当移动结束后，进行散射事件；在光子移动过程中，三维前进方向定义为，探测器方向为/>，散射后的新前进方向为/>。对于任一波长/>下的MC模拟，所有光子都从探测器垂直射入组织，即/>。基于IS确定散射角的流程如下,1、第一次偏置散射事件。若当前光子从未经历过偏置散射，并且远离探测器行进，即/>，它将经历第一次偏置散射。偏置散射后的行进方向与探测器方向的夹角，即从TT-B SPF中采样得到的偏置角为，此时等效的由TT SPF采样得到的真实散射角为。第一次偏置散射后的似然比计算如下，/>；2、分裂。本发明使用和《Lima I T, Kalra A, Hernández-Figueroa H E, et al. Fastcalculation of multipath diffusive reflectance in optical coherencetomography[J]. Biomedical optics express, 2012, 3(4): 692-700.》类似的光子分裂机制：当一个光子已经被第一次偏置散射到行进方向/>，另一个非偏置的光子将在原来的位置生成，并带有较大的似然比表示为/>，两个光子被视为由原始的一个光子分裂而成；随后非偏置光子将经历一次非偏置散射，其散射角度/>由/>采样得到。而与《Lima I T, Kalra A, Hernández-Figueroa H E, et al. Fast calculationof multipath diffusive reflectance in optical coherence tomography[J].Biomedical optics express, 2012, 3(4): 692-700.》不同的是，本发明中该散射角度/>将被持续采样，直到其小于90度，这避免了该非偏置光子被直接散射回探测器，这一点在之前的文献没有被考虑。3、额外偏置散射。当经历过第一次偏置散射的光子沿着远离探测器的方向前进时，即/>；它将会以概率/>进行基于TT-B SPF的额外偏置散射，以的概率进行基于TT SPF的非偏置散射。无论进行何种散射，该次散射的似然比可以计算为，/>；4、权重更新。在经历基于IS的散射过程后，所有事件计算的似然度将依次相乘得到总似然度/>。在该光子最终被探测器探测前，其权重更新为/>。S3.4、判断光子是否存活；
若光子散射后逸出仿体边界，则光子寿命终止，加载下一光子；若光子散射后并未逸出仿体边界，且能量大于事先设定的能量阈值，则光子继续存活，并返回步骤S3.2；若光子散射后并未逸出边界，但能量小于事先设定的能量阈值，则进行俄罗斯轮盘赌，即随机生成一个在[0,1]均匀分布的随机数，若该随机数小于等于预先设定好的机会概率，则视为赌赢，否则视为赌输，若光子赌赢，则能量增大十倍并返回步骤S3.2，若赌输，则光子寿命终止，加载下一光子。S3.5、判断光子是否被探测；
若光子在边界处被探测器接收，记录光子包的剩余能量作为对应于波长的BDR ，记录光子包的行走路径长作为光程/>，随后终止该光子包的寿命，加载下一光子包。S3.6、F-BDR计算与FD-OCT图像生成；
当所有波长的光子包经历步骤S3.1到S3.5，记录下的光子包的剩余能量按照波长和光程进行二维排列，得到全光谱下的F-BDR图：, />，/>为光谱最小波长，/>为光谱最大波长；根据该F-BDR图使用信号计算方法重建FD-OCT图像，并将其作为F-BDR图的另一表现形式。
其中所述的信号计算方法可采用《Mao J, Ling Y, Xue P, et al. MonteCarlo-based full-wavelength simulator of Fourier-domain optical coherencetomography[J]. Biomedical Optics Express, 2022, 13(12): 6317-6334.》中的技术。
其中根据该F-BDR图使用信号计算方法重建FD-OCT图像具体方法为，使用F-BDR值作为强度，光程作为相位，构建干涉条纹并进行傅里叶反变换重建FD-OCT图像。
进一步地，非线性二乘拟合方法可以更换为任一数据拟合方法。进一步地，为了验证本发明的有效性，使用本发明所述算法与《Mao J, Ling Y, Xue P, et al. MonteCarlo-based full-wavelength simulator of Fourier-domain optical coherencetomography[J]. Biomedical Optics Express, 2022, 13(12): 6317-6334.》提出的F-BDR的计算复杂度极高的全光谱模拟算法进行对比，并仿照上述文献计算FD-OCT信号作为进一步对比验证。按照上述文献的设置，本发明选用聚苯乙烯小球的水溶液作为仿体进行模拟，其中小球的折射率为，水的折射率为/>，聚苯乙烯小球的吸收系数为/>。共设计两种仿体：仿体1，体积分数为0.34%的2/> 半径小球的水溶液；仿体2，体积分数为2.83%的3/> 半径小球的水溶液。本发明提出的F-BDR的快速模拟方法称为TT-IS方法，模拟结果将对比三种方法：《Mao J, Ling Y, Xue P, et al. Monte Carlo-based full-wavelength simulator of Fourier-domain optical coherencetomography[J]. Biomedical Optics Express, 2022, 13(12): 6317-6334.》中提出的原始方法，其使用精确的波长依赖的米氏SPF，称之为Mie方法，该方法的结果为F-BDR的地面真值（Ground truth）；第二种方法是将上述文献中的米氏SPF替换为本发明使用的TT SPF，但散射过程不经过偏置，称之为TT方法；第三种方法是复现了《Zhao S. Advanced MonteCarlo simulation and machine learning for frequency domain optical coherencetomography[D]. California Institute of Technology, 2016.》中基于无波长依赖性的Henyey-Greenstein SPF的快速模拟方法，称之为HG-IS方法。全部模拟设置如下表1所示：
表1模拟设置参数表
共进行三种模拟实验进行验证：（1）单条A线的F-BDR计算；（2）面扫描F-BDR以及面扫描FD-OCT图像的模拟；（3）体扫描FD-OCT图像的模拟。仿体图示中颜色从深至浅依次是仿体2、仿体1和背景水。MC模拟的计算时间以中央处理器（CPU）核时计算。
在模拟实验（1）中，使用仿体1和2堆叠而成的多层仿体，对于Mie和TT方法，每波长使用2千万光子进行模拟，在本发明所述的TT-IS方法中，由于后向散射光子的采集效率提升，每波长仅使用1万光子进行模拟。由于HG-IS方法不考虑波长依赖性，故该实验对HG-IS方法不适用。计算最终模拟结果的误差为待比较方法相对于参考方法的rMSE（Relativemean square error，相对均方误差），该误差值计算表达式如下，；其中/>为待比较方法得到的F-BDR图，/>为参考方法得到的F-BDR图，和为F-BDR图的格点数目，由于傅里叶变换的性质，本发明中。在模拟实验（2）中，使用仿体1和2构造而成的复杂仿体。TT-IS方法在横向范围[-1 mm, 1 mm]范围内计算512条线的F-BDR图，并基于此得到全光谱模拟的面扫描FD-OCT图像。HG-IS方法基于单波长参数同样计算512条线，但无法给出F-BDR图，仅能给出面扫描FD-OCT图像。上述两方法的每条线的模拟中，每波长加载5千光子。由于Mie和TT方法计算耗时极大，故本次实验不适用两种方法。
在模拟实验（3）中，使用仿体1和2构造而成的复杂三维仿体。除基于本发明所述的TT-IS方法外的几种方法均无法有效扩展到三维仿体，故本实验仅适用于TT-IS方法。在模拟实验（2）的基础上，计算128个面的FD-OCT图像并合并，共计模拟512×128条线。由于此时F-BDR为四维无法展示，故只展示FD-OCT图像。每条线的每个波长加载2500光子。一维A线的F-BDR模拟结果：
如图2所示，图2(a)给出了多层仿体的图示，(b)给出了Mie方法的地面真值模拟结果，(c)和(d)分别给出了TT和TT-IS方法的模拟结果。可以看到不同层的界面清晰可见，提取的不同深度和波长上的F-BDR变化趋势也都大致相同。表2汇总了三种方法的模拟时长和互相之间的rMSE。可以看到在取得几乎相同的效果时，本发明提出的TT-IS方法实现了两个数量级以上的加速，具体来说，在与地面真值Mie方法结果仅有不到2% rMSE的情况下，实现了373倍的计算加速。
表2三种方法的模拟时长和互相之间的rMSE表
二维面扫描的F-BDR与FD-OCT图像模拟结果：TT-IS方法对应的二维面扫描的F-BDR计算结果如图3所示。图3(a)给出了二维复杂仿体的结构，图3(b)给出了模拟的F-BDR体，图3(c)给出了复杂仿体中心圆部分顶层的F-BDR图示。该F-BDR的模拟仅需204.8 核时，平均每条A线仅需0.4核时。可以看到在F-BDR上，仿体结构清晰可见，同时通过(c)中缩放展示的圆顶层F-BDR可以清晰的看到波长依赖的趋势。该F-BDR对应的FD-OCT图像如图4(b)所示，图4(a)给出了HG-IS 方法计算的FD-OCT图像，可以看到相比于HG-IS方法的结果，本发明所述快速MC模拟算法的结果更具有真实性，更清晰的观察到仿体结构。三维体扫描FD-OCT图像模拟结果：
如图5所示，图5(a)给出了三维仿体的结构模型，图5(b)给出了本发明所述TT-IS方法的模拟结构。四维的F-BDR无法被展示。可以看到仿体的三维结构清晰可见，由于上层遮挡导致的阴影也可以被准确模拟。由于巨大的计算负担，基于全光谱的三维FD-OCT模拟在之前的方法，例如Mie方法中是不可实现的，而本方法仅需5898核时，这一时间在可接受范围内。
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型。

Claims (9)

1.一种全光谱后向散射漫反射率的模拟方法，其特征在于：包括，
S1、根据TT SPF 理论推导偏置散射函数TT-B SPF；
S1.1、定义光子的偏置散射条件；
在光子包根据偏置散射函数TT-B SPF进行偏置散射的假设下，定义偏置散射条件为：偏置散射后的光子行进方向为，探测器方向为/>，则偏置散射夹角/>的余弦值，即/>，偏置散射角/>根据TT-B SPF采样获得，具体采样方式为对TT-B SPF实行拒绝-接受采样方法以得到角度值[Casella G,Robert C P, Wells M T. Generalized accept-reject sampling schemes[J]. LectureNotes-Monograph Series, 2004: 342-347.];
S1.2、根据TT SPF和偏置散射条件推导TT-B SPF的偏置单向归一化系数；
TT SPF由两个分别由前向散射和后向散射主导的单向 SPF线性组合而成，写为，；其中/>是TT SPF,/>是TT SPF的参数集合向量，/>表示前向散射SPF，/>表示后向散射SPF，/>是散射角度，满足，/>和为/>为两个方向的各向异性因子，分别控制了散射向前向和后向偏向的程度，/>和/>为增强因子，分别对/>和/>控制的效应进行增强，/>是权重因子，对前向散射和后向散射进行平衡，并满足/>；对于某一方向散射主导的单向 SPF，函数表达式为，/> ；其中/>分别为各向异性因子和增强因子，x=f代指前向，x=b代指后向，是单向归一化系数，表示为，/>根据偏置散射条件将以上归一化系数进行变形，推导出TT-B SPF中偏置单向归一化系数应写作，/>；其中/>，为偏置各向异性因子，/>为偏置增强因子，两者取值范围同/>与/>, x=f代指前向，x=b代指后向;
S1.3、根据偏置单向归一化系数推导TT-B SPF;将偏置单向归一化系数替换/>定义式中的单向归一化系数，得到偏置单向 SPF表达式为，/> ，，其中/>为偏置单向 SPF的普遍表示，x=f代指前向，x=b代指后向，为偏置散射夹角;随后线性组合前向散射和后向散射的单向 SPF，得到TT-BSPF的表达式为，/> ；其中/>为TT-B SPF，为偏置函数的权重因子，/>是TT-B SPF的参数集合向量；
S2、求解TT SPF中的波长依赖的参数集合；
S3、部署基于TT SPF、TT-B SPF和IS 的模拟方法以计算得到F-BDR；
其中步骤S3中部署模拟方法以得到F-BDR的方法为，
S3.1、加载对应当前波长的光子包并确定波长相关的散射系数；
S3.2、光子移动、吸收与边界判断；
S3.3、基于IS的光子散射过程；
S3.4、判断光子是否存活；
S3.5、判断光子是否被探测；
S3.6、F-BDR计算与FD-OCT图像生成。
2.根据权利要求1所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：所述步骤S2中求解TT SPF中的波长依赖的参数集合的方法为：S2.1、根据TT SPF中波长依赖的参数集合中的五个参数和当前波长处的米氏SPF对TT SPF进行拟合，以获取优化后的TTSPF中波长依赖的参数集合中的五个参数；
使用非线性最小二乘法建立优化方程，通过优化调整TT SPF的五个参数来拟合当前波长下的米氏SPF和TT SPF，当优化方程达到最小值时，输出此时优化后的TT SPF中波长依赖的参数集合中的五个参数；
具体方法为：在波长为时的米氏SPF定义为/>，数学表达式为：；其中/>是散射角度，/>是粒子尺寸参数，被定义为/>，/>为背景介质的折射率，/>为散射粒子半径，/>是散射效率，/>和/>是散射振幅矩阵的分量，TT SPF在波长为/>时的参数通过参数集合向量表示为/>；使用非线性最小二乘法的优化方程表示为，；其中，/>；/>是SPF的缩放因子，/>是缩放因子中的权重函数，表示为/>；其中/>是在散射角/>时的值；通过优化调整参数集合向量/>中的五个参数，使得上述优化方程达到最小值，即此时TT SPF和当前波长的米氏SPF最为接近，使得上述优化方程达到最小值的参数集合向量为最优参数集合向量，该向量的五个参数为最优的TT SPF中波长依赖的五个参数；
S2.2、使用非线性最小二乘优化得到的优化后的TT SPF中波长依赖的五个参数构造波长依赖的TT SPF组和TT-B SPF组；
通过重复步骤S2.1对所有波长对应的TT SPF和米氏SPF进行拟合，得到每个波长下优化后的TT SPF的波长依赖的参数集合，然后将这些参数集合中的五个参数带入TT SPF的表达式构造一组波长依赖的TT SPF；之后将该组TT SPF中的五个参数重用到对应的TT-B SPF中以得到一组波长依赖的TT-B SPF，即该组TT-B SPF中对应于波长的参数集合向量。
3.根据权利要求1所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：关于加载对应当前波长的光子包并确定波长相关的散射系数的方法为，加载对应当前模拟波长的光子包/>，光子包以垂直仿体上表面的方向射入仿体；此时该光子包对应的仿体散射系数/>为/>；其中/>是仿体中粒子的数字密度，/>是粒子尺寸参数，/>为散射粒子半径，/>是散射效率。
4.根据权利要求1所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：关于光子移动、吸收与边界判断的具体操作为，
光子移动：光子包被加载后，在仿体中单次移动一随机步长，步长表示为；其中/>为光子包所在当前仿体层的吸收系数，/>为0到1之间均匀分布的随机数；光子边界判断：移动结束后，判断是否击中下一仿体单元的边界，若击中，则进行反射或透射的判断，若经过反射或透射后被接收，则跳转步骤S3.5；若并未击中边界，则当移动结束后结束，根据/>执行吸收；光子吸收：移动结束后并未击中边界，光子包的能量根据/>被仿体吸收，导致能量衰减。
5.根据权利要求1所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：关于基于IS的散射，
若光子包并未逸出仿体边界，且并未被探测器接收，则当移动结束后，进行散射事件；在光子移动过程中，三维前进方向定义为，探测器方向为/>，散射后的新前进方向为/>，对于任一波长/>下的MC模拟，所有光子都从探测器垂直射入组织，即/>；基于IS确定散射角的流程包括第一次散射事件、分裂、额外偏置散射、权重更新。
6.根据权利要求1所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：关于存活判断的方法，
若光子散射后逸出仿体边界，则光子寿命终止，加载下一光子；若光子散射后并未逸出仿体边界，且能量大于事先设定的能量阈值，则光子继续存活，并返回步骤S3.2；若光子散射后并未逸出边界，但能量小于事先设定的能量阈值，则进行俄罗斯轮盘赌，若光子赌赢，则能量增大十倍并返回步骤S3.2，若赌输，则光子寿命终止，加载下一光子。
7.根据权利要求1所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：关于探测判断的方法为：若光子在边界处被探测器接收，记录光子包的剩余能量作为对应于波长的BDR />，记录光子包的行走路径长作为光程/>，随后终止该光子包的寿命，加载下一光子包。
8.根据权利要求1所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：关于F-BDR计算与FD-OCT图像生成的方法为，
当所有波长的光子包经历步骤S3.1到S3.5，记录下的光子包的剩余能量按照波长和光程进行二维排列，得到全光谱下的F-BDR图：, />，/>为光谱最小波长，/>为光谱最大波长；
根据该F-BDR图使用信号计算方法重建FD-OCT图像，并将其作为F-BDR图的另一表现形式。
9.根据权利要求8所述的一种全光谱后向散射漫反射率的模拟方法，其特征在于：根据该F-BDR图使用信号计算方法重建FD-OCT图像具体方法为，使用F-BDR值作为强度，光程作为相位，构建干涉条纹并进行傅里叶反变换重建FD-OCT图像。
CN202310929532.2A 2023-07-27 2023-07-27 一种全光谱后向散射漫反射率的模拟方法 Active CN116644621B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310929532.2A CN116644621B (zh) 2023-07-27 2023-07-27 一种全光谱后向散射漫反射率的模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310929532.2A CN116644621B (zh) 2023-07-27 2023-07-27 一种全光谱后向散射漫反射率的模拟方法

Publications (2)

Publication Number Publication Date
CN116644621A true CN116644621A (zh) 2023-08-25
CN116644621B CN116644621B (zh) 2023-10-20

Family

ID=87623399

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310929532.2A Active CN116644621B (zh) 2023-07-27 2023-07-27 一种全光谱后向散射漫反射率的模拟方法

Country Status (1)

Country Link
CN (1) CN116644621B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109856092A (zh) * 2018-11-22 2019-06-07 上海无线电设备研究所 一种海洋荧光探测半解析蒙特卡罗模拟方法
CN109995427A (zh) * 2019-03-25 2019-07-09 西安电子科技大学 一种水下上行激光通信的蒙特卡洛仿真方法
CN114662271A (zh) * 2021-12-22 2022-06-24 中国人民解放军空军工程大学 集中式mimo雷达多目标跟踪中的波束-功率-带宽联合分配方法
CN116050230A (zh) * 2022-10-25 2023-05-02 上海交通大学 基于蒙特卡洛的fd-oct的全波长信号模拟方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109856092A (zh) * 2018-11-22 2019-06-07 上海无线电设备研究所 一种海洋荧光探测半解析蒙特卡罗模拟方法
CN109995427A (zh) * 2019-03-25 2019-07-09 西安电子科技大学 一种水下上行激光通信的蒙特卡洛仿真方法
CN114662271A (zh) * 2021-12-22 2022-06-24 中国人民解放军空军工程大学 集中式mimo雷达多目标跟踪中的波束-功率-带宽联合分配方法
CN116050230A (zh) * 2022-10-25 2023-05-02 上海交通大学 基于蒙特卡洛的fd-oct的全波长信号模拟方法

Also Published As

Publication number Publication date
CN116644621B (zh) 2023-10-20

Similar Documents

Publication Publication Date Title
Ren et al. Frequency domain optical tomography based on the equation of radiative transfer
EP1844438B1 (fr) Procede et systeme de simulation ou de synthese numerique d'images echographiques
Jarabo et al. A framework for transient rendering
CN106526674B (zh) 一种三维全波形反演能量加权梯度预处理方法
CN102779350B (zh) 一种锥束ct迭代重建算法投影矩阵构建方法
CN109661684A (zh) 用于现实交互式超声模拟的射线追踪方法
Leino et al. Perturbation Monte Carlo method for quantitative photoacoustic tomography
CN104921706B (zh) 基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法
CN107392977A (zh) 单视图切伦科夫发光断层成像重建方法
Matveyev et al. Computationally efficient model of OCT scan formation by focused beams and its usage to demonstrate a novel principle of OCT-angiography
CN116050230B (zh) 基于蒙特卡洛的fd-oct的全波长信号模拟方法
Potlov et al. Localization of inhomogeneities in diffuse optical tomography based on late arriving photons
CN116644621B (zh) 一种全光谱后向散射漫反射率的模拟方法
CN110826207A (zh) 一种任意形状牙组织的oct成像快速蒙特卡洛模型
Wang et al. Accelerating Monte Carlo simulation of light propagation in tissue mimicking turbid medium based on generative adversarial networks
CN115294300A (zh) 多分支注意力先验参数化有限投影快速荧光断层重建方法
JP6973530B2 (ja) 波数ベクトルを使った体積表面生成器
Qureshi et al. Determination of optimal number of projections and parametric sensitivity analysis of operators for parallel‐ray transmission tomography using hybrid continuous genetic algorithm
Scherleitner et al. Optical tomography imaging based on higher order Born approximation of diffuse photon density waves
Yuan et al. Denoising in Monte Carlo photon transport simulations using GPU-accelerated adaptive non-local mean filter
JP6973531B2 (ja) 可変波数ベクトルを使った体積表面生成器
Gao et al. Modeling and experimental validation of dynamical effects in Bragg coherent x-ray diffractive imaging of finite crystals
Boverman et al. Estimation and statistical bounds for three-dimensional polar shapes in diffuse optical tomography
Soonthornsaratoon Gradient-based methods for quantitative photoacoustic tomography
Seesan et al. Sample and system parameter estimation from local speckle pattern by fully numerically trained deep convolution neural network

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