[go: up one dir, main page]

CN110911007B - 基于成像光谱仪的生物组织光学参数重构方法 - Google Patents

基于成像光谱仪的生物组织光学参数重构方法 Download PDF

Info

Publication number
CN110911007B
CN110911007B CN201911385533.5A CN201911385533A CN110911007B CN 110911007 B CN110911007 B CN 110911007B CN 201911385533 A CN201911385533 A CN 201911385533A CN 110911007 B CN110911007 B CN 110911007B
Authority
CN
China
Prior art keywords
biological tissue
monte carlo
voxel
scattering
imaging spectrometer
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
CN201911385533.5A
Other languages
English (en)
Other versions
CN110911007A (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.)
Hangzhou Keluoma Photoelectric Technology Co ltd
Original Assignee
Hangzhou Keluoma Photoelectric 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 Hangzhou Keluoma Photoelectric Technology Co ltd filed Critical Hangzhou Keluoma Photoelectric Technology Co ltd
Priority to CN201911385533.5A priority Critical patent/CN110911007B/zh
Publication of CN110911007A publication Critical patent/CN110911007A/zh
Application granted granted Critical
Publication of CN110911007B publication Critical patent/CN110911007B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0075—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence by spectroscopy, i.e. measuring spectra, e.g. Raman spectroscopy, infrared absorption spectroscopy
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Medical Informatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Public Health (AREA)
  • Databases & Information Systems (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Pure & Applied Mathematics (AREA)
  • Pathology (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Biomedical Technology (AREA)
  • Epidemiology (AREA)
  • Veterinary Medicine (AREA)
  • Animal Behavior & Ethology (AREA)
  • Computing Systems (AREA)
  • Surgery (AREA)
  • Molecular Biology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Primary Health Care (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

本发明的目的是针对现有技术的不足，提出一种基于成像光谱仪的生物组织参数重构方法，通过本发明的装置与方法，探测生物组织每一个点的反射光谱，通过使用大量光谱数据代替大量空间探测数据，从而提高探测系统的紧凑性与探测效率，在生物组织光学参数重构中，构建三层生物组织模型，同时兼顾了生物组织的复杂性与计算效率，并使用扰动式蒙特卡洛仿真，对光与生物组织的作用进行高效仿真，最后采用雅克比矩阵重构算法，结合光谱探测数据以及扰动式蒙特卡洛仿真，实现生物组织光学参数的重构。

Description

基于成像光谱仪的生物组织光学参数重构方法
技术领域
本发明属于光学技术领域，具体涉及一种基于成像光谱仪的生物组织参数重构方法。
背景技术
近年来，生物医学光子学的研究备受关注，其中一个原因是光学技术在医疗中的应用是无损的、数字量化的，可为医护人员提供明确的检查结果。在生物医学光子学的研究中，生物组织的光学参数重构尤为重要，通过组织光学参数的研究，可对生物组织内部的不同物质的含量进行量化分析。然而，由于生物组织的复杂散射特性，生物组织光学参数的探测系统需要探测大量的数据，才能准确的计算生物组织参数，这对于探测系统的探测模式提出了很高的要求，往往需要布置大量的探测点开展探测工作，以获取详尽数据。
本发明专利中，我们提出一种基于成像光谱仪的生物组织参数重构方法，通过使用成像光谱仪对生物组织进行光谱图像检测，通过不同波段通道的光谱数据，实现探测区域的大量探测数据获取，并结合生物物质中重要的吸收、散射物质的光学特性，建立计算模型，重构生物组织光学参数。
发明内容
本发明的目的是针对现有技术的不足，提出一种基于成像光谱仪的生物组织参数重构方法。
本发明解决其技术问题所采用的技术方案如下：
本发明的算法部分包括三层生物组织模型建模算法、蒙特卡洛正问题建模算法以及雅克比矩阵重构算法，硬件部分包括成像光谱仪、光源与直线扫描电机，光源均匀照射生物组织，成像光谱仪以反射模式检测光源与生物组织作用的后向散射光，直线扫描电机带动成像光谱仪做直线位移扫描，获得二维空间的后向散射光谱，在算法中需要对生物组织进行三层生物组织模型建模，其中第1和第2层生物组织模型定义为非均匀生物组织，并划分网格体素，每个网格体素定义各自的生物组织的光学参数，第3层生物组织模型定义为均匀生物组织，第三层组织拥有相同的生物光学参数，基于三层生物组织模型，使用扰动蒙特卡洛仿真开展正问题建模，扰动蒙特卡洛仿真由双步蒙特卡洛仿真组成，第一步蒙特卡洛仿真与经典仿真方法基本相同，但针对“被探测到的光子”，程序将记录该光子出射时的随机数种子数值，定义为“有效种子数值”，在第二步蒙特卡洛仿真中，仿真程序将读取有效种子数值完成光子随机游走的仿真，并记录后向散射光强，结合扰动蒙特卡洛仿真的数据，构建雅克比矩阵，并结合探测的后向散射光谱数据，使用逆问题算法，求解生物组织的光学参数。
所述三层生物组织模型建模，其中第1和第2层生物组织模型定义为非均匀生物组织，划分为离散的正方体体素(Voxelj,j＝1:2N，N是大于1的整数，表示单层生物组织表面离散化后的体素数目)，第3层生物组织定义为均匀生物组织，为了便于表述，将第三层生物组织定义为Voxel2N+1，基于上述设置，当我们研究组织表面的第j号体素Voxelj时，其体素中心位置定义为rj，其散射系数和吸收系数分别定义为μs(λ,rj)和μa(λ,rj)，其中μs(λ,rj)是体素Voxelj中散射光谱系数，定义体素Voxelj中的各向异性因子为g(λ,rj)，上述两个参数可转换为Mie散射方程与散射体半径crad(rj)、散射体浓度csca(rj)、散射体相对折射率cind(rj)，在生物组织中，起吸收作用的物质主要有：水、含氧血红蛋白、非含氧血红蛋白、脂类以及胶原蛋白等，因此Voxelj的吸收系数光谱可由下式决定：
其中与分别表示水、含氧血红蛋白、非含氧血红蛋白、脂肪以及胶原蛋白的摩尔吸收系数，定义第j号体素中的水、含氧血红蛋白、非含氧血红蛋白、脂肪以及胶原蛋白的浓度：cw(rj),co(rj),cd(rj),cf(rj)与cc(rj)，在研究第Voxelj的情况下，8个未知参数定义为：p(rj)＝[cw(rj),co(rj),cd(rj),cf(rj),cc(rj),crad(rj),csca(rj),cind(rj)]T，将上述8个参数代入Mie散射理论以及式(1)，可以得到散射系数与吸收系数，整个生物组织模型的未知参数数目为8×(2N+1)，在成像光谱仪扫描探测实验中，可获得第一层每一个体素漫反射的漫反射光谱，定义为D(rmi)＝[Dλ1(rmi),Dλ2(rmi),...,DλB(rmi)]T，其中rmi表示第i(i＝1:N)号探测位置，探测位置位于第一层体素的上表面中心。
所述蒙特卡洛正问题建模算法，使用扰动蒙特卡洛仿真实现光子在生物组织中传播的建模，扰动蒙特卡洛仿真的算法由双步蒙特卡洛仿真组成，第一步蒙特卡洛仿真与经典仿真方法基本相同，利用抽样方法追踪生物组织中的光子迁移，但针对“被探测到的光子”，程序将记录该光子出射时的随机数种子数值，定义为“有效种子数值”，在第二步蒙特卡洛仿真中，仿真程序将读取有效种子数值完成光子随机游走的仿真，在第二步蒙特卡洛仿真中不会追踪“无法被探测到的光子”，提高了蒙特卡洛仿真的效率，在第二步蒙特卡洛仿真中，吸收系数的数值可随意改变，在结合实验数据前，预先针对一系列离散的散射数值以及波长数值运行第一步蒙特卡洛仿真，记录不同散射数值、不同波长下的有效种子数值，在雅克比矩阵重构算法中，将根据散射系数以及工作波长的数值，选择对应的有效种子数值，使用第二步蒙特卡洛仿真进行高效率的仿真计算。
所述雅克比矩阵重构算法，核心计算满足下列关系：
在上述式子中，每一个符号均代表一个子矩阵，雅克比矩阵内的子矩阵定义为：
其中，rmi表示蒙特卡洛模拟中的探测位置，rj表示第j(j＝1:2N+1)号体素Voxelj的中心位置，p表示当前整个生物组织的未知参数的数值，可定义为p＝[p(r1),p(r2),...,p(r2N+1)]T，子矩阵δD(rmi)定义为[δDλ1(rmi),δDλ2(rmi),...,δDλB(rmi)]T，其中δDλk(rmi)表示在探测位置rmi处第k号光谱通道的探测数值Dλk(rmi)与蒙特卡洛模拟数值Dλk(rsi)的差值，在式(2)中，M表示探测点的数目，N为探测位置数目，M＝N，将每一个探测位置设置在第一层体素顶面的中心，子矩阵δp(rj)定义为：δp(rj)＝[δcw(rj),δco(rj),δcd(rj),δcf(rj),δcc(rj),δcrad(rj),δcsca(rj),δcind(rj)]T，δp(rj)表示Voxelj的参数更新值，在逆问题的反演过程中，通过求解
可获得δp＝[δp(r1),δp(r2),...,δp(r2N+1)]T，其中是正则化参数，是一个有理数，δD＝[δD(rm1),δD(rm2),...,δD(rmM)]T，完成式(4)的计算后，根据下式更新生物组织的未知参数：
p＝[p(r1),p(r2),...,p(r2N+1)]T+[δp(r1),δp(r2),...,δp(r2N+1)]T (5)，
不断重复上述过程，并根据当前生物组织参数计算雅克比矩阵的所有元素以及漫反射光谱的模拟结果，再利用式(5)计算参数更新值，进而更新生物组织参数，直至收敛为止。
所述雅克比矩阵，其子矩阵的元素定义为Jλk(rsi,cx(rj)；p)，生物组织表面的漫反射光谱的蒙特卡洛仿真数据Dλk(rsi)，雅克比的元素Jλk(rsi,cx(rj)；p)的计算式为：
根据当前的波长λk以及散射数值，选择相应的有效种子数值(已在第一步蒙特卡洛仿真中获取)，运行第二步蒙特卡洛仿真计算雅克比矩阵中的元素，在第二步蒙特卡洛仿真中，当第g号光子包经过体素Voxelj时，记录该体素内的散射事件产生次数以及传输总步长，将数值保存为STg(rj)和Leng(rj)，当g号光子包最终从rsi位置出射并被成像光谱仪探测，探测到的光子包权重为ωg，该光子包的模拟可用于更新元素Jλk(rsi,cx(rj)；p)的数值(其初值设为0)，具体的计算方法为：在体素Voxelj中针对cx(rj)引入扰动Δcx(rj)，基于参数p(rj)+Δcx(rj)，结合Mie散射方法以及式(1)得到新的散射与吸收系数和然后根据下式计算雅克比矩阵元素：
其中在式(7)中，μs(λk,rj)和μa(λk,rj)是使用参数p(rj)计算得到的散射与吸收系数，通过大量的光子包追踪，最终可以准确计算雅克比矩阵中的子矩阵所有元素。
所述成像光谱仪，其核心元件包括成像镜头、狭缝、准直透镜、光栅、聚焦透镜和面阵相机，狭缝位于成像镜头的焦面，空间中的线状物体发出的光线被成像镜头采集，成像在狭缝位置，光线经过狭缝后，由准直透镜准直，准直透镜的前焦面与狭缝位置重合，准直的光线铜鼓光栅后发生衍射，不同波长的光衍射角度不同，并通过聚焦透镜聚焦到面阵相机中，面阵相机的感光面与聚焦透镜的后焦面重合。
所述光源，包含着可见光-红外光的波长信号，光源发出的光线照亮生物组织，被生物组织后向反射后，被成像光谱仪探测。
所述直线扫描电机，是一种电机与导轨的结合，电机的转动，带动导轨上的滑块做直接移动，成像光谱仪与滑块连接，在滑块直线移动的时候，成像光谱仪也同时做直线移动，从而实现一系列直线的光谱图像采集，最后可拼接出一个二维物体的光谱图像数据。
本发明的有益效果：通过本发明的装置与方法，探测生物组织每一个点的反射光谱，通过使用大量光谱数据代替大量空间探测数据，从而提高探测系统的紧凑性与探测效率，在生物组织光学参数重构中，构建三层生物组织模型，同时兼顾了生物组织的复杂性与计算效率，并使用扰动式蒙特卡洛仿真，对光与生物组织的作用进行高效仿真，最后采用雅克比矩阵重构算法，结合光谱探测数据以及扰动式蒙特卡洛仿真，实现生物组织光学参数的重构。
附图说明
图1为基于成像光谱仪的生物组织探测系统示意图。
图2为生物组织光学参数重构方法流程图。
图3为扰动式双步蒙特卡洛算法流程图。
具体实施方式
为了使公众能够更加清楚地理解本发明的技术实质和有益效果，申请人将在下面以实施例的方式作详细说明，但是对实施例的描述均不是对本发明方案的限制，任何依据本发明构思所作出的仅仅为形式上的而非实质性的等效变换都应视为本发明的技术方案范畴。
实施例1
下面结合附图1、附图2、附图3和实施例1对本发明作进一步说明。
一种基于成像光谱仪的生物组织参数重构方法，算法部分包括三层生物组织模型建模算法、蒙特卡洛正问题建模算法以及雅克比矩阵重构算法，硬件部分包括光源、成像光谱仪与直线扫描电机，如附图1所示，光源1和光源2均匀照射生物组织，成像光谱仪3以反射模式检测光源与生物组织作用的后向散射光，直线扫描电机4带动成像光谱仪3做直线位移扫描，获得二维空间的后向散射光谱，生物组织光学参数重构方法流程如附图2所示，在算法中需要对生物组织进行三层生物组织模型建模，其中第1和第2层生物组织模型定义为非均匀生物组织，并划分网格体素，每个网格体素定义各自的生物组织的光学参数，第3层生物组织模型定义为均匀生物组织，第三层组织拥有相同的生物光学参数，基于三层生物组织模型，使用扰动蒙特卡洛仿真开展正问题建模，扰动蒙特卡洛仿真由双步蒙特卡洛仿真组成，第一步蒙特卡洛仿真与经典仿真方法基本相同，但针对“被探测到的光子”，程序将记录该光子出射时的随机数种子数值，定义为“有效种子数值”，在第二步蒙特卡洛仿真中，仿真程序将读取有效种子数值完成光子随机游走的仿真，并记录后向散射光强，结合扰动蒙特卡洛仿真的数据，构建雅克比矩阵，并结合探测的后向散射光谱数据，使用逆问题算法，求解生物组织的光学参数。
所述三层生物组织模型建模，如附图1所示，其中第1和第2层生物组织模型定义为非均匀生物组织，划分为离散的正方体体素(Voxelj,j＝1:2N，N是大于1的整数，表示单层生物组织表面离散化后的体素数目)，第3层生物组织定义为均匀生物组织，为了便于表述，将第三层生物组织定义为Voxel2N+1，基于上述设置，当我们研究组织表面的第j号体素Voxelj时，其体素中心位置定义为rj，其散射系数和吸收系数分别定义为μs(λ,rj)和μa(λ,rj)，其中μs(λ,rj)是体素Voxelj中散射光谱系数，定义体素Voxelj中的各向异性因子为g(λ,rj)，上述两个参数可转换为Mie散射方程与散射体半径crad(rj)、散射体浓度csca(rj)、散射体相对折射率cind(rj)，在生物组织中，起吸收作用的物质主要有：水、含氧血红蛋白、非含氧血红蛋白、脂类以及胶原蛋白等，因此Voxelj的吸收系数光谱可由下式决定：
其中与分别表示水、含氧血红蛋白、非含氧血红蛋白、脂肪以及胶原蛋白的摩尔吸收系数，定义第j号体素中的水、含氧血红蛋白、非含氧血红蛋白、脂肪以及胶原蛋白的浓度：cw(rj),co(rj),cd(rj),cf(rj)与cc(rj)，在研究第Voxelj的情况下，8个未知参数定义为：p(rj)＝[cw(rj),co(rj),cd(rj),cf(rj),cc(rj),crad(rj),csca(rj),cind(rj)]T，将上述8个参数代入Mie散射理论以及式(1)，可以得到散射系数与吸收系数，整个生物组织模型的未知参数数目为8×(2N+1)，在成像光谱仪扫描探测实验中，可获得第一层每一个体素漫反射的漫反射光谱，定义为D(rmi)＝[Dλ1(rmi),Dλ2(rmi),...,DλB(rmi)]T，其中rmi表示第i(i＝1:N)号探测位置，探测位置位于第一层体素的上表面中心。
所述蒙特卡洛正问题建模算法，使用扰动蒙特卡洛仿真实现光子在生物组织中传播的建模，如附图3所示，扰动蒙特卡洛仿真的算法由双步蒙特卡洛仿真组成，第一步蒙特卡洛仿真与经典仿真方法基本相同，利用抽样方法追踪生物组织中的光子迁移，但针对“被探测到的光子”，程序将记录该光子出射时的随机数种子数值，定义为“有效种子数值”，在第二步蒙特卡洛仿真中，仿真程序将读取有效种子数值完成光子随机游走的仿真，在第二步蒙特卡洛仿真中不会追踪“无法被探测到的光子”，提高了蒙特卡洛仿真的效率，在第二步蒙特卡洛仿真中，吸收系数的数值可随意改变，在结合实验数据前，预先针对一系列离散的散射数值以及波长数值运行第一步蒙特卡洛仿真，记录不同散射数值、不同波长下的有效种子数值，在雅克比矩阵重构算法中，将根据散射系数以及工作波长的数值，选择对应的有效种子数值，使用第二步蒙特卡洛仿真进行高效率的仿真计算。
所述雅克比矩阵重构算法，核心计算满足下列关系：
在上述式子中，每一个符号均代表一个子矩阵，雅克比矩阵内的子矩阵定义为：
其中，rmi表示蒙特卡洛模拟中的探测位置，rj表示第j(j＝1:2N+1)号体素Voxelj的中心位置，p表示当前整个生物组织的未知参数的数值，可定义为p＝[p(r1),p(r2),...,p(r2N+1)]T，子矩阵δD(rmi)定义为[δDλ1(rmi),δDλ2(rmi),...,δDλB(rmi)]T，其中δDλk(rmi)表示在探测位置rmi处第k号光谱通道的探测数值Dλk(rmi)与蒙特卡洛模拟数值Dλk(rsi)的差值，在式(2)中，M表示探测点的数目，N为探测位置数目，M＝N，将每一个探测位置设置在第一层体素顶面的中心，子矩阵δp(rj)定义为：δp(rj)＝[δcw(rj),δco(rj),δcd(rj),δcf(rj),δcc(rj),δcrad(rj),δcsca(rj),δcind(rj)]T，δp(rj)表示Voxelj的参数更新值，在逆问题的反演过程中，通过求解
可获得δp＝[δp(r1),δp(r2),...,δp(r2N+1)]T，其中是正则化参数，是一个有理数，δD＝[δD(rm1),δD(rm2),...,δD(rmM)]T，完成式(4)的计算后，根据下式更新生物组织的未知参数：
p＝[p(r1),p(r2),...,p(r2N+1)]T+[δp(r1),δp(r2),...,δp(r2N+1)]T (5)，
不断重复上述过程，并根据当前生物组织参数计算雅克比矩阵的所有元素以及漫反射光谱的模拟结果，再利用式(5)计算参数更新值，进而更新生物组织参数，直至收敛为止。
所述雅克比矩阵，其子矩阵的元素定义为Jλk(rsi,cx(rj)；p)，生物组织表面的漫反射光谱的蒙特卡洛仿真数据Dλk(rsi)，雅克比的元素Jλk(rsi,cx(rj)；p)的计算式为：
根据当前的波长λk以及散射数值，选择相应的有效种子数值(已在第一步蒙特卡洛仿真中获取)，运行第二步蒙特卡洛仿真计算雅克比矩阵中的元素，在第二步蒙特卡洛仿真中，当第g号光子包经过体素Voxelj时，记录该体素内的散射事件产生次数以及传输总步长，将数值保存为STg(rj)和Leng(rj)，当g号光子包最终从rsi位置出射并被成像光谱仪探测，探测到的光子包权重为ωg，该光子包的模拟可用于更新元素Jλk(rsi,cx(rj)；p)的数值(其初值设为0)，具体的计算方法为：在体素Voxelj中针对cx(rj)引入扰动Δcx(rj)，基于参数p(rj)+Δcx(rj)，结合Mie散射方法以及式(1)得到新的散射与吸收系数和然后根据下式计算雅克比矩阵元素：
其中在式(7)中，μs(λk,rj)和μa(λk,rj)是使用参数p(rj)计算得到的散射与吸收系数，通过大量的光子包追踪，最终可以准确计算雅克比矩阵中的子矩阵所有元素。
所述成像光谱仪，其核心元件包括成像镜头、狭缝、准直透镜、光栅、聚焦透镜和面阵相机，狭缝位于成像镜头的焦面，空间中的线状物体发出的光线被成像镜头采集，成像在狭缝位置，光线经过狭缝后，由准直透镜准直，准直透镜的前焦面与狭缝位置重合，准直的光线铜鼓光栅后发生衍射，不同波长的光衍射角度不同，并通过聚焦透镜聚焦到面阵相机中，面阵相机的感光面与聚焦透镜的后焦面重合。
所述光源，包含着可见光-红外光的波长信号，光源发出的光线照亮生物组织，被生物组织后向反射后，被成像光谱仪探测。
所述直线扫描电机，是一种电机与导轨的结合，电机的转动，带动导轨上的滑块做直接移动，成像光谱仪与滑块连接，在滑块直线移动的时候，成像光谱仪也同时做直线移动，从而实现一系列直线的光谱图像采集，最后可拼接出一个二维物体的光谱图像数据。

Claims (8)

1.一种基于成像光谱仪的生物组织参数重构方法，其特征在于算法部分包括三层生物组织模型建模算法、蒙特卡洛正问题建模算法以及雅克比矩阵重构算法，硬件部分包括光源、成像光谱仪与直线扫描电机，光源(1)和光源(2)均匀照射生物组织，成像光谱仪(3)以反射模式检测光源与生物组织作用的后向散射光，直线扫描电机(4)带动成像光谱仪(3)做直线位移扫描，获得二维空间的后向散射光谱，生物组织光学参数重构方法流程在算法中需要对生物组织进行三层生物组织模型建模，其中第一和第二层生物组织模型定义为非均匀生物组织，并划分网格体素，每个网格体素定义各自的生物组织的光学参数，第三层生物组织模型定义为均匀生物组织，第三层组织拥有相同的生物光学参数，基于三层生物组织模型，使用扰动蒙特卡洛仿真开展正问题建模，扰动蒙特卡洛仿真由双步蒙特卡洛仿真组成，第一步蒙特卡洛仿真与经典仿真方法基本相同，但针对“被探测到的光子”，程序将记录该光子出射时的随机数种子数值，定义为“有效种子数值”，在第二步蒙特卡洛仿真中，仿真程序将读取有效种子数值完成光子随机游走的仿真，并记录后向散射光强，结合扰动蒙特卡洛仿真的数据，构建雅克比矩阵，其子矩阵的元素定义为Jλk(rsi,cx(rj)；p)，生物组织表面的漫反射光谱的蒙特卡洛仿真数据Dλk(rsi)，雅克比的元素Jλk(rsi,cx(rj)；p)的计算式为：并结合探测的后向散射光谱数据，使用逆问题算法，求解生物组织的光学参数。
2.根据权利要求1所述的一种基于成像光谱仪的生物组织参数重构方法，其特征在于所述的三层生物组织模型建模，其中第1和第2层生物组织模型定义为非均匀生物组织，划分为离散的正方体体素(Voxelj,j＝1:2N，N是大于1的整数，表示单层生物组织表面离散化后的体素数目)，第3层生物组织定义为均匀生物组织，为了便于表述，将第三层生物组织定义为Voxel2N+1，当我们研究组织表面的第j号体素Voxelj时，其体素中心位置定义为rj，其散射系数和吸收系数分别定义为μs(λ,rj)和μa(λ,rj)，其中μs(λ,rj)是体素Voxelj中散射光谱系数，定义体素Voxelj中的各向异性因子为g(λ,rj)，上述散射系数μs(λ,rj)和各向异性因子g(λ,rj)可转换为Mie散射方程与散射体半径crad(rj)、散射体浓度csca(rj)、散射体相对折射率cind(rj)，在生物组织中，起吸收作用的物质主要有：水、含氧血红蛋白、非含氧血红蛋白、脂类以及胶原蛋白等，因此Voxelj的吸收系数光谱可由下式决定：
其中与分别表示水、含氧血红蛋白、非含氧血红蛋白、脂肪以及胶原蛋白的摩尔吸收系数，定义第j号体素中的水、含氧血红蛋白、非含氧血红蛋白、脂肪以及胶原蛋白的浓度：cw(rj),co(rj),cd(rj),cf(rj)与cc(rj)，在研究第Voxelj的情况下，8个未知参数定义为：p(rj)＝[cw(rj),co(rj),cd(rj),cf(rj),cc(rj),crad(rj),csca(rj),cind(rj)]T，将上述散射体半径crad(rj)、散射体浓度csca(rj)、散射体相对折射率cind(rj)、第j号体素中的水、含氧血红蛋白、非含氧血红蛋白、脂肪以及胶原蛋白的浓度代入Mie散射理论以及式(1)，可以得到散射系数与吸收系数，整个生物组织模型的未知参数数目为8×(2N+1)，在成像光谱仪扫描探测实验中，可获得第一层每一个体素漫反射的漫反射光谱，定义为D(rmi)＝[Dλ1(rmi),Dλ2(rmi),...,DλB(rmi)]T，其中rmi表示第i(i＝1:N)号探测位置，探测位置位于第一层体素的上表面中心。
3.根据权利要求1所述的一种基于成像光谱仪的生物组织参数重构方法，其特征在于所述的蒙特卡洛正问题建模算法，使用扰动蒙特卡洛仿真实现光子在生物组织中传播的建模，扰动蒙特卡洛仿真的算法由双步蒙特卡洛仿真组成，第一步蒙特卡洛仿真与经典仿真方法基本相同，利用抽样方法追踪生物组织中的光子迁移，但针对“被探测到的光子”，程序将记录该光子出射时的随机数种子数值，定义为“有效种子数值”，在第二步蒙特卡洛仿真中，仿真程序将读取有效种子数值完成光子随机游走的仿真，在第二步蒙特卡洛仿真中不会追踪“无法被探测到的光子”，提高了蒙特卡洛仿真的效率，在第二步蒙特卡洛仿真中，吸收系数的数值可随意改变，在结合实验数据前，预先针对一系列离散的散射数值以及波长数值运行第一步蒙特卡洛仿真，记录不同散射数值、不同波长下的有效种子数值，在雅克比矩阵重构算法中，将根据散射系数以及工作波长的数值，选择对应的有效种子数值，使用第二步蒙特卡洛仿真进行高效率的仿真计算。
4.根据权利要求2所述的一种基于成像光谱仪的生物组织参数重构方法，其特征在于所述的雅克比矩阵重构算法，核心计算满足下列关系：
其中每一个符号均代表一个子矩阵，雅克比矩阵内的子矩阵定义为：
其中，rmi表示蒙特卡洛模拟中的探测位置，rj表示第j(j＝1:2N+1)号体素Voxelj的中心位置，p表示当前整个生物组织的未知参数的数值，可定义为p＝[p(r1),p(r2),...,p(r2N+1)]T，子矩阵δD(rmi)定义为[δDλ1(rmi),δDλ2(rmi),...,δDλB(rmi)]T，其中δDλk(rmi)表示在探测位置rmi处第k号光谱通道的探测数值Dλk(rmi)与蒙特卡洛模拟数值Dλk(rsi)的差值，在式(2)中，M表示探测点的数目，N为探测位置数目，M＝N，将每一个探测位置设置在第一层体素顶面的中心，子矩阵δp(rj)定义为：δp(rj)＝[δcw(rj),δco(rj),δcd(rj),δcf(rj),δcc(rj),δcrad(rj),δcsca(rj),δcind(rj)]T，δp(rj)表示Voxelj的参数更新值，在逆问题的反演过程中，通过求解
可获得δp＝[δp(r1),δp(r2),...,δp(r2N+1)]T，其中是正则化参数，是一个有理数，δD＝[δD(rm1),δD(rm2),...,δD(rmM)]T，完成式(4)的计算后，根据下式更新生物组织的未知参数：
p＝[p(r1),p(r2),...,p(r2N+1)]T+[δp(r1),δp(r2),...,δp(r2N+1)]T (5)，
不断重复上述过程，并根据当前生物组织参数计算雅克比矩阵的所有元素以及漫反射光谱的模拟结果，再利用式(5)计算参数更新值，进而更新生物组织参数，直至收敛为止。
5.根据权利要求1所述的一种基于成像光谱仪的生物组织参数重构方法，其特征在于所述的雅克比矩阵，其子矩阵的元素定义为Jλk(rsi,cx(rj)；p)，生物组织表面的漫反射光谱的蒙特卡洛仿真数据Dλk(rsi)，雅克比的元素Jλk(rsi,cx(rj)；p)的计算式为：
根据当前的波长λk以及散射数值，选择相应的有效种子数值运行第二步蒙特卡洛仿真计算雅克比矩阵中的元素，在第二步蒙特卡洛仿真中，当第g号光子包经过体素Voxelj时，记录该体素内的散射事件产生次数以及传输总步长，将数值保存为STg(rj)和leng(rj)，当g号光子包最终从rsi位置出射并被成像光谱仪探测，探测到的光子包权重为ωg，该光子包的模拟可用于更新元素Jλk(rsi,cx(rj)；p)的数值(其初值设为0)，具体的计算方法为：在体素Voxelj中针对cx(rj)引入扰动Δcx(rj)，基于参数p(rj)+Δcx(rj)，结合Mie散射方法以及式(1)得到新的散射与吸收系数和然后根据下式计算雅克比矩阵元素：
其中在式(7)中，μs(λk,rj)和μa(λk,rj)是使用参数p(rj)计算得到的散射与吸收系数，通过大量的光子包追踪，最终可以准确计算雅克比矩阵中的子矩阵所有元素。
6.根据权利要求1所述的一种基于成像光谱仪的生物组织参数重构方法，其特征在于所述的成像光谱仪，其核心元件包括成像镜头、狭缝、准直透镜、光栅、聚焦透镜和面阵相机，狭缝位于成像镜头的焦面，空间中的线状物体发出的光线被成像镜头采集，成像在狭缝位置，光线经过狭缝后，由准直透镜准直，准直透镜的前焦面与狭缝位置重合，准直的光线铜鼓光栅后发生衍射，不同波长的光衍射角度不同，并通过聚焦透镜聚焦到面阵相机中，面阵相机的感光面与聚焦透镜的后焦面重合。
7.根据权利要求1所述的一种基于成像光谱仪的生物组织参数重构方法，其特征在于所述的光源，包含着可见光-红外光的波长信号，光源发出的光线照亮生物组织，被生物组织后向反射后，被成像光谱仪探测。
8.根据权利要求1所述的一种基于成像光谱仪的生物组织参数重构方法，其特征在于所述的直线扫描电机，是一种电机与导轨的结合，电机的转动，带动导轨上的滑块做直接移动，成像光谱仪与滑块连接，在滑块直线移动的时候，成像光谱仪也同时做直线移动，从而实现一系列直线的光谱图像采集，最后可拼接出一个二维物体的光谱图像数据。
CN201911385533.5A 2019-12-29 2019-12-29 基于成像光谱仪的生物组织光学参数重构方法 Active CN110911007B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911385533.5A CN110911007B (zh) 2019-12-29 2019-12-29 基于成像光谱仪的生物组织光学参数重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911385533.5A CN110911007B (zh) 2019-12-29 2019-12-29 基于成像光谱仪的生物组织光学参数重构方法

Publications (2)

Publication Number Publication Date
CN110911007A CN110911007A (zh) 2020-03-24
CN110911007B true CN110911007B (zh) 2023-07-25

Family

ID=69826409

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911385533.5A Active CN110911007B (zh) 2019-12-29 2019-12-29 基于成像光谱仪的生物组织光学参数重构方法

Country Status (1)

Country Link
CN (1) CN110911007B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116050230B (zh) * 2022-10-25 2023-08-22 上海交通大学 基于蒙特卡洛的fd-oct的全波长信号模拟方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103853929A (zh) * 2014-03-17 2014-06-11 东华理工大学 一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析过程及方法
CN104679946A (zh) * 2014-11-14 2015-06-03 华中科技大学 一种基于体素的微扰荧光蒙特卡罗模拟方法
CN105740573A (zh) * 2016-03-02 2016-07-06 苏州网颢信息科技有限公司 一种用于放射线剂量计算的双步蒙特卡洛模拟方法
CN105816151A (zh) * 2016-03-10 2016-08-03 天津大学 一种基于空间频域测量的均匀组织体光学参数重建方法
WO2016150935A1 (fr) * 2015-03-24 2016-09-29 Commissariat A L'energie Atomique Et Aux Energies Alternatives Procede et dispositif pour detecter des radioelements
WO2017193699A1 (zh) * 2016-03-28 2017-11-16 陈威 一种利用数学模型计算人皮肤胶原蛋白相关3个参数的方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6015969A (en) * 1996-09-16 2000-01-18 The Regents Of The University Of California Multiple-wavelength spectroscopic quantitation of light-absorbing species in scattering media
NO325061B1 (no) * 2001-03-06 2008-01-28 Photosense As Fremgangsmate og arrangement for bestemmelse av den optiske egenskap av et multisjiktvev
US6954663B2 (en) * 2003-01-07 2005-10-11 Art Advanced Research Technologies Inc. Continuous wave optical imaging assuming a scatter-law
WO2009043045A1 (en) * 2007-09-28 2009-04-02 Duke University Systems and methods for spectral analysis of a tissue mass using an instrument, an optical probe, and a monte carlo or a diffusion algorithm
US9907494B2 (en) * 2012-04-18 2018-03-06 Hutchinson Technology Incorporated NIRS device with optical wavelength and path length correction
RU2510506C2 (ru) * 2012-04-24 2014-03-27 Белорусский Государственный Университет (Бгу) Способ определения оптических и биофизических параметров биоткани
CN103356170B (zh) * 2013-05-24 2015-02-18 天津大学 用于带异质体组织光学参数重建的快速蒙特卡罗成像方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103853929A (zh) * 2014-03-17 2014-06-11 东华理工大学 一种基于蒙卡响应矩阵的低分辨率γ能谱反演解析过程及方法
CN104679946A (zh) * 2014-11-14 2015-06-03 华中科技大学 一种基于体素的微扰荧光蒙特卡罗模拟方法
WO2016150935A1 (fr) * 2015-03-24 2016-09-29 Commissariat A L'energie Atomique Et Aux Energies Alternatives Procede et dispositif pour detecter des radioelements
CN105740573A (zh) * 2016-03-02 2016-07-06 苏州网颢信息科技有限公司 一种用于放射线剂量计算的双步蒙特卡洛模拟方法
CN105816151A (zh) * 2016-03-10 2016-08-03 天津大学 一种基于空间频域测量的均匀组织体光学参数重建方法
WO2017193699A1 (zh) * 2016-03-28 2017-11-16 陈威 一种利用数学模型计算人皮肤胶原蛋白相关3个参数的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
《Monte Carlo methods for design and analysis of radiation detectors》;William L. Dunn，J. Kenneth Shultis;《Radiation Physics and Chemistry》;852-858 *

Also Published As

Publication number Publication date
CN110911007A (zh) 2020-03-24

Similar Documents

Publication Publication Date Title
US6108576A (en) Time-resolved diffusion tomographic 2D and 3D imaging in highly scattering turbid media
Chen et al. Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates
US20090069674A1 (en) Measurement apparatus
US20090005676A1 (en) Method and Device For Reconstructing a Three-Dimensional Fluorescence Optical Tomography Image by Double Measurement
CN110946553A (zh) 一种基于高光谱图像的在体组织光学参数测量装置及方法
JPH0829329A (ja) 吸収物質濃度の空間分布画像化方法
US20120049088A1 (en) Systems, methods and computer-accessible media for hyperspectral excitation-resolved fluorescence tomography
JP5305818B2 (ja) 生体情報取得装置
CN104921706A (zh) 基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法
CN105559750A (zh) 组织结构引导的复合正则化生物发光断层成像重建方法
Kao et al. Quantifying tissue optical properties of human heads in vivo using continuous-wave near-infrared spectroscopy and subject-specific three-dimensional Monte Carlo models
Qi et al. Experimental research on noninvasive reconstruction of optical parameter fields based on transient radiative transfer equation for diagnosis applications
Flynn et al. White light-informed optical properties improve ultrasound-guided fluorescence tomography of photoactive protoporphyrin IX
Wang et al. Fluorescence molecular tomography reconstruction of small targets using stacked auto-encoder neural networks
CN110911007B (zh) 基于成像光谱仪的生物组织光学参数重构方法
US6975401B2 (en) Method of calculating optical path distribution inside scattering absorber
Li et al. Short wave infrared band Spatial-Temporal-Spectral resolved sensing system and its application in bio-samples measurement
Galeano et al. Unmixing of human skin optical reflectance maps by Non-negative Matrix Factorization algorithm
Borjkhani et al. Performance assessment of high-density diffuse optical topography regarding source-detector array topology
Baranyai et al. Analysis of laser light migration in apple tissue by Monte Carlo simulation
Hani et al. Melanin determination using optimised inverse Monte Carlo for skin—light interaction
JP3730927B2 (ja) 吸収物質濃度空間分布画像化装置
JP3276947B2 (ja) 吸収物質濃度空間分布画像装置
Yuasa et al. Simulation study on penetration depth of light in the source–detector distance using a nine-layered skin tissue model in the visible wavelength range
TW201629490A (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