CN111481168B - 一种光声内窥成像图像重建方法及系统 - Google Patents
一种光声内窥成像图像重建方法及系统 Download PDFInfo
- Publication number
- CN111481168B CN111481168B CN201910080573.2A CN201910080573A CN111481168B CN 111481168 B CN111481168 B CN 111481168B CN 201910080573 A CN201910080573 A CN 201910080573A CN 111481168 B CN111481168 B CN 111481168B
- Authority
- CN
- China
- Prior art keywords
- light absorption
- sound velocity
- absorption coefficient
- iteration
- photoacoustic
- 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
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0093—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy
- A61B5/0095—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy by applying light and detecting acoustic waves, i.e. photoacoustic measurements
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/12—Diagnosis using ultrasonic, sonic or infrasonic waves in body cavities or body tracts, e.g. by using catheters
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR 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/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Analysis (AREA)
- Animal Behavior & Ethology (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- General Health & Medical Sciences (AREA)
- Surgery (AREA)
- Biophysics (AREA)
- Molecular Biology (AREA)
- Pathology (AREA)
- Medical Informatics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Radiology & Medical Imaging (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Operations Research (AREA)
- Acoustics & Sound (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computing Systems (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明公开了一种光声内窥成像图像重建方法及系统。所述重建方法包括:根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值;获取超声探测器在所述生物腔体组织内采集的光声信号的光声信号测量值;根据所述光声信号理论值以及所述光声信号测量值构建目标函数;所述目标函数包括光吸收系数分布的非线性最小二乘函数以及声速分布的非线性最小二乘函数;根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像。采用本发明所提供的重建方法及系统能够提高光声内窥成像图像重建精确度。
Description
技术领域
本发明涉及医学成像领域,特别是涉及一种光声内窥成像图像重建方法及系统。
背景技术
生物光声内窥(Photoacoustic endoscopy,PAE)成像是一种新型的非电离式生物医学功能成像方法,它结合了超声成像的高分辨率和光学成像的高对比度的优势。生物腔体组织在短脉冲激光的照射下产生光声信号,其幅值与入射光的强度成正比,其特性由组织的光吸收特性决定,通过采用合适的算法可从超声探测器采集的声压时间序列中反演重建出腔体横截面上组织的初始声压分布或者光吸收分布图,反映腔体组织的形态结构。在此基础上,还可重建组织的光学参数(主要是光吸收系数和散射系数)分布图,反映组织的功能成分。
在PAE图像重建的过程中,为了简化问题,通常假设超声波在待测组织内的传播速度是恒定的或者是均匀分布的。但在实际应用中,超声波在不同组织中传播的速度存在较大差异,考虑到生物组织的复杂性,超声波在不同生物组织成分中的传播速度是不同的。例如冠状动脉血管包括内腔、管壁内膜/中膜、管壁外膜和粥样硬化斑块等,每一层的成分都不同,而含有病变组织的腔体的组织成份更是多样化的,所以很难准确估计超声波在不同成分组织中传播时的声速分布。因此,声速恒定的假设会导致重建PAE成像图像中存在严重的声学畸变、伪影、模糊以及目标错位等问题,进而导致PAE成像图像重建精确度低。
发明内容
本发明的目的是提供一种光声内窥成像图像重建方法及系统,以解决PAE成像图像重建精确度低的问题。
为实现上述目的,本发明提供了如下方案:
一种光声内窥成像图像重建方法,包括:
根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值;
获取超声探测器在所述生物腔体组织内采集的光声信号的光声信号测量值;
根据所述光声信号理论值以及所述光声信号测量值构建目标函数;所述目标函数包括光吸收系数分布的非线性最小二乘函数以及声速分布的非线性最小二乘函数;
根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像;所述光声内窥成像图像包括光吸收系数分布图以及声速分布图。
可选的,所述根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值,具体包括:
利用准直光源模型确定在所述生物腔体组织的边界处含有光源项的辐射传输方程的扩散近似方程;
根据所述扩散近似方程,利用时域有限差分算法确定光吸收能量理论值;
根据所述光吸收能量理论值确定所述生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值。
可选的,所述根据所述光声信号理论值以及所述光声信号测量值构建目标函数,具体包括:
根据公式 以及公式 构建目标函数;其中,F1(r,μa,k-1(r),cs,k-1(r))为给定声速的前提下求解的光吸收系数分布的非线性最小二乘函数;F2(r,μa,k-1(r),cs,k-1(r))为给定光吸收系数的前提下求解的声速分布的非线性最小二乘函数;f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2,r为生物腔体组织的横截面所在的θ-l平面极坐标系中的一点;pm(r,t)为位置r处、时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;||μa,k-1(r)||TV为TV正则化项;η为TV正则化参数;||·||为2-范数; α为同伦参数;cs,0(r)为位置r处声速的初始值。
可选的,所述根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像,具体包括:
根据所述二次近似函数确定第k次迭代后位置r处的光吸收系数;
根据所述声速分布的非线性最小二乘函数计算第k-1次迭代后得到的位置r处的声速的修正值;
根据所述修正值确定第k次迭代后位置r处的声速;
根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
可选的,所述根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像之前,还包括:
获取光吸收系数的收敛公差以及声速分布的收敛公差;
计算第k-1次迭代后位置r处的光吸收系数以及第k-1次迭代后位置r处的声速;
计算所述第k次迭代后位置r处的光吸收系数以及所述第k-1次迭代后位置r处的光吸收系数之差的光吸收系数绝对值;
计算所述第k次迭代后位置r处的声速以及所述第k-1次迭代后位置r处的声速之差的声速绝对值;
判断所述光吸收系数绝对值是否小于等于所述光吸收系数的收敛公差,得到第一判断结果;
若所述第一判断结果表示为所述光吸收系数绝对值小于等于所述光吸收系数的收敛公差,判断所述声速绝对值是否小于等于所述声速分布的收敛公差,得到第二判断结果;
若所述第二判断结果表示为所述声速绝对值小于等于所述声速分布的收敛公差,根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
一种光声内窥成像图像重建系统,包括:
光声信号理论值确定模块,用于根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值;
光声信号测量值采集模块,用于获取超声探测器在所述生物腔体组织内采集的光声信号的光声信号测量值;
目标函数构建模块,用于根据所述光声信号理论值以及所述光声信号测量值构建目标函数;所述目标函数包括光吸收系数分布的非线性最小二乘函数以及声速分布的非线性最小二乘函数;
光声内窥成像图像重建模块,用于根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像;所述光声内窥成像图像包括光吸收系数分布图以及声速分布图。
可选的,所述光声信号理论值确定模块具体包括:
扩散近似方程确定单元,用于利用准直光源模型确定在所述生物腔体组织的边界处含有光源项的辐射传输方程的扩散近似方程;
光吸收能量理论值确定单元,用于根据所述扩散近似方程,利用时域有限差分算法确定光吸收能量理论值;
光声信号理论值确定单元,用于根据所述光吸收能量理论值确定所述生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值。
可选的,所述目标函数构建模块具体包括:
目标函数构建单元,用于根据公式 以及公式 构建目标函数;其中,F1(r,μa,k-1(r),cs,k-1(r))为给定声速的前提下求解的光吸收系数分布的非线性最小二乘函数;F2(r,μa,k-1(r),cs,k-1(r))为给定光吸收系数的前提下求解的声速分布的非线性最小二乘函数;f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2,r为生物腔体组织的横截面所在的θ-l平面极坐标系中的一点;pm(r,t)为位置r处、时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;||μa,k-1(r)||TV为TV正则化项;η为TV正则化参数;||·||为2-范数; α为同伦参数;cs,0(r)为位置r处声速的初始值。
可选的,所述光声内窥成像图像重建模块具体包括:
二次近似函数确定单元,用于根据公式确定所述光吸收系数分布的非线性最小二乘函数的二次近似函数;其中,X、Y为二次近似函数QL(X,Y)中的定点;f(Y)为用定点Y表示的f(μa,k-1(r));为f(Y)的梯度;L为迭代步长;
第k次迭代后位置r处的光吸收系数确定单元,用于根据所述二次近似函数确定第k次迭代后位置r处的光吸收系数;
修正值确定单元,用于根据所述声速分布的非线性最小二乘函数计算第k-1次迭代后得到的位置r处的声速的修正值;
第k次迭代后位置r处的声速确定单元,用于根据所述修正值确定第k次迭代后位置r处的声速;
光声内窥成像图像重建单元,用于根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
可选的,还包括:
收敛公差获取单元,用于获取光吸收系数的收敛公差以及声速分布的收敛公差;
计算单元,用于计算第k-1次迭代后位置r处的光吸收系数以及第k-1次迭代后位置r处的声速;
光吸收系数绝对值计算单元,用于计算所述第k次迭代后位置r处的光吸收系数以及所述第k-1次迭代后位置r处的光吸收系数之差的光吸收系数绝对值;
声速绝对值计算单元,用于计算所述第k次迭代后位置r处的声速以及所述第k-1次迭代后位置r处的声速之差的声速绝对值;
第一判断单元,用于判断所述光吸收系数绝对值是否小于等于所述光吸收系数的收敛公差,得到第一判断结果;
第二判断单元,用于若所述第一判断结果表示为所述光吸收系数绝对值小于等于所述光吸收系数的收敛公差,判断所述声速绝对值是否小于等于所述声速分布的收敛公差,得到第二判断结果;若所述第二判断结果表示为所述声速绝对值小于等于所述声速分布的收敛公差,根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
根据本发明提供的具体实施例,本发明公开了以下技术效果:本发明提供了一种光声内窥成像图像重建方法及系统,根据光声信号的光声信号理论值以及光声信号测量值构建目标函数,根据目标函数重建光生内窥成像图像;本发明根据光声信号直接重建组织中的光吸收系数分布图以及声速分布图,相对于现有的假设声速恒定的PAE成像图像重建方法,避免了在进行PAE成像前无法准确估计组织中的声速分布的问题,从而有效避免了重建PAE成像图像中存在严重的声学畸变、伪影、模糊以及目标错位等问题,提高了PAE成像图像重建精确度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为本发明所提供的光声内窥成像图像重建方法流程图;
图2为本发明所提供的XOY平面直角坐标系示意图;
图3为本发明所提供的θ-l平面极坐标系示意图;
图4为本发明所提供的光声内窥成像图像重建系统结构图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本发明的目的是提供一种光声内窥成像图像重建方法及系统,能够提高光声内窥成像图像重建精确度。
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明中各符号为:XOY为平面直角坐标系;θ-l为平面极坐标系;θ为θ-l极坐标系的极角;l为θ-l极坐标系的极径;m为腔体横截面被等角度分割的总份数;θi为成像导管的第i个测量角度,其中,i=1,2,...,m;H(r)为光吸收能量分布的理论值;为哈密顿算子;r为腔体横截面所在的θ-l平面极坐标系中的一点;Ω为待成像组织区域;Φ(r)为时刻t位置r处的光子密度函数;Φ(r,rs,q)为Φ(r)的Laplace变换;q为Laplace变换的复频率因子;δ(r-rs)为位置r处无向点光源函数的Laplace变换;c为光在腔体组织中的传播速度;为组织表面点处的外法向量;Qs为光源强度;rs为θ-l平面极坐标系中光源所在位置;Rf为扩散传输内反射系数;n为组织相对于环境的相对折射率;μa(r)为组织中位置r处的光吸收系数;μs(r)为组织中位置r处的光散射系数;g为组织的各向异性因子;h为普朗克常数;f为入射光的频率;p(r,t)为位置r处时刻t的光声信号的理论值;cs(r)为超声波在组织中位置r处的传播速度;I(t)为入射激光脉冲的时域函数;Cp为组织的比热容;β为组织的体积膨胀温度系数;uθ为组织质点在θ方向的振动速度;ul为组织质点在l方向的振动速度;ρ0为组织密度;μa,0(r)为位置r处的光吸收系数的初始值;cs,0(r)为位置r处的声速的初始值;为光吸收系数的收敛公差;为声速分布的收敛公差;k为迭代次数;pm(r,t)为位置r处时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;||·||TV为TV正则化项;||μa,k-1(r)||TV为第k-1次迭代后得到的位置r处的光吸收系数的TV正则化项;η为TV正则化参数;||·||为2-范数;QL(X,Y)为二次近似函数;L为迭代步长;X、Y为QL(X,Y)中的定点;为f(Y)的梯度;为X-Y和的内积;pL(Y)为QL(X,Y)在定点Y处的极小值点;μa,k(r)为第k次迭代后得到的位置r处的光吸收系数;α为同伦参数;g1(cs(r))为位置r处声速的目标映射;g2(cs(r))为位置r处声速的平凡映射;d(cs,k-1(r))为第k次迭代中cs,k-1(r)的修正值;B(cs,k-1(r))为g1(cs,k-1(r))的雅可比矩阵;BT(cs,k-1(r))为B(cs,k-1(r))的转置矩阵;I为单位矩阵;cs,k(r)为第k次迭代后得到的位置r处的声速;λ为迭代步长;ε1为第k次和第k-1次迭代所得的光吸收系数的绝对差;ε2为第k次和第k-1次迭代所得的声速的绝对差。
图1为本发明所提供的光声内窥成像图像重建方法流程图,如图1所示,一种光声内窥成像图像重建方法,包括:
步骤101:根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值。
采用准直光源模型,将光源视为组织边界处的内向光子流,进而嵌入到Robin边界条件中,则在组织边界处含有光源项的辐射传输方程的扩散近似方程为:
其中,是哈密顿算子;如图2-图3所示,r是腔体横截面所在的θ-l平面极坐标系中的一点;Ω是待成像组织区域;Φ(r)是时刻t位置r处的光子密度函数;Φ(r,rs,q)是Φ(r)的Laplace变换;δ(r-rs)是位置r处无向点光源函数的Laplace变换;q是Laplace变换的复频率因子;c是光在腔体组织中的传播速度;是组织表面点处的外法向量;Qs是光源强度;rs是θ-l平面极坐标系中光源所在位置;Rf是扩散传输内反射系数;n是组织相对于环境的相对折射率;μa(r)和μs(r)分别是组织中位置r处的光吸收系数和光散射系数;g是组织的各向异性因子。
采用时域有限差分(Finite difference time domain,FDTD)算法求解式(1)得到光子密度函数Φ(r),进而得到位置r处的光吸收能量理论值H(r):
H(r)=μa(r)·h·f·Φ(r) (2)
其中,h是普朗克常数,f是入射光的频率。
根据H(r),采用FDTD算法求解如下方程,得到组织产生的光声信号的理论值:
其中,p(r,t)是位置r处、时刻t的光声信号的理论值;cs(r)是超声波在位置r处的传播速度;I(t)是入射激光脉冲的时域函数;Cp是组织的比热容;β是组织的体积膨胀温度系数;uθ和ul分别是组织质点在θ方向和l方向的振动速度;ρ0是组织密度;θ是平面θ-l极坐标系的极角;l是平面θ-l极坐标系的极径。
步骤102:获取超声探测器在所述生物腔体组织内采集的光声信号的光声信号测量值。
步骤103:根据所述光声信号理论值以及所述光声信号测量值构建目标函数;所述目标函数包括光吸收系数分布的非线性最小二乘函数以及声速分布的非线性最小二乘函数。
步骤104:根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像;所述光声内窥成像图像包括光吸收系数分布图以及声速分布图。
初始化参数:
计算第k次迭代后得到的位置r处的光吸收系数:
在给定声速的前提下求解光吸收系数分布的问题表述为如下的非线性最小二乘问题:
其中,
f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2 (5)
pm(r,t)是位置r处、时刻t的光声信号的测量值;μa,k-1(r)是第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)是第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))是光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,根据式(3)前向仿真出的光声信号的理论值;||μa,k-1(r)||TV是TV正则化项;η为TV正则化参数;||·||是2-范数。
式(4)的二次近似函数为:
忽略QL(X,Y)中的常数项,则其在定点Y处的极小值点表示为:
pL(Y)是QL(X,Y)在定点Y处的极小值点。
最终得到第k次迭代后位置r处的光吸收系数:
μa,k(r)=pL(μa,k-1(r)) (8)
其中,μa,k(r)是第k次迭代后得到的位置r处的光吸收系数。
计算第k次迭代后得到的位置r处的声速:
在给定光吸收系数的前提下求解声速分布的问题表述为如下的非线性最小二乘问题:
其中,α是同伦参数;g1(cs(r))是目标映射:
g2(cs(r))是平凡映射:
其中,cs,0(r)是位置r处声速的初始值。
计算第k次迭代中cs,k-1(r)的修正值d(cs,k-1(r)):
其中,d(cs,k-1(r))是第k次迭代中cs,k-1(r)的修正值;B(cs,k-1(r))是g1(cs,k-1(r))的雅可比矩阵;BT(cs,k-1(r))是B(cs,k-1(r))的转置矩阵;I是单位矩阵。
最终得到第k次迭代后位置r处的声速:
cs,k(r)=cs,k-1(r)+λd(cs,k-1(r)) (13)
其中,cs,k(r)是第k次迭代后得到的位置r处的声速;λ是迭代步长。
获取光吸收系数的收敛公差以及声速分布的收敛公差;
计算第k-1次迭代后位置r处的光吸收系数以及第k-1次迭代后位置r处的声速;
计算所述第k次迭代后位置r处的光吸收系数以及所述第k-1次迭代后位置r处的光吸收系数之差的光吸收系数绝对值ε1=|μa,k(r)-μa,k-1(r)|;
计算所述第k次迭代后位置r处的声速以及所述第k-1次迭代后位置r处的声速之差的声速绝对值ε2=|cs,k(r)-cs,k-1(r)|;
图4为本发明所提供的光声内窥成像图像重建系统结构图,如图4所示,一种光声内窥成像图像重建系统,包括:
光声信号理论值确定模块401,用于根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值。
所述光声信号理论值确定模块401具体包括:扩散近似方程确定单元,用于利用准直光源模型确定在所述生物腔体组织的边界处含有光源项的辐射传输方程的扩散近似方程;光吸收能量理论值确定单元,用于根据所述扩散近似方程,利用时域有限差分算法确定光吸收能量理论值;光声信号理论值确定单元,用于根据所述光吸收能量理论值确定所述生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值。
光声信号测量值采集模块402,用于获取超声探测器在所述生物腔体组织内采集的光声信号的光声信号测量值。
目标函数构建模块403,用于根据所述光声信号理论值以及所述光声信号测量值构建目标函数;所述目标函数包括光吸收系数分布的非线性最小二乘函数以及声速分布的非线性最小二乘函数。
所述目标函数构建模块403具体包括:目标函数构建单元,用于根据公式 以及公式 构建目标函数;其中,F1(r,μa,k-1(r),cs,k-1(r))为给定声速的前提下求解的光吸收系数分布的非线性最小二乘函数;F2(r,μa,k-1(r),cs,k-1(r))为给定光吸收系数的前提下求解的声速分布的非线性最小二乘函数;f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2,r为生物腔体组织的横截面所在的θ-l平面极坐标系中的一点;pm(r,t)为位置r处、时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;||μa,k-1(r)||TV为TV正则化项;η为TV正则化参数;||·||为2-范数; α为同伦参数;cs,0(r)为位置r处声速的初始值。
光声内窥成像图像重建模块404,用于根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像;所述光声内窥成像图像包括光吸收系数分布图以及声速分布图。
所述光声内窥成像图像重建模块404具体包括:二次近似函数确定单元,用于根据公式确定所述光吸收系数分布的非线性最小二乘函数的二次近似函数;其中,X、Y为二次近似函数QL(X,Y)中的定点;f(Y)为用定点Y表示的f(μa,k-1(r));为f(Y)的梯度;L为迭代步长;第k次迭代后位置r处的光吸收系数确定单元,用于根据所述二次近似函数确定第k次迭代后位置r处的光吸收系数;修正值确定单元,用于根据所述声速分布的非线性最小二乘函数计算所述第k-1次迭代后得到的位置r处的声速的修正值;第k次迭代后位置r处的声速确定单元,用于根据所述修正值确定第k次迭代后位置r处的声速;光声内窥成像图像重建单元,用于根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
本发明所提供的光声内窥成像图像重建系统还包括:收敛公差获取单元,用于获取光吸收系数的收敛公差以及声速分布的收敛公差;计算单元,用于计算第k-1次迭代后位置r处的光吸收系数以及第k-1次迭代后位置r处的声速;光吸收系数绝对值计算单元,用于计算所述第k次迭代后位置r处的光吸收系数以及所述第k-1次迭代后位置r处的光吸收系数之差的光吸收系数绝对值;声速绝对值计算单元,用于计算所述第k次迭代后位置r处的声速以及所述第k-1次迭代后位置r处的声速之差的声速绝对值;第一判断单元,用于判断所述光吸收系数绝对值是否小于等于所述光吸收系数的收敛公差,得到第一判断结果;第二判断单元,用于若所述第一判断结果表示为所述光吸收系数绝对值小于等于所述光吸收系数的收敛公差,判断所述声速绝对值是否小于等于所述声速分布的收敛公差,得到第二判断结果;若所述第二判断结果表示为所述声速绝对值小于等于所述声速分布的收敛公差,根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
本发明提供了一种光声内窥成像图像重建方法及系统,可同时重建出腔体横截面的光吸收系数分布图和声速分布图,同时反映组织的光学特性和声学特性,有效减少重建图像中由声速分布不均匀所致的失真和伪影等问题,解决了因待测组织内声速分布不均匀引起的误差,提高PAE成像图像的质量。
本说明书中各个实施例采用递进的方式描述,每个实施例重点说明的都是与其他实施例的不同之处,各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言,由于其与实施例公开的方法相对应,所以描述的比较简单,相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处。综上所述,本说明书内容不应理解为对本发明的限制。
Claims (6)
1.一种光声内窥成像图像重建方法,其特征在于,包括:
根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值;
获取超声探测器在所述生物腔体组织内采集的光声信号的光声信号测量值;
根据所述光声信号理论值以及所述光声信号测量值构建目标函数;所述目标函数包括光吸收系数分布的非线性最小二乘函数以及声速分布的非线性最小二乘函数;
根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像;所述光声内窥成像图像包括光吸收系数分布图以及声速分布图;所述根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像,具体包括:
根据公式确定所述光吸收系数分布的非线性最小二乘函数的二次近似函数;其中,X、Y为二次近似函数QL(X,Y)中的定点;f(Y)为用定点Y表示的f(μa,k-1(r)),f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2,r为生物腔体组织的横截面所在的θ-l平面极坐标系中的一点;θ为θ-l平面极坐标系的极角;l为θ-l平面极坐标系的极径;pm(r,t)为位置r处、时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;▽f(Y)为f(Y)的梯度;L为迭代步长;<X-Y,▽f(Y)>为X-Y和▽f(Y)的内积;η为TV正则化参数;
根据所述二次近似函数确定第k次迭代后位置r处的光吸收系数;
根据所述声速分布的非线性最小二乘函数计算第k-1次迭代后得到的位置r处的声速的修正值;
根据所述修正值确定第k次迭代后位置r处的声速;
根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像;
所述根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像之前,还包括:
获取光吸收系数的收敛公差以及声速分布的收敛公差;
计算第k-1次迭代后位置r处的光吸收系数以及第k-1次迭代后位置r处的声速;
计算所述第k次迭代后位置r处的光吸收系数以及所述第k-1次迭代后位置r处的光吸收系数之差的光吸收系数绝对值;
计算所述第k次迭代后位置r处的声速以及所述第k-1次迭代后位置r处的声速之差的声速绝对值;
判断所述光吸收系数绝对值是否小于等于所述光吸收系数的收敛公差,得到第一判断结果;
若所述第一判断结果表示为所述光吸收系数绝对值小于等于所述光吸收系数的收敛公差,判断所述声速绝对值是否小于等于所述声速分布的收敛公差,得到第二判断结果;
若所述第二判断结果表示为所述声速绝对值小于等于所述声速分布的收敛公差,根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
2.根据权利要求1所述的光声内窥成像图像重建方法,其特征在于,所述根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值,具体包括:
利用准直光源模型确定在所述生物腔体组织的边界处含有光源项的辐射传输方程的扩散近似方程;
根据所述扩散近似方程,利用时域有限差分算法确定光吸收能量理论值;
根据所述光吸收能量理论值确定所述生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值。
3.根据权利要求1所述的光声内窥成像图像重建方法,其特征在于,所述根据所述光声信号理论值以及所述光声信号测量值构建目标函数,具体包括:
根据公式 以及公式 构建目标函数;其中,F1(r,μa,k-1(r),cs,k-1(r))为给定声速的前提下求解的光吸收系数分布的非线性最小二乘函数;F2(r,μa,k-1(r),cs,k-1(r))为给定光吸收系数的前提下求解的声速分布的非线性最小二乘函数;f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2,r为生物腔体组织的横截面所在的θ-l平面极坐标系中的一点;θ为θ-l平面极坐标系的极角;l为θ-l平面极坐标系的极径;pm(r,t)为位置r处、时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;||μa,k-1(r)||TV为TV正则化项;η为TV正则化参数; α为同伦参数;cs,0(r)为位置r处声速的初始值。
4.一种光声内窥成像图像重建系统,其特征在于,包括:
光声信号理论值确定模块,用于根据前向仿真法确定生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值;
光声信号测量值采集模块,用于获取超声探测器在所述生物腔体组织内采集的光声信号的光声信号测量值;
目标函数构建模块,用于根据所述光声信号理论值以及所述光声信号测量值构建目标函数;所述目标函数包括光吸收系数分布的非线性最小二乘函数以及声速分布的非线性最小二乘函数;
光声内窥成像图像重建模块,用于根据所述目标函数重建所述生物腔体组织的横截面上的光声内窥成像图像;所述光声内窥成像图像包括光吸收系数分布图以及声速分布图;所述光声内窥成像图像重建模块具体包括:
二次近似函数确定单元,用于根据公式确定所述光吸收系数分布的非线性最小二乘函数的二次近似函数;其中,X、Y为二次近似函数QL(X,Y)中的定点;f(Y)为用定点Y表示的f(μa,k-1(r)),f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2,r为生物腔体组织的横截面所在的θ-l平面极坐标系中的一点;θ为θ-l平面极坐标系的极角;l为θ-l平面极坐标系的极径;pm(r,t)为位置r处、时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;为f(Y)的梯度;L为迭代步长;为X-Y和的内积;η为TV正则化参数;
第k次迭代后位置r处的光吸收系数确定单元,用于根据所述二次近似函数确定第k次迭代后位置r处的光吸收系数;
修正值确定单元,用于根据所述声速分布的非线性最小二乘函数计算第k-1次迭代后得到的位置r处的声速的修正值;
第k次迭代后位置r处的声速确定单元,用于根据所述修正值确定第k次迭代后位置r处的声速;
光声内窥成像图像重建单元,用于根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像;
还包括:
收敛公差获取单元,用于获取光吸收系数的收敛公差以及声速分布的收敛公差;
计算单元,用于计算第k-1次迭代后位置r处的光吸收系数以及第k-1次迭代后位置r处的声速;
光吸收系数绝对值计算单元,用于计算所述第k次迭代后位置r处的光吸收系数以及所述第k-1次迭代后位置r处的光吸收系数之差的光吸收系数绝对值;
声速绝对值计算单元,用于计算所述第k次迭代后位置r处的声速以及所述第k-1次迭代后位置r处的声速之差的声速绝对值;
第一判断单元,用于判断所述光吸收系数绝对值是否小于等于所述光吸收系数的收敛公差,得到第一判断结果;
第二判断单元,用于若所述第一判断结果表示为所述光吸收系数绝对值小于等于所述光吸收系数的收敛公差,判断所述声速绝对值是否小于等于所述声速分布的收敛公差,得到第二判断结果;若所述第二判断结果表示为所述声速绝对值小于等于所述声速分布的收敛公差,根据所述第k次迭代后位置r处的光吸收系数以及所述第k次迭代后位置r处的声速重建所述生物腔体组织的横截面上的光声内窥成像图像。
5.根据权利要求4所述的光声内窥成像图像重建系统,其特征在于,所述光声信号理论值确定模块具体包括:
扩散近似方程确定单元,用于利用准直光源模型确定在所述生物腔体组织的边界处含有光源项的辐射传输方程的扩散近似方程;
光吸收能量理论值确定单元,用于根据所述扩散近似方程,利用时域有限差分算法确定光吸收能量理论值;
光声信号理论值确定单元,用于根据所述光吸收能量理论值确定所述生物腔体组织在短脉冲激光照射下产生的光声信号的光声信号理论值。
6.根据权利要求4所述的光声内窥成像图像重建系统,其特征在于,所述目标函数构建模块具体包括:
目标函数构建单元,用于根据公式 以及公式 构建目标函数;其中,F1(r,μa,k-1(r),cs,k-1(r))为给定声速的前提下求解的光吸收系数分布的非线性最小二乘函数;F2(r,μa,k-1(r),cs,k-1(r))为给定光吸收系数的前提下求解的声速分布的非线性最小二乘函数;f(μa,k-1(r))=||pm(r,t)-p(r,μa,k-1(r),cs,k-1(r))||2,r为生物腔体组织的横截面所在的θ-l平面极坐标系中的一点;θ为θ-l平面极坐标系的极角;l为θ-l平面极坐标系的极径;pm(r,t)为位置r处、时刻t的光声信号的测量值;μa,k-1(r)为第k-1次迭代后得到的位置r处的光吸收系数;cs,k-1(r)为第k-1次迭代后得到的位置r处的声速;p(r,μa,k-1(r),cs,k-1(r))为光吸收系数为μa,k-1(r)、声速为cs,k-1(r)的情况下,前向仿真出的光声信号的理论值;||μa,k-1(r)||TV为TV正则化项;η为TV正则化参数; α为同伦参数;cs,0(r)为位置r处声速的初始值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910080573.2A CN111481168B (zh) | 2019-01-28 | 2019-01-28 | 一种光声内窥成像图像重建方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910080573.2A CN111481168B (zh) | 2019-01-28 | 2019-01-28 | 一种光声内窥成像图像重建方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111481168A CN111481168A (zh) | 2020-08-04 |
CN111481168B true CN111481168B (zh) | 2023-04-07 |
Family
ID=71791252
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910080573.2A Active CN111481168B (zh) | 2019-01-28 | 2019-01-28 | 一种光声内窥成像图像重建方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111481168B (zh) |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107007259A (zh) * | 2017-03-29 | 2017-08-04 | 华北电力大学(保定) | 一种用于生物光声内窥成像的光吸收系数重建方法 |
CN108095690A (zh) * | 2017-12-17 | 2018-06-01 | 北京工业大学 | 基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2012138965A2 (en) * | 2011-04-08 | 2012-10-11 | University Of Florida Research Foundation, Inc. | Enhanced image reconstruction in photoacoustic tomography |
-
2019
- 2019-01-28 CN CN201910080573.2A patent/CN111481168B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107007259A (zh) * | 2017-03-29 | 2017-08-04 | 华北电力大学(保定) | 一种用于生物光声内窥成像的光吸收系数重建方法 |
CN108095690A (zh) * | 2017-12-17 | 2018-06-01 | 北京工业大学 | 基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法 |
Non-Patent Citations (2)
Title |
---|
Chao Huang et al..Joint Reconstruction of Absorbed Optical Energy Density and Sound Speed Distributions in Photoacoustic Computed Tomography: A Numerical Investigation.《IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING》.2016,第2卷(第2期),第136-149页. * |
Joint Reconstruction of Absorbed Optical Energy Density and Sound Speed Distributions in Photoacoustic Computed Tomography: A Numerical Investigation;Chao Huang et al.;《IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING》;20160630;第2卷(第2期);第136-149页 * |
Also Published As
Publication number | Publication date |
---|---|
CN111481168A (zh) | 2020-08-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Choi et al. | Practical photoacoustic tomography: realistic limitations and technical solutions | |
JP5837115B2 (ja) | 被検体情報取得装置 | |
US8457440B1 (en) | Method and system for background subtraction in medical optical coherence tomography system | |
Jin et al. | Fast and high-resolution three-dimensional hybrid-domain photoacoustic imaging incorporating analytical-focused transducer beam amplitude | |
EP3427645B1 (en) | Display data obtaining apparatus and display data obtaining method | |
CN107007259B (zh) | 一种用于生物光声内窥成像的光吸收系数重建方法 | |
Zheng et al. | 2-D image reconstruction of photoacoustic endoscopic imaging based on time-reversal | |
Cox et al. | Quantitative photoacoustic imaging: fitting a model of light transport to the initial pressure distribution | |
CN111956180B (zh) | 一种重建光声内窥层析图像的方法 | |
CN109924949A (zh) | 一种基于卷积神经网络的近红外光谱断层成像重建方法 | |
Lu et al. | Full-frequency correction of spatial impulse response in back-projection scheme using space-variant filtering for optoacoustic mesoscopy | |
Kretzek et al. | GPU-based 3D SAFT reconstruction including attenuation correction | |
Qi et al. | Cross-sectional photoacoustic tomography image reconstruction with a multi-curve integration model | |
Wang et al. | Accelerating image reconstruction in ultrasound transmission tomography using L-BFGS algorithm | |
Zheng et al. | Reconstruction of optical absorption coefficient distribution in intravascular photoacoustic imaging | |
CN111481168B (zh) | 一种光声内窥成像图像重建方法及系统 | |
Zhang et al. | A reconstruction algorithm for thermoacoustic tomography with compensation for acoustic speed heterogeneity | |
JP5645637B2 (ja) | 被検体情報取得装置および被検体情報取得方法 | |
JP7356504B2 (ja) | 物質の非線形バルク弾性の超音波推定 | |
Sun et al. | An iterative gradient convolutional neural network and its application in endoscopic photoacoustic image formation from incomplete acoustic measurement | |
Sun et al. | Simultaneous reconstruction of optical absorption property and speed of sound in intravascular photoacoustic tomography | |
Willemink et al. | Imaging of acoustic attenuation and speed of sound maps using photoacoustic measurements | |
CN111820868B (zh) | 一种生物光声内窥图像重建方法及系统 | |
Kruizinga et al. | Ultrasound-guided photoacoustic image reconstruction: image completion and boundary suppression | |
Sun et al. | Image reconstruction for endoscopic photoacoustic tomography including effects of detector responses |
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 |