CN113724351A - 一种光声图像衰减校正方法 - Google Patents

一种光声图像衰减校正方法 Download PDF

Info

Publication number
CN113724351A
CN113724351A CN202110974322.6A CN202110974322A CN113724351A CN 113724351 A CN113724351 A CN 113724351A CN 202110974322 A CN202110974322 A CN 202110974322A CN 113724351 A CN113724351 A CN 113724351A
Authority
CN
China
Prior art keywords
photoacoustic image
absorption coefficient
tissue absorption
phi
luminous flux
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
CN202110974322.6A
Other languages
English (en)
Other versions
CN113724351B (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.)
Southern Medical University
Original Assignee
Southern Medical University
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 Southern Medical University filed Critical Southern Medical University
Priority to CN202110974322.6A priority Critical patent/CN113724351B/zh
Publication of CN113724351A publication Critical patent/CN113724351A/zh
Application granted granted Critical
Publication of CN113724351B publication Critical patent/CN113724351B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

本发明的一种光声图像衰减校正方法，包括如下步骤：S1：将重建得到的光声图像P转为列向量形式p，将p代入线性系统方程p＝φμa中；S2：采用最小二乘方式将线性系统方程构造为优化问题模型
Figure DDA0003226764460000011
S3：通过分步迭代梯度下降方法对优化问题模型进行迭代求解，获得组织吸收系数μa；S4：将步骤S3获得的组织吸收系数μa变换为矩阵形式，并生成光通量分布矩阵Φ，将光声图像P和光通量分布矩阵Φ代入
Figure DDA0003226764460000012
得到校正后的光声图像P′。该方法采用逐点迭代的方式进行，不需要人工介入就能够实现光通量分布的校正，还能够提高深处光声成像能力，提高了光声成像的质量，且校正速度快、能够直接集成到各光声成像系统中。

Description

一种光声图像衰减校正方法
技术领域
本发明涉及生物医学成像领域，特别是涉及一种光声图像衰减校正方法。
背景技术
光声断层成像(PAT)是一种新兴的生物医学成像技术，具有良好的超声分辨率和光学对比度，可以提供组织结构及功能信息。传统的PAT方法是定性成像，只能提供被吸收的激光能量密度的分布，而被吸收的激光能量密度是光学吸收系数与激光能量分布的乘积。光学吸收系数是组织的本质特征，它直接关系到组织的结构和功能信息。因此，从传统光声图像中消除激光能量分布的影响，进而恢复出组织光学吸收系数的分布具有重大意义。
基于动物不同组织对光的吸收系数与散射系数的不同，现有技术中从传统光声图像中消除激光能量分布的方式是对光声图像直接分割以获取不同区域，并假定各个分割区域内光学吸收系数、散射系数是一致的，再进行衰减校正计算，能够减少计算量。但这种方法需要人工进行手动图像分割，且光声图像对比度和结构信息较差，在现实中应用是非常困难的。此外，动物各器官内的光学参数是空间变化的，现有的对光声图像进行分割校正的方法不适应于这种变化。
因此针对现有技术不足，提供一种光声图像衰减校正方法以解决现有技术不足甚为必要。
发明内容
本发明的目的在于避免现有技术的不足之处而提供一种光声图像衰减校正方法，采用逐点迭代的方式进行，能够解决光声图像初始光能量分布不均匀的问题，且不需要人工介入，能够直接集成到各光声成像系统中。
本发明的上述目的通过以下技术措施实现：
提供一种光声图像衰减校正方法，包括步骤如下：
S1：将重建得到的光声图像P转为列向量形式p，将p代入线性系统方程中，所述线性系统方程为
p＝φμa 式(Ⅰ)；
其中，μa是组织吸收系数，p和μa的尺寸为mn×1，φ为mn×mn的稀疏光通量矩阵；
S2：采用最小二乘方式将线性系统方程构造为优化问题模型，所述优化问题模型为
Figure BDA0003226764440000021
S3：通过分步迭代梯度下降方法对优化问题模型进行迭代求解，获得组织吸收系数μa；
S4：将步骤S3获得的组织吸收系数μa变换为矩阵形式，并生成光通量分布矩阵Φ，将光声图像P和光通量分布矩阵Φ代入式(Ⅲ)得到校正后的光声图像P′；
Figure BDA0003226764440000022
其中，P的尺寸为m×n，Φ的尺寸为m×n。
优选的，所述步骤S3具体包括：
S3.1：定义总迭代次数为K，误差为ε，散射系数为μ′s，当前迭代次数为k，且0≤k≤K；
S3.2：初始化k＝1，组织吸收系数μa为0；
S3.3：使用组织吸收系数μa和散射系数μ′s更新稀疏光通量矩阵φ；
S3.4：固定φ不变，使用梯度下降求解优化问题模型，并更新组织吸收系数μa；
S3.5：判断第k次迭代时组织吸收系数的中间值μa的误差是否小于ε或k是否小于最大迭代次数K，若是，则输出k次迭代的组织吸收系数μa，否则，返回S3.3。
优选的，步骤S1中光声图像P使用基于模型的迭代方法或者反投影法重建获得。
优选的，步骤S1中光声图像P等于光通量分布矩阵Φ和组织吸收系数μa之积，具体是
P＝Φ×μa 式(Ⅳ)
其中，将P和μa转为列向量形式，将Φ转为稀疏矩阵形式，即可得到式(Ⅰ)所述的线性系统方程。
本发明的一种光声图像衰减校正方法，包括如下步骤：S1：将重建得到的的光声图像P转为列向量形式p，将p代入线性系统方程p＝φμa中；S2：采用最小二乘方式将线性系统方程构造为优化问题模型
Figure BDA0003226764440000031
S3：通过分步迭代梯度下降方法对优化问题模型进行迭代求解，获得组织吸收系数μa；S4：将步骤S3获得的组织吸收系数μa变换为矩阵形式，并生成光通量分布矩阵Φ，将光声图像P和光通量分布矩阵Φ代入
Figure BDA0003226764440000041
得到校正后的光声图像P′。该方法采用逐点迭代的方式进行，不需要人工介入就能够实现光通量分布的校正，进而实现光声图像衰减的校正，且校正速度快、能够直接集成到各光声成像系统中。
附图说明
利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。
图1是本发明一种光声图像衰减校正方法的流程图。
图2是实施例2中多光谱光声图像衰减校正的流程图。
图3是实施例3中活体实验不同波长下光通量校正结果示意图，其中，图3(a)是原始光声图像，图3(b)是吸收系数图，图3(c)是光通量分布图，图3(d)是校正后的光声图像。
图4是实施例3中活体实验波谱分离结果示意图，其中，图4(a)是原始光声图像波谱分离结果，图4(b)是校正后光声图像波谱分离结果，图4(c)是二者的差。
具体实施方式
结合以下实施例对本发明的技术方案作进一步说明。
实施例1
一种光声图像衰减校正方法，如图1所示，包括步骤如下：
S1：将重建得到的光声图像P转为列向量形式p，将p代入线性系统方程中，线性系统方程为
p＝φμa 式(Ⅰ)；
其中，μa是组织吸收系数，p和μa的尺寸为mn×1，φ为mn×mn的稀疏光通量矩阵。
光声图像是通过重建吸收激光脉冲产生的超声波而形成的，图像中
Figure BDA0003226764440000051
处的像素值表示为：
Figure BDA0003226764440000052
其中，α是系统响应系数，Г是热弹性膨胀系数，
Figure BDA0003226764440000053
是组织吸收系数，
Figure BDA0003226764440000054
是组织散射系数，Φ是光通量分布矩阵，
Figure BDA0003226764440000055
是空间位置向量。当消除掉常数项α和Г之后，光声图像等于光通量分布和组织吸收系数之积。
S2：采用最小二乘方式将线性系统方程构造为优化问题模型，优化问题模型为
Figure BDA0003226764440000056
S3：通过分步迭代梯度下降方法对优化问题模型进行迭代求解，获得组织吸收系数μa；
S4：将步骤S3获得的组织吸收系数μa变换为矩阵形式，并生成光通量分布矩阵Φ，将光声图像P和光通量分布矩阵Φ代入式(Ⅲ)得到校正后的光声图像P′；
Figure BDA0003226764440000057
其中，P的尺寸为m×n，Φ的尺寸为m×n。
本实施例中，步骤S3具体包括：
S3.1：定义总迭代次数为K，误差为ε，散射系数为μ′s，当前迭代次数为k，且0≤k≤K；
S3.2：初始化k＝1，组织吸收系数μa为0；
S3.3：使用组织吸收系数μa和散射系数μ′s更新稀疏光通量矩阵φ；
S3.4：固定φ不变，使用梯度下降求解优化问题模型，并更新组织吸收系数μa；
S3.5：判断第k次迭代时组织吸收系数的中间值μa的误差是否小于ε或k是否小于最大迭代次数K，若是，则输出k次迭代的组织吸收系数μa，否则，返回S3.3。采用逐点迭代的方式，能够适用于空间上任意点的光学参数都是变化的动物各器官的光声图像。
该光声图像衰减校正方法，采用逐点迭代的方式，能够校正光通量分布，且校正过程不需要人工介入，校正速度快，校正结果准确。
实施例2
一种光声图像衰减校正方法，其它特征与实施例1相同，不同之处在于：如图2所示，本实施例中进行衰减校正的光声图像为多光谱光声图像。需要说明的是，本实施例中多光谱光声图像具体由四种不同波长的光声图像构成，但实际应用中不限于四种，还可以是两种、三种或者五种等。
本实施例中，步骤S1中光声图像P使用基于模型的迭代方法或者反投影法重建获得。步骤S1中光声图像P等于光通量分布矩阵Φ和组织吸收系数μa之积，具体是
P＝Φ×μa 式(Ⅳ)
其中，将P和μa为列向量形式，将Φ转为稀疏矩阵形式，即可得到式(Ⅰ)所述的线性系统方程。
需要说明的是，光声图像P的重建方法不限于使用基于模型的迭代方法或者反投影法，还可以是其它方法，适应不同的光学成像系统。
该光声图像衰减校正方法，能够校正多光谱光声图像光通量分布，且能够集成到各光学成像系统中。
实施例3
一种光声图像衰减校正方法，其它特征与实施例1相同，不同之处在于：本实施例具体以活体小鼠的肾脏部位为实验对象进行光声图像衰减校正，使用商用小动物多光谱光声断层成像系统(MSOT inVision128,iThera Medical,Germany)进行实验数据采集。重建得到的原始光声图像如图3(a)所示，获得的吸收系数图如图3(b)所示，光通量分布图如图3(c)所示，校正后的光声图像如图3(d)所示。从图3中可以看出，光通量随着成像深度的增加而逐渐衰减，且校正后的光声图像与原始光声图像相比，在所有波长下都显示出增强的深层组织可见性，特别是对血管。原始光声图像的波谱分离结果如图4(a)所示，校正后光声图像的波谱分离结果如图4(b)所示，二者的差如图4(c)所示。从图4中可以看出，校正后的光声图像可以更准确的分离出深处含氧血和缺氧血的分布，增强了肾脏内部及中心动脉的血氧可见度。
该光声图像衰减校正方法，能够校正光通量分布，且能够提高深处组织可见度，提高光声成像的质量。另外，针对多波长光声图像的衰减校正该方法能够增加波谱分离结果的准确性，提高深处组织血氧分布的可见度。
最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (4)

1.一种光声图像衰减校正方法，其特征在于，包括如下步骤：
S1：将重建得到的光声图像P转为列向量形式p，将p代入线性系统方程中，所述线性系统方程为
p＝φμa式(I)；
其中，μa是组织吸收系数，p和μa的尺寸为mn×1，φ为mn×mn的稀疏光通量矩阵；
S2：采用最小二乘方式将线性系统方程构造为优化问题模型，所述优化问题模型为
Figure FDA0003226764430000011
S3：通过分步迭代梯度下降方法对优化问题模型进行迭代求解，获得组织吸收系数μa；
S4：将步骤S3获得的组织吸收系数μa变换为矩阵形式，并生成光通量分布矩阵Φ，将光声图像P和光通量分布矩阵Φ代入式(III)得到校正后的光声图像P′；
Figure FDA0003226764430000012
其中，P的尺寸为m×n，Φ的尺寸为m×n。
2.根据权利要求1所述的一种光声图像衰减校正方法，其特征在于，所述步骤S3具体包括：
S3.1：定义总迭代次数为K，误差为ε，散射系数为μ′s，当前迭代次数为k，且0≤k≤K；
S3.2：初始化k＝1，组织吸收系数μa为0；
S3.3：使用组织吸收系数μa和散射系数μ′s更新稀疏光通量矩阵φ；
S3.4：固定φ不变，使用梯度下降求解优化问题模型，并更新组织吸收系数μa；
S3.5：判断第k次迭代时组织吸收系数的中间值μa的误差是否小于ε或k是否小于最大迭代次数K，若是，则输出k次迭代的组织吸收系数μa，否则，返回S3.3。
3.根据权利要求1所述的一种光声图像衰减校正方法，其特征在于：步骤S1中光声图像P使用基于模型的迭代方法或者反投影法重建获得。
4.根据权利要求1所述的一种光声图像衰减校正方法，其特征在于：步骤S1中光声图像P等于光通量分布矩阵Φ和组织吸收系数μa之积，具体是
P＝Φ×μa式(IV)
其中，将P和μa转为列向量形式，将Φ转为稀疏矩阵形式，即可得到式(I)所述的线性系统方程。
CN202110974322.6A 2021-08-24 2021-08-24 一种光声图像衰减校正方法 Active CN113724351B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110974322.6A CN113724351B (zh) 2021-08-24 2021-08-24 一种光声图像衰减校正方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110974322.6A CN113724351B (zh) 2021-08-24 2021-08-24 一种光声图像衰减校正方法

Publications (2)

Publication Number Publication Date
CN113724351A true CN113724351A (zh) 2021-11-30
CN113724351B CN113724351B (zh) 2023-12-01

Family

ID=78677597

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110974322.6A Active CN113724351B (zh) 2021-08-24 2021-08-24 一种光声图像衰减校正方法

Country Status (1)

Country Link
CN (1) CN113724351B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115024739A (zh) * 2022-08-11 2022-09-09 之江实验室 生物体内格留乃森参数分布的测量方法、应用

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH09257694A (ja) * 1996-01-18 1997-10-03 Hamamatsu Photonics Kk 光ｃｔ装置及び光ｃｔによる画像再構成方法
CN102393969A (zh) * 2011-06-02 2012-03-28 西安电子科技大学 基于生物组织特异性的光学三维成像方法
WO2013159250A1 (zh) * 2012-04-28 2013-10-31 清华大学 一种动态荧光分子断层图像的重建方法
EP2660618A1 (en) * 2012-05-04 2013-11-06 Esaote S.p.A. Biomedical image reconstruction method
WO2018099321A1 (zh) * 2016-11-30 2018-06-07 华南理工大学 一种基于广义树稀疏的权重核范数磁共振成像重建方法
WO2018120329A1 (zh) * 2016-12-28 2018-07-05 深圳市华星光电技术有限公司 基于稀疏域重构的单帧图像超分辨重建方法及装置
CN110610528A (zh) * 2019-05-31 2019-12-24 南方医科大学 基于模型的双约束光声断层图像重建方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH09257694A (ja) * 1996-01-18 1997-10-03 Hamamatsu Photonics Kk 光ｃｔ装置及び光ｃｔによる画像再構成方法
CN102393969A (zh) * 2011-06-02 2012-03-28 西安电子科技大学 基于生物组织特异性的光学三维成像方法
WO2013159250A1 (zh) * 2012-04-28 2013-10-31 清华大学 一种动态荧光分子断层图像的重建方法
EP2660618A1 (en) * 2012-05-04 2013-11-06 Esaote S.p.A. Biomedical image reconstruction method
WO2018099321A1 (zh) * 2016-11-30 2018-06-07 华南理工大学 一种基于广义树稀疏的权重核范数磁共振成像重建方法
WO2018120329A1 (zh) * 2016-12-28 2018-07-05 深圳市华星光电技术有限公司 基于稀疏域重构的单帧图像超分辨重建方法及装置
CN110610528A (zh) * 2019-05-31 2019-12-24 南方医科大学 基于模型的双约束光声断层图像重建方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115024739A (zh) * 2022-08-11 2022-09-09 之江实验室 生物体内格留乃森参数分布的测量方法、应用

Also Published As

Publication number Publication date
CN113724351B (zh) 2023-12-01

Similar Documents

Publication Publication Date Title
Guan et al. Limited-view and sparse photoacoustic tomography for neuroimaging with deep learning
US11170258B2 (en) Reducing noise in an image
Rosenthal et al. Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography
Jose et al. Speed‐of‐sound compensated photoacoustic tomography for accurate imaging
Jeon et al. A deep learning-based model that reduces speed of sound aberrations for improved in vivo photoacoustic imaging
Deng et al. Deep learning in photoacoustic imaging: a review
CN109598722B (zh) 基于递归神经网络的图像分析方法
US20120302880A1 (en) System and method for specificity-based multimodality three- dimensional optical tomography imaging
Rajendran et al. Photoacoustic imaging aided with deep learning: a review
Gehrung et al. Co-registration of optoacoustic tomography and magnetic resonance imaging data from murine tumour models
CN102306385B (zh) 任意扫描方式下光声成像的图像重建方法
CN108451508B (zh) 基于多层感知机的生物自发荧光三维成像方法
Guo et al. AS-Net: fast photoacoustic reconstruction with multi-feature fusion from sparse data
CN110738701A (zh) 一种肿瘤三维定位系统
Steinberg et al. Superiorized photo-acoustic non-negative reconstruction (SPANNER) for clinical photoacoustic imaging
EP2780890B1 (en) System for creating a tomographic object image based on multiple imaging modalities
Lan et al. Y-Net: a hybrid deep learning reconstruction framework for photoacoustic imaging in vivo
Liu et al. Dictionary learning sparse-sampling reconstruction method for in-vivo 3D photoacoustic computed tomography
Zhang et al. MRI information-based correction and restoration of photoacoustic tomography
KR102037117B1 (ko) 구조유사도를 이용한 아티펙트 저감방법 및 프로그램 및 의료영상획득장치
CN113724351A (zh) 一种光声图像衰减校正方法
Yang et al. Recent advances in deep-learning-enhanced photoacoustic imaging
US11854124B2 (en) Real-time photoacoustic imaging using a precise forward model and fast iterative inverse
US7242004B2 (en) Image correction method, image correction apparatus, and image correction program
Lan et al. Cross-domain unsupervised reconstruction with equivariance for photoacoustic computed tomography

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