CN107180442A - 一种基于Renyi熵的光声图像重建前置滤波器 - Google Patents

一种基于Renyi熵的光声图像重建前置滤波器 Download PDF

Info

Publication number
CN107180442A
CN107180442A CN201710240539.8A CN201710240539A CN107180442A CN 107180442 A CN107180442 A CN 107180442A CN 201710240539 A CN201710240539 A CN 201710240539A CN 107180442 A CN107180442 A CN 107180442A
Authority
CN
China
Prior art keywords
photoacoustic signal
photoacoustic
signal
image
mrow
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.)
Pending
Application number
CN201710240539.8A
Other languages
English (en)
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.)
Taiyuan University of Technology
Original Assignee
Taiyuan University of Technology
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 Taiyuan University of Technology filed Critical Taiyuan University of Technology
Priority to CN201710240539.8A priority Critical patent/CN107180442A/zh
Publication of CN107180442A publication Critical patent/CN107180442A/zh
Pending legal-status Critical Current

Links

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/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/30—Noise filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
    • G06V10/443—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components by matching or filtering
    • 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/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/416—Exact reconstruction

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Multimedia (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明属于信号处理领域，具体涉及利用Renyi熵对光声信号进行杂波的过滤，从而提高光声信号的信噪比和光声成像质量。一种基于Renyi熵的光声图像重建前置滤波器，包括如下步骤：第一步，光声信号获取；第二步，光声信号的时频分布求解；第三步，各光声信号点的Renyi熵求解；第四步，阈值确定；第五步，滤波处理；第六步，用经过滤波处理的光声信号进行图像重建。与现有技术相比，本发明具有以下优点：（1）利用本发明前置滤波器处理后的光声信号进行光声图像重建时，重建图像的信噪比有了明显提高。（2）利用本发明前置滤波器处理后的光声信号进行光声图像重建时，重建图像的均方误差明显减小。

Description

一种基于Renyi熵的光声图像重建前置滤波器
技术领域
本发明属于信号处理领域，具体涉及利用Renyi熵对光声信号进行杂波的过滤，从而提高光声信号的信噪比和光声成像质量。
背景技术
在光声图像重建的过程中，影响成像质量的因素有很多，其中最重要的两个因素为原始光声信号的信噪比和图像重建算法。然而，目前有关光声图像重建的研究更多的偏向对算法的不断优化与改进，却忽略了光声信号本身质量的影响。光声信号采自于超声换能器，它是在激光脉冲产生的同时开始工作的。理论上讲，超声换能器接收到的光声信号只包含成像组织的信息，组织以外的区域不应该有光声信号被接收到。但是在实际实验中，由于测试环境、传输媒质以及散射等影响，传感器接收到的光声信号往往会携带大量的噪音信号，造成后续重构图像的分辨率低下，对比度不足。因此，在图像重建之前对原始光声信号进行滤波处理是非常有必要的。
目前，光声信号滤波处理的基本方法包括：滤波反投影算法，正则化方法以及模极大值法等。如，华南师范大学的曾亚光等在2005 年提出的滤波反投影重建算法，该算法可以适当提高重建光声信号的强度和信噪比，但同时也会带来图像伪迹问题；天津大学的宋智源等在2008年提出利用正则化方法对光声信号进行还原，该算法为更好更准确地重现组织内部的光吸收分布结构提供了一条行之有效的途径，但正则化参数一般为噪声的函数，在噪声未知的情况下，光声信号还原过程存在较大的主观性；东北大学肖倩等在2011年提出利用模极大值法对光声信号中的杂波进行滤除，该算法可以很好地保留间断点的信息，但其输出信噪比受算法分解尺度的影响较大。
发明内容
本发明的目的在于利用Renyi熵对原始光声信号进行滤波处理来提高光声信号质量，从而能够更好地进行图像重建。
本发明是采用以下技术方案实现的：一种基于Renyi熵的光声图像重建前置滤波器，包括如下步骤：
第一步光声信号获取
在二维光声成像中，用短脉冲激光垂直于成像平面照射生物组织，超声换能器在成像平面上采集光声信号并且沿着圆周平面以步进为n进行扫描，共采样360/n个点，在每一个采样点也即每个光声信号点上采集m个光声信号；
第二步光声信号的时频分布求解
采用短时傅里叶变换即加窗的傅里叶变换来求解每一个光声信号点的时频分布：
s(t,f)＝∫x(θ)h(t-θ)e-j2πfθdθ (1)
式中s(t，f)代表某个光声信号点的时频分布，x(θ)为待处理的该光声信号点采集的m个光声信号，h(t-θ)为选取的窗函数；通过式(1) 可以得到所有光声信号点的时频分布；这一步利用该窗函数和光声信号的乘积来定义局部相关函数，然后对这个局部相关函数进行傅里叶变换；
第三步各光声信号点的Renyi熵求解
Renyi熵是信息熵的一种形式，它所表征的是信号的复杂度，各光声信号点α阶的Renyi熵可以表示为：
其中s(t，f)表示某个光声信号点的时频分布，α为整数且≥2；根据(2)式就可以得到各个光声信号点的Renyi熵；
第四步阈值确定
根据第三步所得各光声信号点的Renyi熵，勾画熵值直方图，并选取直方图中出现频次最小的熵值作为分割阈值；
第五步滤波处理
利用前一步所得分割阈值对所采集的所有光声信号进行分割，将 Renyi熵高于分割阈值的光声信号滤除；
第六步用经过滤波处理的光声信号进行图像重建。
Renyi熵，作为信息学理论中信号复杂性的一个测度，已被用于多个领域信号杂波的去除。Renyi熵用于表征信号复杂度，在信号的随机处(杂波)取值较大，在固定频率统计性较好的地方(即，有用信号)取值较小。利用上述特性，不需要满足其他条件，就可以确定一个用于分离光声信号杂波的阈值，进而对信号进行阈值分割，提高信噪比。所述窗函数可采用汉宁窗函数或海明窗函数。
进一步的，还包括第七步对重建后的图像计算信噪比，将重建后图像的信噪比与设定的信噪比X进行比较，如果小于X，则返回到第二步对光声信号重新进行处理，直到信噪比符合设定要求。
重新处理时，各步骤所有的参数都是不变的，对经过阈值分割后的光声信号再进行Renyi熵计算，再在这些熵值所对应的直方图选取出现频次最小的熵值作为新的分割阈值，对反馈回的光声信号利用阈值进行进一步的滤波。
与现有技术相比，本发明具有以下优点：
(1)利用本发明前置滤波器处理后的光声信号进行光声图像重建时，重建图像的信噪比有了明显提高。
(2)利用本发明前置滤波器处理后的光声信号进行光声图像重建时，重建图像的均方误差明显减小。
总体而言，基于Renyi熵的光声图像重建前置滤波器设计能有效改善光声图像重建的质量。
附图说明
图1为本发明基于光声图像重建前置滤波器的设计流程图。
图2为本发明基于光声图像重建前置滤波器对某一采样位置头发丝光声信号的滤波处理效果。
图3为某一采样位置头发丝光声信号Renyi熵对应的直方图。
图4为本发明基于光声图像重建前置滤波器对多个采样位置小鼠大脑皮层血管光声信号滤波处理前后的对比图。
图5为本发明基于光声图像重建前置滤波器对小鼠大脑皮层血管光声信号进行滤波处理前后图像重建的对比结果。其中，(a)为滤波器处理前的小鼠大脑皮层血管光声图像，(b)为滤波器处理后的小鼠大脑皮层血管光声图像。下插图为滤波器处理前后小鼠大脑皮层血管光声图像的信噪比和均方误差附表。
具体实施方式
下面结合附图和具体实施方式对本发明的应用原理做进一步的阐述。
实验平台采用的激光源为Surelite I-20 YAG激光器，波长为 532nm，脉冲重复频率为20Hz，脉冲宽度为5-7ns，单脉冲能量为 420mJ。采用浸入式非聚焦超声换能器接收光声信号，中心频率为 3.5MHz，有效带宽为1.12MHz，有效直径为9.525nm。采用前置放大器来对接收到的光声信号进行放大，放大器采用的是PREAMP2-D 前置放大器。光声信号经过前置放大器放大后进入示波器进行采样，示波器采用的是DPO 3012型数字示波器，最高采样率为2.5GS/s。步进电机的型号为RSA100型数控转台，精确控制换能器采样的角度。激光器、步进电机和示波器分别通过RS232，USB和PCI-GPIB 接口卡和计算机相连。
实施例1
利用本发明的前置滤波器对头发丝的光声信号进行处理。超声换能器成像平面上采集光声信号时选择以步进为1°进行扫描，探测器在每一个位置采集10000个点后旋转至下一个位置，进行全角度的圆周扫描。扫描结果进入示波器之后，用本发明的前置滤波器对采集到的光声信号进行滤波。首先，在公式(1)中利用汉宁窗函数来计算头发丝光声信号的时频分布。根据公式(2)并选取α＝2来计算头发丝各光声信号点对应的2阶Renyi熵值，以每个熵值出现的频数为依据，得到头发丝光声信号Renyi熵对应的直方图，如图2插图所示。选取直方图中低熵范围内出现频数最小的点对应的Renyi熵(0.6J)作为头发丝光声信号的分割阈值。利用该分割阈值对头发丝的光声信号进行分割，得到滤波后的头发丝光声信号结果(图2)。经对比计算，利用本发明的前置滤波器滤波处理之后，头发丝光声信号的信噪比较未处理之前提高了48.26％，均方误差减少了31％。
实施例2
利用本发明的前置滤波器对小鼠大脑皮层血管的光声信号进行处理。超声换能器成像平面上采集超声信号时选择以步进为2°进行扫描，探测器在每一个位置采集5000个点后旋转至下一个位置，进行全角度的圆周扫描。扫描结果进入示波器之后，用本发明的前置滤波器对采集到的小鼠大脑皮层血管光声信号进行滤波。首先，在公式 (1)中利用海明窗函数来计算小鼠大脑皮层血管光声信号的时频分布。根据公式(2)选取α＝3来计算小鼠大脑皮层血管各光声信号点对应的3 阶Renyi熵值，以每个熵值出现的频数为依据，得到小鼠大脑皮层血管光声信号Renyi熵对应的直方图。选取直方图中低熵范围内出现频数最小的点对应的Renyi熵(0.8J)作为小鼠大脑皮层血管光声信号的分割阈值。利用该分割阈值对小鼠大脑皮层血管的光声信号进行分割，得到滤波前后小鼠大脑皮层血管光声信号的对比结果，如图3所示。最后利用延迟叠加算法来对处理后的小鼠大脑皮层血管的光声信号进行图像重建(图4)。结果显示，利用本发明的前置滤波器对小鼠大脑皮层血管的光声信号进行滤波处理之后，重建图像的信噪比较滤波之前提高了51.63％，均方误差减少了69％。
以上所述仅为本发明的实施例而已，并不用于限制本发明。本发明可以有各种合适的更改和变化。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之。

Claims (2)

1.一种基于Renyi熵的光声图像重建前置滤波器，其特征在于，包括如下步骤：
第一步 光声信号获取
在二维光声成像中，用短脉冲激光垂直于成像平面照射生物组织，超声换能器在成像平面上采集光声信号并且沿着圆周平面以步进为n进行扫描，共采样360/n个点，在每一个采样点也即每个光声信号点上采集m个光声信号；
第二步 光声信号的时频分布求解
采用短时傅里叶变换即加窗的傅里叶变换来求解每一个光声信号点的时频分布：
s(t,f)＝∫x(θ)h(t-θ)e-j2πfθdθ (1)
式中s(t，f)代表某个光声信号点的时频分布，x(θ)为待处理的该光声信号点采集的m个光声信号，h(t-θ)为选取的窗函数；通过式(1)可以得到所有光声信号点的时频分布；
第三步 各光声信号点的Renyi熵求解
各光声信号点α阶的Renyi熵可以表示为：
<mrow> <msub> <mi>R</mi> <mi>&amp;alpha;</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>1</mn> <mo>-</mo> <mi>&amp;alpha;</mi> </mrow> </mfrac> <msub> <mi>log</mi> <mn>2</mn> </msub> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mi>s</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>t</mi> <mi>d</mi> <mi>f</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
其中s(t，f)表示某个光声信号点的时频分布，α为整数且≥2；根据(2)式就可以得到各个光声信号点的Renyi熵；
第四步 阈值确定
根据第三步所得各光声信号点的Renyi熵，勾画熵值直方图，并选取直方图中出现频次最小的熵值作为分割阈值；
第五步 滤波处理
利用前一步所得分割阈值对所采集的所有光声信号进行分割，将Renyi熵高于分割阈值的光声信号滤除；
第六步 用经过滤波处理的光声信号进行图像重建。
2.如权利要求1所述的一种基于Renyi熵的光声图像重建前置滤波器，其特征在于，还包括第七步对重建后的图像计算信噪比，将重建后图像的信噪比与设定的信噪比X进行比较，如果小于X，则返回到第二步对光声信号重新进行处理，直到信噪比符合设定要求。
CN201710240539.8A 2017-04-13 2017-04-13 一种基于Renyi熵的光声图像重建前置滤波器 Pending CN107180442A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710240539.8A CN107180442A (zh) 2017-04-13 2017-04-13 一种基于Renyi熵的光声图像重建前置滤波器

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710240539.8A CN107180442A (zh) 2017-04-13 2017-04-13 一种基于Renyi熵的光声图像重建前置滤波器

Publications (1)

Publication Number Publication Date
CN107180442A true CN107180442A (zh) 2017-09-19

Family

ID=59830949

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710240539.8A Pending CN107180442A (zh) 2017-04-13 2017-04-13 一种基于Renyi熵的光声图像重建前置滤波器

Country Status (1)

Country Link
CN (1) CN107180442A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110503044A (zh) * 2018-12-03 2019-11-26 神盾股份有限公司 指纹传感器及其指纹感测方法
WO2020082934A1 (zh) * 2018-10-25 2020-04-30 南京大学 一种抑制伪影的光声图像重建方法
CN112842264A (zh) * 2020-12-31 2021-05-28 哈尔滨工业大学（威海） 多模态成像中数字滤波方法、装置和多模态成像技术***

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040022438A1 (en) * 2002-08-02 2004-02-05 Hibbard Lyndon S. Method and apparatus for image segmentation using Jensen-Shannon divergence and Jensen-Renyi divergence
CN103584835A (zh) * 2013-09-24 2014-02-19 南京大学 一种基于压缩感知的光声图像重建方法
CN104899902A (zh) * 2015-04-14 2015-09-09 华北电力大学（保定） 一种血管内光声二维图像的重建方法
CN106056618A (zh) * 2016-06-08 2016-10-26 湖南文理学院 基于Renyi交叉熵与高斯分布的图像阈值分割方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040022438A1 (en) * 2002-08-02 2004-02-05 Hibbard Lyndon S. Method and apparatus for image segmentation using Jensen-Shannon divergence and Jensen-Renyi divergence
CN103584835A (zh) * 2013-09-24 2014-02-19 南京大学 一种基于压缩感知的光声图像重建方法
CN104899902A (zh) * 2015-04-14 2015-09-09 华北电力大学（保定） 一种血管内光声二维图像的重建方法
CN106056618A (zh) * 2016-06-08 2016-10-26 湖南文理学院 基于Renyi交叉熵与高斯分布的图像阈值分割方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张莉等: "基于S变换和Renyi熵的舒张期心音信号分析", 《杭州电子科技大学学报》 *
李伟等: "提高光声成像质量的Renyi熵方法", 《生物医学工程学进展》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020082934A1 (zh) * 2018-10-25 2020-04-30 南京大学 一种抑制伪影的光声图像重建方法
US11763500B2 (en) 2018-10-25 2023-09-19 Nanjing University Photoacoustic image reconstruction method for suppressing artifacts
CN110503044A (zh) * 2018-12-03 2019-11-26 神盾股份有限公司 指纹传感器及其指纹感测方法
CN112842264A (zh) * 2020-12-31 2021-05-28 哈尔滨工业大学（威海） 多模态成像中数字滤波方法、装置和多模态成像技术***
CN112842264B (zh) * 2020-12-31 2023-04-25 哈尔滨工业大学（威海） 多模态成像中数字滤波方法、装置和多模态成像技术***

Similar Documents

Publication Publication Date Title
CN102551810B (zh) 一种多通道同步实时数字化光声成像装置及方法
CN1160022C (zh) 计算机x射线断层造影机
CN107180442A (zh) 一种基于Renyi熵的光声图像重建前置滤波器
WO2021232653A1 (zh) 一种结合滤波反投影算法和神经网络的pet图像重建算法
Zeng et al. High antinoise photoacoustic tomography based on a modified filtered backprojection algorithm with combination wavelet
CN101685158B (zh) 基于隐马尔科夫树模型的sar图像去噪方法
Zhou et al. A noise reduction method for photoacoustic imaging in vivo based on EMD and conditional mutual information
CN109803588A (zh) 流体路径的增强分辨率超声成像
CN102073994B (zh) 基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法
CN107346541A (zh) 一种基于超声射频时间序列小波分析的组织定征方法
CN111415408A (zh) 一种超声空化的微秒级多尺度时空成像及特征图谱计算方法与***
Zhou et al. Wavelet de-noising method with adaptive threshold selection for photoacoustic tomography
CN111248858B (zh) 一种基于频域波数域的光声断层成像重建方法
TWI652048B (zh) 超音波影像均勻度檢測方法及其系統
Kong et al. Investigation on reconstruction for frequency domain photoacoustic imaging via TVAL3 regularization algorithm
CN117237473A (zh) 一种基于扩散模型的光声断层成像稀疏重建方法
US20220028128A1 (en) Photoacoustic image reconstruction method for suppressing artifacts
CN112861731B (zh) 一种基于参数寻优的超声信号去噪方法
US20110044516A1 (en) Contrast improvement method and system for photoacoustic imaging
CN111462077B (zh) 一种利用非线性信息熵定征生物组织的方法
CN115049753B (zh) 一种基于无监督深度学习的锥束ct伪影校正方法
CN113974560A (zh) 环形光声层析***稀疏阵元优化选择及压缩感知成像方法
CN108573474B (zh) 一种采用逆卷积运算的光声图像优化方法
CN110084772B (zh) 基于弯曲波的mri/ct融合方法
Merino et al. Shear Wave Speed Estimation for Crawling Wave Sonoelastography Using the Short-Time Fourier Transform

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
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20170919