CN116050230B - 基于蒙特卡洛的fd-oct的全波长信号模拟方法 - Google Patents
基于蒙特卡洛的fd-oct的全波长信号模拟方法 Download PDFInfo
- Publication number
- CN116050230B CN116050230B CN202211308649.0A CN202211308649A CN116050230B CN 116050230 B CN116050230 B CN 116050230B CN 202211308649 A CN202211308649 A CN 202211308649A CN 116050230 B CN116050230 B CN 116050230B
- Authority
- CN
- China
- Prior art keywords
- wavelength
- signal
- diffuse reflectance
- simulation
- oct
- 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
- 238000004088 simulation Methods 0.000 title claims abstract description 81
- 238000000034 method Methods 0.000 title claims abstract description 35
- 238000003325 tomography Methods 0.000 title 1
- 230000003287 optical effect Effects 0.000 claims abstract description 55
- 238000004364 calculation method Methods 0.000 claims abstract description 33
- 230000001419 dependent effect Effects 0.000 claims abstract description 24
- 238000007781 pre-processing Methods 0.000 claims abstract description 18
- 230000015572 biosynthetic process Effects 0.000 claims abstract description 15
- 238000005070 sampling Methods 0.000 claims abstract description 15
- 238000003786 synthesis reaction Methods 0.000 claims abstract description 15
- 238000001228 spectrum Methods 0.000 claims abstract description 13
- 238000004458 analytical method Methods 0.000 claims abstract description 5
- 238000012014 optical coherence tomography Methods 0.000 claims description 44
- 238000009826 distribution Methods 0.000 claims description 36
- 230000003595 spectral effect Effects 0.000 claims description 23
- 238000004422 calculation algorithm Methods 0.000 claims description 13
- 239000002245 particle Substances 0.000 claims description 9
- 238000010521 absorption reaction Methods 0.000 claims description 8
- 238000010606 normalization Methods 0.000 claims description 8
- 230000005540 biological transmission Effects 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 5
- 238000003384 imaging method Methods 0.000 claims description 4
- 238000001514 detection method Methods 0.000 claims description 3
- 230000007613 environmental effect Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000013139 quantization Methods 0.000 claims description 3
- 235000013405 beer Nutrition 0.000 claims description 2
- 238000000295 emission spectrum Methods 0.000 claims description 2
- 239000011159 matrix material Substances 0.000 claims description 2
- 238000012216 screening Methods 0.000 claims description 2
- 238000001308 synthesis method Methods 0.000 claims description 2
- 238000011002 quantification Methods 0.000 claims 1
- 238000010586 diagram Methods 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 4
- 239000004793 Polystyrene Substances 0.000 description 3
- 239000007864 aqueous solution Substances 0.000 description 3
- 239000011324 bead Substances 0.000 description 3
- 229920002223 polystyrene Polymers 0.000 description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 241001325209 Nama Species 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 239000008188 pellet Substances 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 239000000243 solution Substances 0.000 description 1
- 238000010408 sweeping Methods 0.000 description 1
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
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/08—Probabilistic or stochastic CAD
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
一种基于蒙特卡洛的FD‑OCT的全波长信号模拟方法，包括：正向MC模拟阶段、信号预处理阶段、干涉条纹合成和FD‑OCT信号计算阶段。本发明运用MC引擎和米氏散射理论进行波段内全波长采样点的模拟，并结合信号预处理进行干涉条纹合成以完成FD‑OCT信号模拟、光谱OCT模拟、波长依赖的特征参数计算等任务，在最终的FD‑OCT信号中保留了全部的波长相关的信息，相比于单波长的模拟方法，可以模拟更精确的信号、计算更精确的光学参数，同时通过光谱OCT相关的时频分析可以精确捕捉随波长变化的参数谱，例如后向散射系数。
Description
技术领域
本发明涉及的是一种信号处理领域的技术，具体是一种基于蒙特卡洛的傅里叶域光学相干断层成像(FD-OCT)的全波长信号模拟方法。
背景技术
傅里叶域光学相干断层成像(FD-OCT)是一种高分辨率的非侵入性宽谱成像方法。其使用宽谱光源，在一次测量中获取生物组织全深度的结构信息。FD-OCT在例如眼科诊疗和血管内成像的生物医学领域中具有广泛的应用。蒙特卡洛(MC)模拟方法长期以来被认为是模拟光与组织相互作用的“金标准”，近些年来也被应用于FD-OCT的理论建模。然而基于MC的模拟方法始终无法得到令人满意的FD-OCT信号，并且仅仅将应用局限于简单的信号模拟，无法模拟特殊的FD-OCT应用，例如光谱OCT。究其原因，是因为传统的基于MC的FD-OCT模拟方法违背了FD-OCT理论中的宽谱特性，使用单个波长对应的参数进行模拟，并将单波长模拟的结果复制到整个波段以使用逆离散傅里叶变换(IDFT)方法进行信号计算，从而丢弃了波长依赖的相关特征，导致信号模拟的不准确以及无法模拟波长依赖的光谱OCT等应用。
现有改进技术使用体素来建模组织，结合传统的基于单波长参数的MC模拟方法获取所需信息，并将模拟结果复制到整个所用波段进行FD-OCT干涉条纹的合成。但此类技术损失了波长相关的特征，无法扩展到信号模拟之外的应用。
发明内容
本发明针对现有技术存在的上述不足，提出一种基于蒙特卡洛的FD-OCT的全波长信号模拟方法，运用MC引擎和米氏散射理论进行波段内全波长采样点的模拟，并结合信号预处理进行干涉条纹合成以完成FD-OCT信号模拟、光谱OCT模拟、波长依赖的特征参数计算等任务，在最终的FD-OCT信号中保留了全部的波长相关的信息，相比于单波长的模拟方法，可以模拟更精确的信号、计算更精确的光学参数，同时通过光谱OCT相关的时频分析可以精确捕捉随波长变化的参数谱，例如后向散射系数。
本发明是通过以下技术方案实现的：
本发明涉及一种基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，包括：正向MC模拟阶段、信号预处理阶段、干涉条纹合成和FD-OCT信号计算阶段。
所述的正向MC模拟是指：对于每一波长，根据波长依赖的随机散射行走步长判断光子包是否击中体素块，在击中时进行反射或透射的判断，并进一步判断是否接收；在未击中边界并当步长行走结束时执行吸收并进行波长依赖的随机散射角度的采样，然后进一步判断光子包是否逃脱、是否能量过低、是否通过俄罗斯轮盘赌的判断，直至所有波长的光子包全部模拟完毕。
所述的散射系数μs(λm)＝πρr2Qsca(λm，r，ns，nb)，其中：λm为当前光子包对应波长，m＝0，1，…，M-1，ρ为仿体或组织中粒子的数字密度，r为粒子的半径，Qsca为散射效率，ns为粒子的折射率，nb为背景介质的折射率，该散射系数规定了MC模拟时的散射步长s(λm)与波长相关。
所述的散射角度由米氏散射相函数采样得到，该函数具体为： 其中：θ为[0，π]之间的散射角度，粒子尺寸参数/> S1(xm，θ)和S2(xm，θ)为散射振幅矩阵的分量。在使用波长依赖的散射步长与散射角度进行模拟之后，由MC引擎得到对应于每个波长的模拟结果：波长依赖的漫反射率R(λm)和光程l(λm)。
所述的信号预处理是指：漫反射率沿光程量化，合并量化窗格内的漫反射率，量化后的漫反射率沿波长拟合，经归一化后对每一波长，查找漫反射率沿波长的分布，对每一波长，重新分配能量。
所述的干涉条纹合成是指：按照FD-OCT物臂信号的合成方式，对每个波长对应的二维漫反射率强度分布图进行积分，并求和所有波长的结果以得到包含全波长信息的物臂信号，具体为其中：波数/>S(λm)为光源的发射光谱，N(λm)是对当前波长的MC模拟中探测到的光子数量，Rm(λm)是信号预处理阶段得到的，由当前波长的漫反射率R(λm)扩展到全波段而来的随波长变化的二维漫反射率图，/>为Rm(λm)沿光程轴的第n行，随后将物臂信号与参考臂信号r(km)进行干涉，得到FD-OCT信号计算所需的干涉条纹ID(km)＝|Is(km)+r(km)|2-|Is(km)-r(km)|2。
所述的FD-OCT信号计算，包括：
①对干涉条纹进行k域线性矫正之后进行IDFT，计算深度分辨的全波长FD-OCT信号或
②对干涉条纹进行k域线性矫正之后进行STFT的时频分析，计算仿体对应FD-OCT模拟系统的光谱信息。
优选地，根据经验证成熟的散射信号模型对该信号进行拟合，并计算所需衰减系数和波长依赖的后向散射系数。
本发明涉及一种实现上述方法的系统，包括：正向MC模拟单元、信号预处理单元、干涉条纹合成单元和信号计算单元，其中：正向MC模拟单元根据环境条件，在波段内每一个波长进行MC模拟，得到对应每一波长的光程和漫反射率信息；信号预处理单元根据光程和漫反射率信息，进行信号量化、合并、拟合、归一化、分布查找和能量重分配的预处理，得到每一波长扩展到全波段的二维漫反射率图和光子被探测的概率分布图；干涉条纹合成单元对二维漫反射率图进行包含幅值和相位的积分，得到物臂信号，再将物臂信号与参考臂信号干涉得到包含全波长信息的干涉条纹；信号计算单元根据干涉条纹进行IDFT操作，得到包含全波长信息的FD-OCT信号并进行参数计算，或根据干涉条纹进行STFT操作得到光谱信息并进行参数计算。
所述的环境条件包括：仿体、光源、波段等条件。
技术效果
本发明使用MC方法在全波长采样点使用米氏散射理论规定的随波长变化的散射系数和米氏SPF模拟FD-OCT过程，进而计算包含全波长信息的FD-OCT信号和波长依赖的光学参数，不仅保留了波长相关的特征，同时可以计算更加准确的光学参数，可以扩展到例如模拟光谱OCT等波长依赖的应用。相比于单波长的模拟方法，本发明可以模拟更精确的信号、计算更精确的光学参数，同时通过光谱OCT相关的时频分析可以精确捕捉随波长变化的参数谱，例如后向散射系数。
附图说明
图1为本发明流程图；
图2为正向MC模拟阶段示意图；
图中：(a)为2层的生物组织仿体中移动轨迹的示例；(b)为详细流程图；
图3为信号预处理示意图；
图4为双层仿体的FD-OCT信号模拟结果示意图；
图中：(a)为双层介质的示意图，(b)为该双层仿体的拟合后的漫反射率强度分布图，(c)为使用本实施例所计算的FD-OCT信号和使用传统单波长算法计算的FD-OCT信号；
图5为光谱OCT仿真及光学参数计算结果示意图；
图中：(a)(b)为使用传统单波长算法和本实施例计算的单散射信号进行衰减系数计算的结果，(c)(d)为使用传统单波长算法和本实施例合成的干涉条纹进行光谱OCT的归一化光谱计算的结果，(e)为将传统单波长算法和本实施例的光谱计算的归一化μb(λ)和米氏散射理论规定的真值进行对比的结果。
具体实施方式
如图1所示，为本实施例涉及一种基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，包括：正向MC模拟阶段、信号预处理阶段、干涉条纹合成和信号计算阶段。
如图2所示，所述的正向MC模拟阶段，具体为：采用基于单波长的经典算法(Hartinger,Nam,Chico-Calero等人《Monte Carlo modeling ofangiographic opticalcoherence tomography》[J].Biomedical optics express,2014,5(12):4338-4349.)，将生物组织仿体以体素块建模，并对每一波长迁移其MC模拟步骤。如图可见，该阶段首先加载对应不同波长λm的光子包P(λm)，并根据米氏理论推导波长相关的散射系数μs(λm)，根据比尔定律得到光子随机的散射行走步长其中：μa为光子所在当前仿体层的吸收系数，ξ为0到1之间均匀分布的随机数。随后根据该步长判断是否击中下一体素块，若击中，则进行反射或透射的判断，若经过反射或透射后被接收，则记录光子包的剩余能量作为漫反射率R(λm)，记录行走路径长作为光程l(λm)；若并未击中边界，则当步长行走结束，根据μa执行吸收。随后根据pMie(θ；xm)进行散射角度θ(λm)的采样，本实施例采取拒绝-接受采样方法进行采样(Lux I,Koblinger L.Monte Carlo particle transportmethods:neutronand photon calculations[M].CRC press,2018.)。散射过后，对光子包进行是否逃脱、是否能量过低、是否通过俄罗斯轮盘赌等判断。当所有波长的光子包全部模拟完毕，将得到对应于每个波长的模拟结果：波长依赖的漫反射率R(λm)和光程l(λm)。
如图3所示，所述的信号预处理具体包括：
第一步，漫反射率沿光程量化：将FD-OCT模拟的理论成像光程2π/δk划分为M个窗格，每个窗格长度为其中δk是在k域上的线性采样间隔，不失一般性，以窗格的上界为光程的离散格点[0，Δl，2Δl，…，(M-1)Δl]。对应于λ0，…，λm，…，λM-1的沿光程随机连续分布的漫反射率R(λ0)，…，R(λm)，…，R(λM-1)将被上述M个窗格量化。
第二步，合并量化窗格内的漫反射率：将光程的量化窗格内的漫反射率相加，视为在离散的光程格点[0，Δl，2Δl，…，L]上的漫反射率大小，得到量化后的以离散光程和波长为正交坐标的漫反射率强度分布图。
第三步，量化后的漫反射率沿波长拟合：使用四阶多项式，沿波长拟合在离散格点上的漫反射率，得到拟合后的以离散光程和波长为正交坐标的漫反射率强度分布图。
第四步，归一化处理：将上一步骤的漫反射率强度分布图除以全图求和的总能量进行归一化，得到光子以(离散光程，波长)为坐标的被探测概率分布图Φ(l，λ)。
第五步，对每一波长，查找漫反射率沿波长的分布：对于任一波长λm，漫反射率分布于l(λm)，可看作一条一维的漫反射率线。在被探测概率分布图Φ(l，λ)中定位、查找、提取该漫反射率线对应的沿波长变化的概率分布Φ(l(λm)，λ)。由于l(λm)精度较高，且MC的散射极具随机性，l(λm)往往无法落在离散的光程格点[0，Δl，2Δl，…，L]上，并由于Φ(l，λ)仅在离散光程格点上有值，因此需要插值出对应的概率分布。考虑采集到的某个光子的漫反射率Ri(λm)，其对应的光程为li(λm)不落在[0，Δl，2Δl，…，L]中的某个格点上，其前后两个格点分别为i-Δl和i+Δl，对应其前后格点的沿波长变化的概率分布Φ(i-Δl，λ)和Φ(i+Δl，λ)，则使用样条插值，通过两条分布曲线内插出Φ(li(λm)，λ)。
第六步，对每一波长，重新分配能量：对于任一波长λm的模拟结果，将R(λm)根据Φ(l(λm)，λ)的权重沿波长进行能量重新分配，即将对应光程l(λm)的一条一维的漫反射率线R(λm)，扩展为对应光程l(λm)和全波长的二维漫反射率强度分布图Rm(λ)。
所述的干涉条纹合成具体为：首先采用基于单波长的经典算法(Hartinger,Nam,Chico-Calero等人《Monte Carlo modeling ofangiographic optical coherencetomography》[J].Biomedical optics express,2014,5(12):4338-4349.)中对扫频光源FD-OCT的设置，将S(λm)设置为1，r(km)设置为常数，随后对每个波长对应的二维漫反射率强度分布图进行积分，并求和所有波长的结果，合成包含全波长信息的物臂信号，并与r(km)干涉生成干涉条纹ID(km)并在k域上进行线性矫正。
所述的信号计算阶段包括：
a)当应用于FD-OCT信号模拟任务，则对ID(km)进行IDFT计算得到FD-OCT信号；
b)应用于光谱OCT的模拟，则对ID(km)进行STFT计算得到光谱OCT的光谱。
以上两种应用可分别扩展计算光学参数的应用，通过数据筛选，可以删除多次散射带来的影响，得到仅有单次散射信号贡献的干涉条纹根据单散射模型，则单散射的FD-OCT信号/>其中：z是深度，μb是后向散射系数。通过拟合对/>进行IDFT之后的FD-OCT信号的对数域的斜率，得到衰减系数(μs+μa)；通过拟合对/>进行STFT之后的光谱OCT计算光谱的最大值，可以获得随波长变化的μb(λ)。由米氏散射理论可推得聚苯乙烯小球的后向散射系数/> 其中：NA为数值孔径，可由光子采集时的滤波器推得。
经过具体实际实验，本实施例选用聚苯乙烯小球的水溶液作为仿体进行模拟，其中小球的折射率为ns＝1.58，水的折射率为nb＝1.33，聚苯乙烯小球的吸收系数为μa＝0。共设计两种仿体：仿体1，体积分数为0.21％的1μm半径小球的水溶液；仿体2，体积分数为0.67％的2μm半径小球的水溶液。关于其余的模拟参数如下表所示：
对于FD-OCT信号模拟的应用，选取仿体1和仿体2堆叠的双层仿体，其中仿体1的物理厚度为0.3mm，仿体2的厚度为0.5mm，由于背景介质都是水，则折射率匹配，没有内反射产生。每个波长使用10亿光子进行模拟，1024个波长采样，总模拟光子数为10240亿。同时使用以Hartinger,Nam,Chico-Calero等人所提出的基于单波长的经典算法为代表的FD-OCT模拟引擎进行对比，选取中心波长1260nm作为模拟波长，模拟光子数为10240亿。
对于光谱OCT的仿真以及光学参数计算的应用，选取半无限大的仿体2作为模拟仿体。每个波长使用2亿光子进行模拟，1024个波长采样，总模拟光子数为2048亿。同时使用与HartingerA E,NamA S,Chico-Calero I,等人在.《Monte Carlo modelingofangiographic optical coherence tomography》([J].Biomedical optics express,2014,5(12):4338-4349.)中公开的基于单波长的经典算法为代表的FD-OCT模拟引擎进行对比，选取中心波长1260nm作为模拟波长，模拟光子数为2048亿。
如图4所示，为上述方法的双层仿体的信号模拟结果：如图可见拟合后的漫反射率强度分布图为上下两层仿体完全不同的波长依赖特征，同时准确的为两个交界面的位置0.8mm＝0.3mm×1.33×2与2.13mm＝0.8mm×1.33×2，其中因子2是由于光程是光子往返的路径长。在模拟的信号中，本实施例计算的全波长FD-OCT信号为准确的三个界面0mm，0.4mm＝0.3mm×1.33以及1.06mm＝0.8mm×1.33，而单波长方法计算的FD-OCT信号并未给出准确的界面，尤其是第二个界面：两层介质的交界面。
如图5所示，为上述方法的光谱OCT仿真及光学参数计算结果：如图可见，由于在该实施例中，吸收系数为0，则衰减系数等于散射系数。(a)中单波长的模拟信号拟合出的散射系数(7.0mm-1)与使用的散射系数真值(8.3mm-1)相差较大，而(b)中本实施例的模拟信号拟合出的散射系数(8.2mm-1)落在使用的散射系数真值区间([8mm-1，8.7mm-1])，吻合较好，表明本实施例的算法更加准确。(c)和(d)中展示的本实施例的光谱计算结果的趋势和拟合后的漫反射率强度分布图的趋势吻合。(e)中进一步计算了传统算法和本实施例的光谱信息中后向散射系数计算量的均方误差。本实施例的均方误差为0.11，远小于传统算法的7.67的均方误差，这表明本实施例可以扩展到计算波长依赖的光学参数的应用，而传统算法无法做到。
综上，与HartingerA E,Nam A S,Chico-Calero I,等人在.《Monte Carlomodeling of angiographic optical coherence tomography》([J].Biomedical opticsexpress,2014,5(12):4338-4349.)中公开的技术相比较，本方法更符合FD-OCT的宽谱特性，将得到更准确的模拟结果。与WangY,Bai L.在《Accurate Monte Carlo simulationoffrequency-domain optical coherence tomography》([J].International journalfor numerical methods in biomedical engineering,2019,35(4):e3177.)中公开的技术相比，本方法可以扩展到信号模拟之外的应用，例如光谱OCT的模拟，同时可以计算波长依赖的光学参数。
上述具体实施可由本领域技术人员在不背离本发明原理和宗旨的前提下以不同的方式对其进行局部调整，本发明的保护范围以权利要求书为准且不由上述具体实施所限，在其范围内的各个实现方案均受本发明之约束。
Claims (7)
1.一种基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，其特征在于，包括：正向MC模拟阶段、信号预处理阶段、干涉条纹合成和信号计算阶段；
所述的正向MC模拟是指：对于每一波长，根据波长依赖的随机散射行走步长判断光子包是否击中体素块，在击中时进行反射或透射的判断，并进一步判断是否接收；在未击中边界并当步长行走结束时执行吸收并进行波长依赖的随机散射角度的采样，然后进一步判断光子包是否逃脱、是否能量过低、是否通过俄罗斯轮盘赌的判断，直至所有波长的光子包全部模拟完毕；
所述的干涉条纹合成是指：按照FD-OCT物臂信号的合成方式，对每个波长对应的二维漫反射率强度分布图进行积分，并求和所有波长的结果以得到包含全波长信息的物臂信号；
所述的FD-OCT信号计算，包括：
①对干涉条纹进行k域线性矫正之后进行IDFT，计算深度分辨的全波长FD-OCT信号，或
②对干涉条纹进行k域线性矫正之后进行STFT的时频分析，计算仿体对应FD-OCT模拟系统的光谱信息；
所述的散射角度由米氏散射相函数采样得到，该函数具体为： 其中：θ为[0，π]之间的散射角度，粒子尺寸参数/>λm为当前光子包对应波长，m＝0，1，…，M-1，Qsca为散射效率，nb为背景介质的折射率，r为粒子的半径，S1(xm，θ)和S2(xm，θ)为散射振幅矩阵的分量，在使用波长依赖的散射步长与散射角度进行模拟之后，由MC引擎得到对应于每个波长的模拟结果：波长依赖的漫反射率R(λm)和光程l(λm)；
所述的干涉条纹合成具体为 其中：波数/>S(λm)为光源的发射光谱，N(λm)是对当前波长的MC模拟中探测到的光子数量，Rm(λm)是信号预处理阶段得到的，由当前波长的漫反射率R(λm)扩展到全波段而来的随波长变化的二维漫反射率图，/>为Rm(λm)沿光程轴的第n行，随后将物臂信号与参考臂信号r(km)进行干涉，得到FD-OCT信号计算所需的干涉条纹ID(km)＝|Is(km)+r(km)|2-|Is(km)-r(km)|2。
2.根据权利要求1所述的基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，其特征是，所述的信号预处理是指：漫反射率沿光程量化，合并量化窗格内的漫反射率，量化后的漫反射率沿波长拟合，经归一化后对每一波长，查找漫反射率沿波长的分布，对每一波长，重新分配能量。
3.根据权利要求1所述的基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，其特征是，具体包括：
所述的正向MC模拟阶段，具体为：采用基于单波长的经典算法将生物组织仿体以体素块建模，并对每一波长迁移其MC模拟步骤；该阶段首先加载对应不同波长λm的光子包P(λm)，并根据米氏理论推导波长相关的散射系数μs(λm)，根据比尔定律得到光子随机的散射行走步长其中：μa为光子所在当前仿体层的吸收系数，ξ为0到1之间均匀分布的随机数；随后根据该步长判断是否击中下一体素块，若击中，则进行反射或透射的判断，若经过反射或透射后被接收，则记录光子包的剩余能量作为漫反射率R(λm)，记录行走路径长作为光程l(λm)；若并未击中边界，则当步长行走结束，根据μa执行吸收；随后根据pMie(θ；xm)进行散射角度θ(λm)的采样，采取拒绝-接受采样方法进行采样；散射过后，对光子包进行是否逃脱、是否能量过低、是否通过俄罗斯轮盘赌判断；当所有波长的光子包全部模拟完毕，将得到对应于每个波长的模拟结果：波长依赖的漫反射率R(λm)和光程l(λm)。
4.根据权利要求1或2所述的基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，其特征是，所述的信号预处理具体包括：
第一步，漫反射率沿光程量化：将FD-OCT模拟的理论成像光程2π/δk划分为M个窗格，每个窗格长度为其中δk是在k域上的线性采样间隔，不失一般性，以窗格的上界为光程的离散格点[0，Δl，2Δl，…，(M-1)Δl]；对应于λ0，…，λm，…，λM-1的沿光程随机连续分布的漫反射率R(λ0)，…，R(λm)，…，R(λM-1)将被上述M个窗格量化；
第二步，合并量化窗格内的漫反射率：将光程的量化窗格内的漫反射率相加，视为在离散的光程格点[0，Δl，2Δl，…，L]上的漫反射率大小，得到量化后的以离散光程和波长为正交坐标的漫反射率强度分布图；
第三步，量化后的漫反射率沿波长拟合：使用四阶多项式，沿波长拟合在离散格点上的漫反射率，得到拟合后的以离散光程和波长为正交坐标的漫反射率强度分布图；
第四步，归一化处理：将上一步骤的漫反射率强度分布图除以全图求和的总能量进行归一化，得到光子以(离散光程，波长)为坐标的被探测概率分布图Φ(l，λ)；
第五步，对每一波长，查找漫反射率沿波长的分布：对于任一波长λm，漫反射率分布于l(λm)，可看作一条一维的漫反射率线；在被探测概率分布图Φ(l，λ)中定位、查找、提取该漫反射率线对应的沿波长变化的概率分布Φ(l(λm)，λ)；由于l(λm)精度较高，且MC的散射极具随机性，l(λm)往往无法落在离散的光程格点[0，Δl，2Δl，…，L]上，并由于Φ(l，λ)仅在离散光程格点上有值，因此需要插值出对应的概率分布；考虑采集到的某个光子的漫反射率Ri(λm)，其对应的光程为li(λm)不落在[0，Δl，2Δl，…，L]中的某个格点上，其前后两个格点分别为i-Δl和i+Δl，对应其前后格点的沿波长变化的概率分布Φ(i-Δl，λ)和Φ(i+Δl，λ)，则使用样条插值，通过两条分布曲线内插出Φ(li(λm)，λ)；
第六步，对每一波长，重新分配能量：对于任一波长λm的模拟结果，将R(λm)根据Φ(l(λm)，λ)的权重沿波长进行能量重新分配，即将对应光程l(λm)的一条一维的漫反射率线R(λm)，扩展为对应光程l(λm)和全波长的二维漫反射率强度分布图Rm(λ)。
5.根据权利要求3所述的基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，其特征是，所述的散射系数μs(λm)＝πρr2Qsca（λm，r，ns，nb)，其中：ρ为仿体或组织中粒子的数字密度，ns为粒子的折射率。
6.根据权利要求1所述的基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法，其特征是，所述的信号计算阶段包括：
a)当应用于FD-OCT信号模拟任务，则对ID(km)进行IDFT计算得到FD-OCT信号；
b)应用于光谱OCT的模拟，则对ID(km)进行STFT计算得到光谱OCT的光谱；
以上两种应用可分别扩展计算光学参数的应用，通过数据筛选，删除多次散射带来的影响，得到仅有单次散射信号贡献的干涉条纹根据单散射模型，则单散射的FD-OCT信号/>其中：z是深度，μb是后向散射系数；通过拟合对进行IDFT之后的FD-OCT信号的对数域的斜率，得到衰减系数(μs+μa)；通过拟合对进行STFT之后的光谱OCT计算光谱的最大值，获得随波长变化的μb(λ)；由米氏散射理论可推得仿体的后向散射系数/>其中：NA为数值孔径，可由光子采集时的滤波器推得。
7.一种实现权利要求1-6中任一所述基于蒙特卡洛的傅里叶域光学相干断层成像的全波长信号模拟方法的系统，其特征在于，包括：正向MC模拟单元、信号预处理单元、干涉条纹合成单元和信号计算单元，其中：正向MC模拟单元根据环境条件，在波段内每一个波长进行MC模拟，得到对应每一波长的光程和漫反射率信息；信号预处理单元根据光程和漫反射率信息，进行信号量化、合并、拟合、归一化、分布查找和能量重分配的预处理，得到每一波长扩展到全波段的二维漫反射率图和光子被探测的概率分布图；干涉条纹合成单元对二维漫反射率图进行包含幅值和相位的积分，得到物臂信号，再将物臂信号与参考臂信号干涉得到包含全波长信息的干涉条纹；信号计算单元根据干涉条纹进行IDFT操作，得到包含全波长信息的FD-OCT信号并进行参数计算，或根据干涉条纹进行STFT操作，得到光谱信息并进行参数计算。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211308649.0A CN116050230B (zh) | 2022-10-25 | 2022-10-25 | 基于蒙特卡洛的fd-oct的全波长信号模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211308649.0A CN116050230B (zh) | 2022-10-25 | 2022-10-25 | 基于蒙特卡洛的fd-oct的全波长信号模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116050230A CN116050230A (zh) | 2023-05-02 |
| CN116050230B true CN116050230B (zh) | 2023-08-22 |
Family
ID=86116976
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211308649.0A Active CN116050230B (zh) | 2022-10-25 | 2022-10-25 | 基于蒙特卡洛的fd-oct的全波长信号模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116050230B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116644621B (zh) * | 2023-07-27 | 2023-10-20 | 始终(无锡)医疗科技有限公司 | 一种全光谱后向散射漫反射率的模拟方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008103486A1 (en) * | 2007-02-23 | 2008-08-28 | Duke University | Scaling method for fast monte carlo simulation of diffuse reflectance spectra |
| JP2011053565A (ja) * | 2009-09-03 | 2011-03-17 | Nippon Telegr & Teleph Corp <Ntt> | 信号分析装置、信号分析方法、プログラム、及び記録媒体 |
| CN103344569A (zh) * | 2013-06-21 | 2013-10-09 | 中国科学院上海光学精密机械研究所 | 偏振复频域光学相干层析成像方法和系统 |
| CN105334172A (zh) * | 2015-10-23 | 2016-02-17 | 浙江科技学院 | 一种水果果肉组织光学特性参数的重构方法 |
| CN107328740A (zh) * | 2017-06-06 | 2017-11-07 | 中国科学院上海光学精密机械研究所 | 偏振频域光学相干层析成像光谱校准方法 |
| CN110911007A (zh) * | 2019-12-29 | 2020-03-24 | 杭州科洛码光电科技有限公司 | 基于成像光谱仪的生物组织光学参数重构方法 |
| CN111289470A (zh) * | 2020-02-06 | 2020-06-16 | 上海交通大学 | 基于计算光学的oct测量成像方法 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7391520B2 (en) * | 2005-07-01 | 2008-06-24 | Carl Zeiss Meditec, Inc. | Fourier domain optical coherence tomography employing a swept multi-wavelength laser and a multi-channel receiver |
| WO2011091369A1 (en) * | 2010-01-22 | 2011-07-28 | Duke University | Multiple window processing schemes for spectroscopic optical coherence tomography (oct) and fourier domain low coherence interferometry |
| US8921767B2 (en) * | 2010-08-02 | 2014-12-30 | The Johns Hopkins University | Automatic calibration of fourier-domain optical coherence tomography systems |
| US9192294B2 (en) * | 2012-05-10 | 2015-11-24 | Carl Zeiss Meditec, Inc. | Systems and methods for faster optical coherence tomography acquisition and processing |
| US20200284572A1 (en) * | 2019-03-06 | 2020-09-10 | Axsun Technologies, Inc. | Method for improving phase stability of phase unstable optical coherence tomography |
-
2022
- 2022-10-25 CN CN202211308649.0A patent/CN116050230B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008103486A1 (en) * | 2007-02-23 | 2008-08-28 | Duke University | Scaling method for fast monte carlo simulation of diffuse reflectance spectra |
| JP2011053565A (ja) * | 2009-09-03 | 2011-03-17 | Nippon Telegr & Teleph Corp <Ntt> | 信号分析装置、信号分析方法、プログラム、及び記録媒体 |
| CN103344569A (zh) * | 2013-06-21 | 2013-10-09 | 中国科学院上海光学精密机械研究所 | 偏振复频域光学相干层析成像方法和系统 |
| CN105334172A (zh) * | 2015-10-23 | 2016-02-17 | 浙江科技学院 | 一种水果果肉组织光学特性参数的重构方法 |
| CN107328740A (zh) * | 2017-06-06 | 2017-11-07 | 中国科学院上海光学精密机械研究所 | 偏振频域光学相干层析成像光谱校准方法 |
| CN110911007A (zh) * | 2019-12-29 | 2020-03-24 | 杭州科洛码光电科技有限公司 | 基于成像光谱仪的生物组织光学参数重构方法 |
| CN111289470A (zh) * | 2020-02-06 | 2020-06-16 | 上海交通大学 | 基于计算光学的oct测量成像方法 |
Non-Patent Citations (1)
| Title |
|---|
| A Mixed-Signal Framework for Modelling Fourier-Domain Optical Coherence Tomography;Yuye Ling, Mengyuan Wang, Yu Gan, Xinwen Yao, Leopold Schmetterer, Chuanqing Zhou, Yikai Su;Asia Communications and Photonics Conference/International Conference on Information Photonics and Optical Communications 2020 (ACP/IPOC);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116050230A (zh) | 2023-05-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN100576224C (zh) | 光在生物组织中传输特性的定量蒙特卡罗模拟方法 | |
| CN101526465A (zh) | 快速多波长组织光学参数测量装置及反构方法 | |
| CN105987881B (zh) | 光谱数据干扰抑制方法、建模方法、预测方法和处理装置 | |
| Zhang et al. | Denoising method based on CNN-LSTM and CEEMD for LDV signals from accelerometer shock testing | |
| CN116050230B (zh) | 基于蒙特卡洛的fd-oct的全波长信号模拟方法 | |
| CN105816151B (zh) | 一种基于空间频域测量的均匀组织体光学参数重建方法 | |
| JP2014055939A (ja) | 物体内部推定装置およびその方法 | |
| CN113191073B (zh) | 一种基于bp神经网络的光学特性参数反演方法 | |
| Wang et al. | Accurate Monte Carlo simulation of frequency‐domain optical coherence tomography | |
| Xu et al. | Photon migration in turbid media using a cumulant approximation to radiative transfer | |
| Abdurashitov et al. | A robust model of an OCT signal in a spectral domain | |
| Mao et al. | Monte Carlo-based full-wavelength simulator of Fourier-domain optical coherence tomography | |
| CN119620102A (zh) | 基于光子计数激光雷达的开放水体叶绿素a浓度剖面探测方法 | |
| CN108871595A (zh) | 超时间分辨冲击波速度计算方法 | |
| Premru et al. | Monte Carlo simulation of radiation transfer in human skin with geometrically correct treatment of boundaries between different tissues | |
| JP5783564B2 (ja) | 三次元光散乱体の実効散乱係数の算定方法 | |
| Frolov et al. | Monte-Carlo simulation of OCT structural images of human skin using experimental B-scans and voxel based approach to optical properties distribution | |
| Lin et al. | High‐accuracy optical coherence elastography digital volume correlation methods to measure depth regions with low correlation | |
| Subramanian et al. | Penetration depth of low-coherence enhanced backscattered light in subdiffusion regime | |
| Alexander et al. | Light propagation in turbid tissue-like scattering media | |
| Golabchi et al. | Refractive effects on optical measurement of alveolar volume: a 2-D ray-tracing approach | |
| Graber et al. | Evaluation of steady-state, time-and frequency-domain data for the problem of optical diffusion tomography | |
| CN118941658A (zh) | 一种光束比拟合法提取共聚焦曲面进行衰减系数成像方法 | |
| CN116644621B (zh) | 一种全光谱后向散射漫反射率的模拟方法 | |
| Watté et al. | Monte Carlo modeling of light transfer in food |
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 |