CN116167948B - 一种基于空变点扩散函数的光声图像复原方法及系统 - Google Patents

一种基于空变点扩散函数的光声图像复原方法及系统 Download PDF

Info

Publication number
CN116167948B
CN116167948B CN202310432349.1A CN202310432349A CN116167948B CN 116167948 B CN116167948 B CN 116167948B CN 202310432349 A CN202310432349 A CN 202310432349A CN 116167948 B CN116167948 B CN 116167948B
Authority
CN
China
Prior art keywords
image
psf
deconvolution
space
representing
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
CN202310432349.1A
Other languages
English (en)
Other versions
CN116167948A (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.)
Institute of Artificial Intelligence of Hefei Comprehensive National Science Center
Original Assignee
Institute of Artificial Intelligence of Hefei Comprehensive National Science Center
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 Institute of Artificial Intelligence of Hefei Comprehensive National Science Center filed Critical Institute of Artificial Intelligence of Hefei Comprehensive National Science Center
Priority to CN202310432349.1A priority Critical patent/CN116167948B/zh
Publication of CN116167948A publication Critical patent/CN116167948A/zh
Application granted granted Critical
Publication of CN116167948B publication Critical patent/CN116167948B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • 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
    • G06T5/00—Image enhancement or restoration
    • G06T5/80—Geometric correction
    • 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/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
    • G06V10/774—Generating sets of training patterns; Bootstrap methods, e.g. bagging or boosting
    • 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/10004—Still image; Photographic image
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Multimedia (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于空变点扩散函数的光声图像复原方法及系统，包括S1：获取部分不同空间位置处的PSF图像，组成PSF图像样本集，基于PSF图像样本集构建空变系统PSF的表征模型，获得空变系统所有空间位置处的PSF图像；S2：构建RL解卷积算法，利用空变系统PSF的表征模型和RL解卷积算法对所采集到的光声图像进行解卷积，得到解卷积图像；S3：将解卷积图像绕图像中心旋转设定角度，得到旋转后图像，将旋转后图像与解卷积图像相加取平均，得到去除了条纹伪影的平均图像，利用RL解卷积算法对所述平均图像进行复原，得到最终的复原图像；该光声图像复原方法及系统提升了复原图像的质量。

Description

一种基于空变点扩散函数的光声图像复原方法及系统
技术领域
本发明涉及图像复原技术领域，尤其涉及一种基于空变点扩散函数的光声图像复原方法及系统。
背景技术
光声断层成像是一种新兴的生物医学成像技术，它结合了光学成像高对比度和声学成像在深层组织中高分辨率的优点。传统的光声断层图像重建算法主要是基于理想的探测器（全视角、无限带宽的点状探测器），然而，由于制造工艺、经济成本和成像空间的限制，实际的探测器很难满足这一要求，从而导致原本清晰的重建图像被模糊；在图像处理领域，图像被模糊的过程可以表述为清晰图像先与点扩散函数（Point Spread Function，PSF）卷积再叠加噪声，因此其逆过程被称为图像复原或图像解卷积，在光声断层成像中，图像复原的方法分为以下两类：信号域的复原方法和图像域的复原方法。
信号域的复原方法主要是对探测器采集到的光声信号进行解卷积，该方法将探测器采集到的光声信号作为原始光声信号与探测器脉冲响应的卷积，于是当探测器的脉冲响应已知时，就可以先利用解卷积算法得到原始光声信号，再通过滤波反投影、迭代、时间反演等传统图像重建算法将原始光声信号重建为图像。
部分文献【C. Zhang, and Y. Wang, "Deconvolution reconstruction offull-view and limited-view photoacoustic tomography: a simulation study," J.Opt. Soc. Am. A 25, 2436-2443 (2008)或者C. Zhang, C. Li, and L. V. Wang, "Fast and robust deconvolution-based image reconstruction for photoacoustictomography in circular geometry: experimental validation," IEEE Photonics J.2, 57-66 (2010)】还通过利用更准确的脉冲响应从而提升图像质量，如所述文献提出的解卷积图像重建算法（Deconvolution Reconstruction, DR），该方法在有限视角的图像重建中，有更强的抵抗噪声的能力，但重建误差受被检测物体的尺寸所影响。
图像域的复原方法是对重建得到的光声图像解卷积，该方法的关键集中于获取准确的PSF和设计解卷积算法。如文献【T. Jetzfellner, and V. Ntziachristos, "Performance of Blind Deconvolution in Optoacoustic Tomography," J. Innov.Opt. Health Sci. 04, 385-393 (2011)】选用了盲解卷积的方法来降低图像复原的难度，即从重建图像中估计出PSF，并在复原过程不断地更新估计出的PSF，盲解卷积算法可以提高图像的分辨率，但无法保证估计出的PSF的准确性，其复原效果具有局限性。
尽管上述方法提升了重建图像的质量，但在实际应用过程中，探测器的脉冲响应通常并不知道，而准确地测量出脉冲响应极具挑战。另外，受到探测器口径、几何形状等因素的影响，系统的PSF通常是空间变化的，而以上解卷积方法认为PSF是空间不变的，即各个位置处的PSF是相同的。正因为忽略了探测器的PSF在不同空间位置也不同的特性，即空变性，上述方法提升图像质量的能力十分有限。
发明内容
基于背景技术存在的技术问题，本发明提出了一种基于空变点扩散函数的光声图像复原方法及系统，利用空间变化的点扩散函数对重建图像进行解卷积，提升了图像的质量，能更好地恢复图像细节。
本发明提出的一种基于空变点扩散函数的光声图像复原方法，包括如下步骤：
S1：获取部分不同空间位置的PSF图像，组成PSF图像样本集，基于PSF图像样本集构建空变系统PSF的表征模型，获得空变系统所有位置处的PSF图像；
构建空变系统PSF的表征模型的目的是为了利用部分PSF图像得到所有不同空间位置处的PSF图像；
S2：构建基于稀疏对数梯度正则约束的RL解卷积算法，利用所述空变系统PSF的表征模型和所述稀疏对数梯度正则约束的RL解卷积算法对所采集到的光声图像进行解卷积，得到解卷积图像；
S3：将所述解卷积图像绕图像中心旋转设定角度后，得到旋转后图像，将所述旋转后图像与所述解卷积图像相加取平均，得到去除了条纹伪影的平均图像，利用所述稀疏对数梯度正则约束的RL解卷积算法对所述平均图像进行复原，得到最终的复原图像；
具体为，旋转平均操作可以去除条纹伪影，因此得到的平均图像中条纹伪影大部分已经被去除；相比于解卷积图像，平均图像中虽然消除了伪影，但旋转平均操作会导致平均图像出现轻微模糊；因此再次利用所述稀疏对数梯度正则约束的RL解卷积算法对所述平均图像进行复原以消除旋转平均操作带来的模糊；
基于稀疏对数梯度正则约束的RL解卷积算法的构建过程如下：
其中，和/>分别表示傅里叶变换和傅里叶逆变换，k表示迭代次数，符号“^”表示频域形式的结果，/>表示第i个特征-点扩散函数，/>表示将/>逆序排列的结果，/>表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，I表示退化图像，O表示清晰图像，/>和/>为常参数；/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代时的正则约束，/>表示利用第/>次迭代估计出的清晰图像和空变系统PSF的表征模型预测的退化图像，并在频域表示出该结果，/>表示第/>次迭代预测出的退化图像与真实退化图像的残差，/>表示频域形式的修正项，用来优化/>以得到/>，·表示乘积。
进一步地，在步骤S1：获取部分不同空间位置的PSF图像，组成PSF图像样本集，基于PSF图像样本集构建空变系统PSF的表征模型，获取空变系统所有位置处的PSF图像中中，具体包括：
利用主成分分析拟合已知位置的PSF图像，通过一系列特征-点扩散函数的线性组合表征所述已知位置的PSF图像；
利用径向基函数插值法按列对所述特征-点扩散函数对应的系数矩阵进行插值，得到所有位置的PSF图像，基于所述空变系统PSF的表征模型简化空变系统的成像模型。
进一步地，在所述利用主成分分析拟合已知位置的PSF图像，通过一系列特征-点扩散函数的线性组合表征所述已知位置的PSF图像中，具体包括：
记录已知的N个位置处的PSF图像及其对应的位置坐标，以坐标为中心对N个位置处的PSF图像进行裁剪，得到裁剪后的PSF图像；
对所述裁剪后的PSF图像进行归一化处理，得到归一化后的PSF图像；
计算所述归一化后的PSF图像之间的协方差，得到协方差矩阵；
对所述协方差矩阵进行特征值分解，得到所述协方差矩阵的N个特征值和N个特征向量/>；
利用所述特征向量构造特征-点扩散函数，通过所述特征-点扩散函数的线性组合表征所述已知位置的PSF图像。
进一步地，基于所述空变系统PSF的表征模型构建空变系统的成像模型，所述空变系统PSF的表征模型P简化如下：
所述空变系统的成像模型简化为下式：
其中，I表示退化图像，O表示清晰图像，表示噪声项，*表示卷积，/>表示第i个特征-点扩散函数，/>表示表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，/>表示物面的坐标，/>表示像面的坐标，/>表示物面和像面坐标。
进一步地，在步骤S2：构建基于稀疏对数梯度正则约束的RL解卷积算法中，具体公式还包括如下：
采用贝叶斯最大后验估计表述图像复原问题，选取O函数使得后验密度函数最大，等价于最小化后验密度函数的负对数；
假设似然密度函数服从泊松分布，得到经典RL算法；
根据Markov随机场理论和Hammersley-Clifford定理，采用重尾分布函数对清晰图像的先验密度函数建模得到的正则项称为稀疏对数梯度正则项；
基于所述稀疏对数梯度正则项构建稀疏对数梯度正则约束的RL解卷积算法。
进一步地，所述O函数公式如下：
所述等价于最小化后验密度函数的负对数公式如下：
所述似然密度函数公式如下：
所述清晰图像的先验密度函数公式如下：
所述稀疏对数梯度正则项具体如下：
其中，表示清晰图像的先验密度函数，·表示乘积，/>表示退化图像中坐标所在位置对应的像素值求阶乘。
进一步地，在步骤S3：将所述解卷积图像绕图像中心旋转设定角度后，得到旋转后图像，将所述旋转后图像与所述解卷积图像相加取平均，得到去除了条纹伪影的平均图像，利用所述稀疏对数梯度正则约束的RL解卷积算法对所述平均图像进行复原以消除旋转平均操作带来的模糊，得到最终的复原图像中，具体包括：
将所述解卷积图像绕图像中心按顺时针旋转设定角度后，得到顺时针旋转图像，将所述顺时针旋转图像与所述解卷积图像相加取平均，得到第一旋转图像；
将所述解卷积图像绕图像中心按逆时针旋转设定角度后，得到逆时针旋转图像，将所述逆时针旋转图像与所述解卷积图像相加取平均，得到第二旋转图像；
将第一旋转图像和第二旋转图像相加取平均，得到去除了条纹伪影的平均图像；
利用所述稀疏对数梯度正则约束RL解卷积算法对所述平均图像进行复原以消除旋转平均操作带来的模糊，得到最终的复原图像。
一种基于空变点扩散函数的光声图像复原系统，其特征在于，包括构建模块、解卷积模块和去伪影模块；
构建模块用于获取部分不同空间位置的PSF图像，组成PSF图像样本集，基于PSF图像样本集构建空变系统PSF的表征模型，获得空变系统所有位置处的PSF图像；
解卷积模块用于构建基于稀疏对数梯度正则约束的RL解卷积算法，利用所述空变系统PSF的表征模型和所述稀疏对数梯度正则约束的RL解卷积算法对所采集到的光声图像进行解卷积，得到解卷积图像；
去伪影模块用于将所述解卷积图像绕图像中心旋转设定角度后，得到旋转后图像，将所述旋转后图像与所述解卷积图像相加取平均，得到去除了条纹伪影的平均图像，利用所述稀疏对数梯度正则约束的RL解卷积算法对所述平均图像进行复原，得到最终的复原图像；
基于稀疏对数梯度正则约束的RL解卷积算法的构建过程如下：
其中，和/>分别表示傅里叶变换和傅里叶逆变换，k表示迭代次数，符号“^”表示频域形式的结果，/>表示第i个特征-点扩散函数，/>表示将/>逆序排列的结果，/>表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，I表示退化图像，O表示清晰图像，/>和/>为常参数，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代时的正则约束，/>表示利用第/>次迭代估计出的清晰图像和空变系统PSF的表征模型预测的退化图像，并在频域表示出该结果，/>表示第/>次迭代预测出的退化图像与真实退化图像的残差，/>表示频域形式的修正项，用来优化/>以得到/>，·表示乘积。
本发明提供的一种基于空变点扩散函数的光声图像复原方法及系统的优点在于：本发明结构中提供的一种基于空变点扩散函数的光声图像复原方法及系统，使用SVPSF对图像进行复原更符合实际情况，相比使用空不变的PSF的复原算法，该方法进一步提升了图像的质量，能更好地恢复图像细节；提出的SLG正则项，其抑制噪声的能力和复原效果优于传统的Tikhonov正则项；所提出的去伪影方法，适用于使用环形阵列换能器采集得到的光声图像，且该方法简单易操作，去除条纹伪影的同时能较好地保留图像细节。
附图说明
图1为本发明的结构示意图；
图2为图像复原的流程示意图；
图3为SVPSF建模的流程图；
图4为去除光声图像中的条纹伪影的流程图；
图5为光声图像的梯度分布图，其中（a）表示实际采集的光声图像，（b）表示光声图像真实的水平梯度分布曲线和不同正则项对应的位势函数曲线；
图6为小鼠腹部血管图像的复原示意图，其中，（a）表示光声图像；（b）表示解卷积图像；（c）去伪影图像；（d）表示（a）中实线框的放大图，（e）表示（b）中实线框的放大图，（f）表示（c）中实线框的放大图，（a）、（b）、（c）中的实线框对应于光声图像的同一位置；
图7为成年男性手指横截面图像的复原示意图，其中，（a）表示光声图像；（b）表示解卷积图像；（c）去伪影图像；（d）表示（a）中实线框的放大图，（e）表示（b）中实线框的放大图，（f）表示（c）中实线框的放大图，（a）、（b）、（c）中的实线框对应于光声图像的同一位置。
具体实施方式
下面，通过具体实施例对本发明的技术方案进行详细说明，在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其他方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似改进，因此本发明不受下面公开的具体实施的限制。
目前为止，主流的图像复原方法都基于PSF空间不变的假设进行解卷积。为了得到更高质量的复原图像，本实施例提出了一种利用空间变化的点扩散函数（SpatiallyVariant Point Spread Function，SVPSF）对重建图像进行解卷积的复原方法。本实施例解决的主要技术问题包括：
（1）准确获取不同空间位置处的PSF图像；
光声图像复原需要知道每个位置的PSF图像，当成像面尺寸过大时，测得每个位置对应的PSF图像非常繁琐且耗时。本实施例利用主成分分析的方法将部分已知位置处的PSF图像表征为一系列基函数的线性组合，根据提出的SVPSF表征模型和插值方法准确计算出其他位置处的PSF图像。
（2）提出新的正则项以抑制噪声，提升复原图像质量；
利用已得到的以及估计出的SVPSF和经典的Richardson-Lucy（RL）算法，在频域进行图像复原，由于解卷积本身是一个逆问题，具有病态性，在复原过程中，噪声会被放大，因此必须选用合适的正则项以抑制噪声。鉴于自然图像的梯度的稀疏性，本实施例利用重尾分布函数对图像的先验分布建模，提出一种全新的正则项，稀疏对数梯度（SparseLogarithmic Gradient，SLG）正则项，本实施例成功利用该正则项抑制了噪声，得到了更优的复原图像。
（3）去除光声图像中的条纹伪影；
由于非理想的成像条件，使用环形阵列换能器重建出的光声图像本身存在比较严重的条纹伪影，本实施例中的解卷积方法虽然可以提高图像分辨率和抑制噪声，但无法有效去除条纹伪影，因此，本实施例进一步根据条纹伪影的分布特性，通过旋转平均的方式去除了条纹伪影；同时再次利用本实施例中的解卷积算法消除旋转平均带来的图像模糊，实现了在消除伪影的同时也保留了图像细节。
需要说明的是，本实施例提出的解卷积方法是对重建出的光声图像进行后处理，因此本申请中提及的PSF图像皆表示用重建算法获得的PSF图像。另外，去伪影方法仅适用于处理稀疏阵元的环阵换能器采得的光声图像。可以预见的是，当使用阵元密集的换能器时，所获得的光声信号数据越丰富，重建得到的光声图像的质量也会越高，条纹伪影几乎不可见。
如图1至7所示，为了具体解决以上主要技术问题，本发明提出的一种基于空变点扩散函数的光声图像复原方法，包括如下步骤S1至S3：
S1：获取部分不同空间位置的PSF图像，组成PSF图像样本集，基于PSF图像样本集构建空变系统PSF的表征模型，获取空变系统所有空间位置处的PSF；
SVPSF建模模块利用部分已知位置处的PSF构造出了其他所有位置对应的PSF；将SVPSF建模得到的结果与采集到的光声图像输入解卷积模块获得解卷积图像。
其中SVPSF建模模块的主要目的是获取整个重建位置处的PSF图像，如图3所示，分为以下两个步骤：一是利用主成分分析拟合已知位置的PSF图像，将其表示为一系列特征-点扩散函数的线性组合，得到PSF的表征模型；二是对模型参数进行插值，从而估计出未知位置处的PSF图像。
其中，利用主成分分析拟合已知位置的PSF图像的具体过程如下
①记录已知的N个位置的PSF及其对应的位置坐标，以坐标为中心对N个位置处的PSF图像进行裁剪，得到裁剪后的PSF图像，裁剪时要保证多数像素值为零的区域被裁掉，同时又能保留PSF图像的有效区域；
②对裁剪后的PSF进行归一化处理，得到归一化后的PSF图像，记为；
③计算归一化PSF之间的协方差，得到协方差矩阵；
（1）
④对协方差矩阵进行特征值分解，得到协方差矩阵阵的N个特征值{λ 1,λ 2,…λ N},和N个特征向量{,/>,…/>}；
⑤利用特征向量构造正交基函数，也称为“特征-点扩散函数”；
（2）
其中，x ij表示第i个特征向量xi的第j个元素；表示第i个特征-点扩散函数。
⑥利用特征-点扩散函数来表征已知位置的PSF图像；
(3)
上述模型实现了对已知位置的PSF图像的拟合，并将已知的PSF图像表示为一系列特征-点扩散函数的线性组合，下一步则是估计其余位置处的PSF图像。
从上述PSF图像的拟合可以看出，要估计其余位置处的PSF图像只需要估计其余位置处特征-点扩散函数对应的系数，即对已知位置的PSF的特征-点扩散函数对应的系数矩阵进行插值，获取其余位置处的模型参数。
本实施例采用径向基函数插值法按列对系数矩阵进行插值。每次选取系数矩阵的一列并对该列系数重新排列，将其写成二维矩阵的形式，保证系数在二维矩阵中的位置与已知的PSF图像对应的位置坐标相匹配；然后插值获得其余位置处的系数值，插值后获得的二维矩阵与成像面尺寸相同，记为/>。结合拟合得到的PSF表征模型和插值获得的模型参数，就可以得到其余位置处的PSF图像，从而得到了所有位置的PSF图像。
需要说明的是，插值前的系数矩阵仅包含已知位置处PSF图像对应的模型参数，它是不完整的系数矩阵，而插值获得的系数矩阵/>既包含了已知位置处PSF图像对应的模型参数，又包含了其余未知位置处PSF图像的模型参数，是完整的系数矩阵，在实际应用中，插值的方法是可选择的，例如双三次插值、多项式拟合、反距离加权法。
通过SVPSF建模，成像面每个位置处的PSF可以写成特征-点扩散函数的线性组合，以构建空变系统PSF的表征模型，如下形式：
（4）
基于所述空变系统PSF的表征模型简化空变系统的成像模型，因而基于上式（4），空变系统的成像模型可以简化表示成以下形式：
（5）
式中，I和O分别表示退化图像和清晰图像，P表示成像系统的点扩散函数，ƞ表示噪声项，*表示卷积，表示第i个特征-点扩散函数，/>表示表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，/>表示物面的坐标，/>表示像面的坐标，表示物面和像面坐标。
S2：构建基于稀疏对数梯度正则约束的RL解卷积算法，利用所述空变系统PSF和所述稀疏对数梯度正则约束的RL解卷积算法对所采集到的光声图像进行解卷积，得到解卷积图像；
从统计角度考虑，图像复原问题可以用贝叶斯最大后验估计来建模，即选取O使得后验密度函数最大：
（6）
其等价于最小化后验密度函数的负对数：
（7）
本实施例假设似然密度函数p(I|O)服从泊松分布，即经典的RL算法：
（8）
其中，表示清晰图像的先验密度函数，·表示乘积，/>表示退化图像中坐标所在位置对应的像素值求阶乘。
由于图像复原具有病态性，需要用图像的先验模型信息进行约束，考虑到自然图像梯度的概率分布是稀疏的，即梯度值集中在零附近，根据Markov随机场理论和Hammersley-Clifford定理，本申请采用重尾分布函数对建模，即
（9）
其中μ和γ为两个常参数，在复原过程中，可根据复原情况自由决定其大小，不同的先验分布对应了不同的正则项，在本实施例中，利用重尾分布函数对先验图像建模得到的正则项称为稀疏对数梯度(SLG)正则项：
（10）
结合式（7）-（9）和最大期望算法，本实施例给出了适用于空变图像复原的基于稀疏对数梯度正则约束的RL解卷积算法（为了方便表示，以下将该解卷积算法简称为SLG-RL算法）的具体实现步骤：
（11）
（12）
（13）
（14）
（15）
需要说明的是，和/>表示傅里叶变换和傅里叶逆变换，k表示迭代次数，符号“^”表示频域形式的结果，/>表示第i个特征-点扩散函数，/>表示将/>逆序排列的结果。/>表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，I表示退化图像，O表示清晰图像，/>和/>为常参数，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代时的正则约束，/>表示利用第/>次迭代估计出的清晰图像和空变系统PSF的表征模型预测的退化图像，并在频域表示出该结果，表示第/>次迭代预测出的退化图像与真实退化图像的残差。/>表示频域形式的修正项，用来优化/>以得到/>。
在具体的实施过程中，初次迭代时，O 0=I；迭代次数视复原效果而定，一般情况下，迭代20-30次便能获得较好的复原效果。
S3：将所述解卷积图像绕图像中心旋转设定角度后，得到旋转后图像，将所述旋转后图像与所述解卷积图像相加取平均，得到消除了条纹伪影平均图像，利用所述RL解卷积算法对所述平均图像进行复原以消除旋转平均操作带来的模糊，得到最终的复原图像。
使用环形阵列换能器采集得到的光声图像存在条纹伪影，步骤S2中的SLG-RL解卷积算法虽然可以提高图像分辨率但无法消除条纹伪影，由于光声图像中的条纹伪影呈周期分布，本实施例根据条纹伪影的分布特点提出了一种消除条纹伪影的方法，需要说明的是，本方法处理的是解卷积后的图像，另外本方法仅适用于处理使用稀疏阵元的环阵换能器采集得到的光声图像；当使用的换能器阵元个数越多时，阵元排布越紧密，所获得的光声信号数据越丰富，重建得到的光声图像的质量也会越高，条纹伪影几乎不可见。
图4给出了去伪影的具体步骤，首先将解卷积图像绕图像中心按顺时针、逆时针方向分别旋转特定的角度，并将其与解卷积图像相加后取平均，如下式所示：
（16）
式中O decov表示解卷积后的图像，和/>分别表示将解卷积图像分别按顺时针、逆时针方向旋转特定角度后的图像，O ave表示平均图像。
通过该对解卷积图像进行旋转取平均，条纹伪影的峰值和谷值相互叠加并抵消，因此得到的平均图像中大部分的条纹伪影已被消除；但是，上述的旋转平均操作类似于一个退化过程，平均图像会产生轻微的模糊，因此，本申请通过MATLAB仿真得到旋转平均过程对应的PSF图像，并再次利用上述提出的SLG-RL解卷积算法对平均图像进行复原；由于大部分条纹伪影波峰波谷已经相互抵消，平均图像再次解卷积后条纹伪影不会再出现，最终得到的复原图像既消除了伪影又有效地保留了图像细节。
其中，利用MATLAB仿真旋转平均过程对应PSF图像的步骤如下：①创建与解卷积图像大小相等的零矩阵；②设定某一像素的值为1，其余像素值为0，按式（16）进行旋转平均得到该像素所在位置处对应的旋转平均PSF图像；以此类推，得到部分像素对应的PSF图像；③利用SVPSF建模所述方法获取所有像素位置对应的PSF图像。
需要说明的是，利用主成分分析分解SVPSF的方法并非本实施例的特有技术方案，目前已有多篇文献利用该方法对SVPSF建模。本实施例推导了基于稀疏对数梯度正则约束（受SLG正则约束）的RL解卷积算法，其中SLG正则项为本实施例特有的技术特征，该正则项由本实施例首次提出，此外，本实施例中根据条纹伪影的分布规律提出的去伪影方法为本实施例有的技术特征，该去伪影方法基于本实施例提出的SLG-RL解卷积算法（基于稀疏对数梯度正则约束的RL解卷积算法），简单易行；该去伪影方法在消除伪影的同时也能做到保留图像细节。
因而根据步骤S1至S3，与现有的光声图像复原算法相比较，该光声图像复原方法具有以下优势：
（一）本方法使用了SVPSF对光声图像进行复原。从理论上来说，使用SVPSF对图像复原更符合实际情况，相比使用空不变的PSF的复原算法，该方法进一步提升了图像的质量，能更好地恢复图像细节。
（二）本方法中提出的SLG正则项，其抑制噪声的能力和复原效果优于传统的Tikhonov正则项，这是因为自然图像的梯度分布是稀疏的，梯度值都趋近于0；SLG正则项是利用重尾分布函数对图像先验建模得到的，其更符合自然图像的概率模型；为了更直观地说明SLG正则项的优势，图5给出了一幅真实采集的光声图像及其水平梯度分布曲线，同时也给出了Tikhonov正则项和SLG正则项对应的位势函数曲线；从图5中可以看出，光声图像的水平梯度都集中在零值附近，SLG正则项对应的位势函数能更好地拟合真实图像的梯度分布情况。
（三）本实施例提出的去伪影方法，适用于使用环形阵列换能器采集得到的光声图像，且该方法简单易操作，去除条纹伪影的同时能较好地保留图像细节。
（四）本实施例中提出的SLG-RL解卷积算法不仅适用于光声图像复原，也可应用于其他具有SVPSF的成像系统，例如荧光成像、显微成像、天文成像等领域。
为了对本方法的技术方案、复原效果有更好的了解，以下结合实施例，对本方法的性能效果进行详细描述。以下两个实施例皆为活体动物实验：
其中设定，实验中选用了具有256个阵元的非聚焦环形阵列换能器接收光声信号，该换能器的成像半径为25mm，中心频率为7.5MHz，带宽约为73%。实验过程中，利用可调谐激光系统发射出的近红外光脉冲照亮样品以激发光声信号，并使用Verasonics系统采集光声信号，在所有的实验研究中，光声图像的重建都采用了滤波反投影算法，在以下的实施例中，去除条纹伪影过程中的旋转角度都设置为1°。
实施例1：小鼠腹部血管图像的复原
本实施例主要是对采集到的小鼠腹部血管图像进行复原。如图6所示，其中，（a）图表示采集到的光声图像，（b）图表示利用本方法中的SLG-RL算法得到的解卷积图像，（c）图表示利用本方法中的去伪影技术得到去伪影图像，即最终的复原图像。为了更直观地看出本方法的复原效果，（d）-（f）给出了（a）-（c）图实线框内的放大细节图。如图中箭头1所示，解卷积后小鼠血管变得更锐利，图像的细节得到恢复；但解卷积并不能消除图像中的条纹伪影。因此，使用去伪影技术对解卷积图像进行处理；从（f）图可以看出，使用去伪影技术后，图像中的条纹伪影大部分被抑制，如箭头2所示。
实施例2：成年男性手指图像的复原
本实施例主要是对手指横截面血管图像进行复原的结果。同样地，如图7所示，（a）图表示采集到的光声图像，（b）图表示利用本方法中的SLG-RL算法得到的解卷积图像，（c）图表示利用本方法中的去伪影技术得到去伪影图像，即最终的复原图像。（d）-（f）给出了（a）-（c）图实线框内的放大细节图。如图中箭头所示，解卷积后图像整体变得更锐利，手指血管的轮廓变得更清晰，相邻血管之间易于区分。相较于（e）图，（f）图中的条纹伪影得到了有效地抑制，如箭头所示。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (8)

1.一种基于空变点扩散函数的光声图像复原方法，包括如下步骤：
S1：获取部分不同空间位置的PSF图像，组成PSF图像样本集，基于所述PSF图像样本集构建空变系统PSF的表征模型，获得空变系统所有位置处的PSF图像；
S2：构建基于稀疏对数梯度正则约束的RL解卷积算法，利用所述空变系统PSF的表征模型和所述稀疏对数梯度正则约束的RL解卷积算法对所采集到的光声图像进行解卷积，得到解卷积图像；
S3：将所述解卷积图像绕图像中心旋转设定角度后，得到旋转后图像，将所述旋转后图像与所述解卷积图像相加取平均，得到去除了条纹伪影的平均图像，利用所述稀疏对数梯度正则约束的RL解卷积算法对所述平均图像进行复原，得到最终的复原图像；
基于稀疏对数梯度正则约束的RL解卷积算法的构建过程如下：
其中，和/>分别表示傅里叶变换和傅里叶逆变换，k表示迭代次数，符号“^”表示频域形式的结果，/>表示第i个特征-点扩散函数，/>表示将/>逆序排列的结果，/>表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，I表示退化图像，O表示清晰图像，/>和/>为常参数，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代时的正则约束，/>表示利用第/>次迭代估计出的清晰图像和空变系统PSF的表征模型预测的退化图像，并在频域表示出该结果，/>表示第/>次迭代预测出的退化图像与真实退化图像的残差，/>表示频域形式的修正项，用来优化/>以得到/>，·表示乘积。
2.根据权利要求1所述的基于空变点扩散函数的光声图像复原方法，其特征在于，在步骤S1：获取部分不同空间位置的PSF图像，组成PSF图像样本集，基于PSF图像样本集构建空变系统PSF的表征模型，获取空变系统所有位置处的PSF图像中，具体包括：
利用主成分分析拟合已知位置的PSF图像，通过一系列特征-点扩散函数的线性组合表征所述已知位置的PSF图像；
利用径向基函数插值法按列对所述特征-点扩散函数对应的系数矩阵进行插值，得到所有位置的PSF图像，基于所述空变系统PSF的表征模型简化空变系统的成像模型。
3.根据权利要求2所述的基于空变点扩散函数的光声图像复原方法，其特征在于，在所述利用主成分分析拟合已知位置的PSF图像，通过一系列特征-点扩散函数的线性组合表征所述已知位置的PSF图像中，具体包括：
记录已知的N个位置处的PSF图像及其对应的位置坐标，以坐标为中心对N个位置处的PSF图像进行裁剪，得到裁剪后的PSF图像；
对所述裁剪后的PSF图像进行归一化处理，得到归一化后的PSF图像，记为；
计算所述归一化后的PSF图像之间的协方差，得到协方差矩阵，其中，，/>；
对所述协方差矩阵进行特征值分解，得到所述协方差矩阵的N个特征值和N个特征向量/>；
利用所述特征向量构造特征-点扩散函数，通过所述特征-点扩散函数的线性组合表征所述已知位置的PSF图像。
4.根据权利要求3所述的基于空变点扩散函数的光声图像复原方法，其特征在于，基于所述空变系统PSF的表征模型构建空变系统的成像模型，所述空变系统PSF的表征模型P简化如下：
所述空变系统的成像模型简化为下式：
其中，I表示退化图像，O表示清晰图像，表示噪声项，*表示卷积，/>表示第i个特征-点扩散函数，/>表示表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，/>表示物面的坐标，/>表示像面的坐标，/>表示物面和像面坐标。
5.根据权利要求4所述的基于空变点扩散函数的光声图像复原方法，其特征在于，在步骤S2：构建基于稀疏对数梯度正则约束的RL解卷积算法中，具体公式还包括如下：
采用贝叶斯最大后验估计表述图像复原问题，选取O函数使得后验密度函数最大，等价于最小化后验密度函数的负对数；
假设似然密度函数服从泊松分布，得到经典RL算法；
根据Markov随机场理论和Hammersley-Clifford定理，采用重尾分布函数对清晰图像的先验密度函数建模得到的正则项称为稀疏对数梯度正则项；
基于所述稀疏对数梯度正则项构建稀疏对数梯度正则约束的RL解卷积算法。
6.根据权利要求5所述的基于空变点扩散函数的光声图像复原方法，其特征在于，所述O函数公式如下：
所述等价于最小化后验密度函数的负对数公式如下：
所述似然密度函数公式如下：
所述清晰图像的先验密度函数公式如下：
所述稀疏对数梯度正则项具体如下：
其中，·表示乘积，和/>为常参数，/>表示退化图像中/>坐标所在位置对应的像素值求阶乘。
7.根据权利要求1所述的基于空变点扩散函数的光声图像复原方法，其特征在于，在步骤S3：将所述解卷积图像绕图像中心旋转设定角度后，得到旋转后图像，将所述旋转后图像与所述解卷积图像相加取平均，得到去除了条纹伪影的平均图像，利用所述稀疏对数梯度正则约束的RL解卷积算法对所述平均图像进行复原以消除旋转平均操作带来的模糊，得到最终的复原图像中，具体包括：
将所述解卷积图像绕图像中心按顺时针旋转设定角度后，得到顺时针旋转图像，将所述顺时针旋转图像与所述解卷积图像相加取平均，得到第一旋转图像；
将所述解卷积图像绕图像中心按逆时针旋转设定角度后，得到逆时针旋转图像，将所述逆时针旋转图像与所述解卷积图像相加取平均，得到第二旋转图像；
将第一旋转图像和第二旋转图像相加取平均，得到去除了条纹伪影的平均图像；
利用所述稀疏对数梯度正则约束RL解卷积算法对所述平均图像进行复原以消除旋转平均操作带来的模糊，得到最终的复原图像。
8.一种基于空变点扩散函数的光声图像复原系统，其特征在于，包括构建模块、解卷积模块和去伪影模块；
构建模块用于获取部分不同空间位置的PSF图像，组成PSF图像样本集，基于PSF图像样本集构建空变系统PSF的表征模型，获得空变系统所有位置处的PSF图像；
解卷积模块用于构建基于稀疏对数梯度正则约束的RL解卷积算法，利用所述空变系统PSF的表征模型和所述稀疏对数梯度正则约束的RL解卷积算法对所采集到的光声图像进行解卷积，得到解卷积图像；
去伪影模块用于将所述解卷积图像绕图像中心旋转设定角度后，得到旋转后图像，将所述旋转后图像与所述解卷积图像相加取平均，得到去除了条纹伪影的平均图像，利用所述稀疏对数梯度正则约束的RL解卷积算法对所述平均图像进行复原，得到最终的复原图像；
基于稀疏对数梯度正则约束的RL解卷积算法的构建过程如下：
其中，和/>分别表示傅里叶变换和傅里叶逆变换，k表示迭代次数，符号“^”表示频域形式的结果，/>表示第i个特征-点扩散函数，/>表示将/>逆序排列的结果，/>表示对所述特征-点扩散函数的系数矩阵进行插值后获得的完整的系数矩阵，I表示退化图像，O表示清晰图像，/>和/>为常参数，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代估计出的清晰图像，/>表示第/>次迭代时的正则约束，/>表示利用第/>次迭代估计出的清晰图像和空变系统PSF的表征模型预测的退化图像，并在频域表示出该结果，/>表示第/>次迭代预测出的退化图像与真实退化图像的残差，/>表示频域形式的修正项，用来优化/>以得到/>，·表示乘积。
CN202310432349.1A 2023-04-21 2023-04-21 一种基于空变点扩散函数的光声图像复原方法及系统 Active CN116167948B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310432349.1A CN116167948B (zh) 2023-04-21 2023-04-21 一种基于空变点扩散函数的光声图像复原方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310432349.1A CN116167948B (zh) 2023-04-21 2023-04-21 一种基于空变点扩散函数的光声图像复原方法及系统

Publications (2)

Publication Number Publication Date
CN116167948A CN116167948A (zh) 2023-05-26
CN116167948B true CN116167948B (zh) 2023-07-18

Family

ID=86416674

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310432349.1A Active CN116167948B (zh) 2023-04-21 2023-04-21 一种基于空变点扩散函数的光声图像复原方法及系统

Country Status (1)

Country Link
CN (1) CN116167948B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117274094B (zh) * 2023-09-22 2024-05-28 哈尔滨工业大学 一种用于可穿戴超声成像质量提升的反卷积重构方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009070120A1 (en) * 2007-11-30 2009-06-04 Sidec Technologies Ab Lp-regularization of sparse representations applied to structure determination methods in molecular biology/structural chemistry
WO2019174068A1 (zh) * 2018-03-15 2019-09-19 华中科技大学 基于距离加权稀疏表达先验的图像复原与匹配一体化方法
US10502593B1 (en) * 2018-06-07 2019-12-10 Philip M. Johnson Linear and rotary multitrack absolute position encoder and methods using the same

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2574216B1 (en) * 2010-09-28 2015-11-25 Max-Planck-Gesellschaft zur Förderung der Wissenschaften e.V. Method and device for recovering a digital image from a sequence of observed digital images
CN102005031B (zh) * 2010-11-03 2012-05-02 宁波鑫高益磁材有限公司 一种mri系统中k空间采样数据的运动伪影消除方法及装置
CN102063716B (zh) * 2011-01-13 2012-07-04 耿则勋 一种基于各向异性约束的多帧迭代盲解卷积图像复原方法
CN102236887A (zh) * 2011-03-11 2011-11-09 贵州大学 基于旋转差分和加权总变分的运动模糊图像复原方法
US9307208B2 (en) * 2014-03-07 2016-04-05 The University Of British Columbia System and method for solving inverse imaging problems
US10416087B2 (en) * 2016-01-01 2019-09-17 Kla-Tencor Corporation Systems and methods for defect detection using image reconstruction
CN105976332B (zh) * 2016-05-03 2019-03-01 北京大学深圳研究生院 基于图像中亮条纹信息的图像去模糊方法
WO2018146562A1 (en) * 2017-02-09 2018-08-16 Technion Research & Development Foundation Ltd. Sparsity-based super-resolution correlation microscopy
CN107507135B (zh) * 2017-07-11 2020-04-24 天津大学 基于编码光圈和靶标的图像重构方法
CN112581378B (zh) * 2019-09-30 2022-09-13 河海大学常州校区 基于显著性强度和梯度先验的图像盲去模糊方法和装置
CN113592728B (zh) * 2021-07-01 2024-04-05 温州理工学院 一种图像复原方法、系统、处理终端及计算机介质
CN114757826A (zh) * 2022-03-24 2022-07-15 湘潭大学 一种基于多特征的pocs图像超分辨率重建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009070120A1 (en) * 2007-11-30 2009-06-04 Sidec Technologies Ab Lp-regularization of sparse representations applied to structure determination methods in molecular biology/structural chemistry
WO2019174068A1 (zh) * 2018-03-15 2019-09-19 华中科技大学 基于距离加权稀疏表达先验的图像复原与匹配一体化方法
US10502593B1 (en) * 2018-06-07 2019-12-10 Philip M. Johnson Linear and rotary multitrack absolute position encoder and methods using the same

Also Published As

Publication number Publication date
CN116167948A (zh) 2023-05-26

Similar Documents

Publication Publication Date Title
Yu et al. Solving inverse problems with piecewise linear estimators: From Gaussian mixture models to structured sparsity
Pustelnik et al. Wavelet-based image deconvolution and reconstruction
Satish et al. A Comprehensive Review of Blind Deconvolution Techniques for Image Deblurring.
US8971599B2 (en) Tomographic iterative reconstruction
Pantin et al. Deconvolution and blind deconvolution in astronomy
CN116167948B (zh) 一种基于空变点扩散函数的光声图像复原方法及系统
Shen et al. A blind restoration method for remote sensing images
Michailovich et al. Deconvolution of medical images from microscopic to whole body images
Thapa et al. Comparison of super-resolution algorithms applied to retinal images
Cascarano et al. Plug-and-Play gradient-based denoisers applied to CT image enhancement
Wang et al. Blurred image restoration using knife-edge function and optimal window Wiener filtering
González et al. Compressive optical deflectometric tomography: A constrained total-variation minimization approach
González et al. Non-parametric PSF estimation from celestial transit solar images using blind deconvolution
EP2791907A1 (fr) Procédé de reconstruction d'un signal en imagerie médicale à partir de mesures expérimentales perturbées, et dispositif d'imagerie médicale mettant en uvre ce procédé
CN112488920A (zh) 一种基于类高斯模糊核的图像正则化超分辨重建方法
Hou et al. Image deconvolution with multi-stage convex relaxation and its perceptual evaluation
Solanki et al. An efficient satellite image super resolution technique for shift-variant images using improved new edge directed interpolation
Zhang et al. Restoration of images and 3D data to higher resolution by deconvolution with sparsity regularization
Petrovici et al. Maximum entropy principle in image restoration
EP2377093A1 (en) Method and an apparatus for deconvoluting a noisy measured signal obtained from a sensor device
Sciacchitano Image reconstruction under non-Gaussian noise
CN107886478A (zh) 一种ct图像重建方法及系统、终端及可读存储介质
Kathiravan et al. An overview of sr techniques applied to images, videos and magnetic resonance images
Wang et al. L_1-Norm Regularization for Short-and-Sparse Blind Deconvolution: Point Source Separability and Region Selection
Peace Applications of Deep Learning for Ill-Posed Inverse Problems Within Optical 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