CN110610528A - 基于模型的双约束光声断层图像重建方法 - Google Patents
基于模型的双约束光声断层图像重建方法 Download PDFInfo
- Publication number
- CN110610528A CN110610528A CN201910758438.9A CN201910758438A CN110610528A CN 110610528 A CN110610528 A CN 110610528A CN 201910758438 A CN201910758438 A CN 201910758438A CN 110610528 A CN110610528 A CN 110610528A
- Authority
- CN
- China
- Prior art keywords
- image
- regularization
- model
- reconstruction
- nlm
- 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.)
- Withdrawn
Links
- 238000000034 method Methods 0.000 title claims abstract description 78
- 238000003325 tomography Methods 0.000 title claims abstract description 25
- 239000011159 matrix material Substances 0.000 claims abstract description 13
- 239000013598 vector Substances 0.000 claims abstract description 6
- 230000006870 function Effects 0.000 claims description 46
- 238000012417 linear regression Methods 0.000 claims description 11
- 238000004458 analytical method Methods 0.000 claims description 9
- 238000006467 substitution reaction Methods 0.000 claims description 6
- 238000011478 gradient descent method Methods 0.000 claims description 3
- 230000002401 inhibitory effect Effects 0.000 abstract description 4
- 239000000126 substance Substances 0.000 abstract description 2
- 238000004422 calculation algorithm Methods 0.000 description 10
- 238000003384 imaging method Methods 0.000 description 10
- 238000002474 experimental method Methods 0.000 description 6
- 210000001519 tissue Anatomy 0.000 description 5
- 238000001514 detection method Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 238000001727 in vivo Methods 0.000 description 2
- 210000003734 kidney Anatomy 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 210000004204 blood vessel Anatomy 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 239000002872 contrast media Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 210000003128 head Anatomy 0.000 description 1
- 210000004185 liver Anatomy 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 239000012528 membrane Substances 0.000 description 1
- 238000011580 nude mouse model Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000012634 optical imaging Methods 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 210000000813 small intestine Anatomy 0.000 description 1
- 210000000278 spinal cord Anatomy 0.000 description 1
- 210000000952 spleen Anatomy 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000002604 ultrasonography Methods 0.000 description 1
- 238000012285 ultrasound imaging Methods 0.000 description 1
Classifications
-
- 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
- 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
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/90—Dynamic range modification of images or parts thereof
- G06T5/92—Dynamic range modification of images or parts thereof based on global image properties
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/40—Analysis of texture
-
- 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
- G06T2207/10101—Optical tomography; Optical coherence tomography [OCT]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种基于模型的双约束光声断层图像重建方法，包括如下步骤：S1、引入λnlm和λs对最小化代价函数进行优化得到目标函数其中，为最小二乘项，p是大小为的列向量表示投影数据，W的大小为表示系统矩阵，R(x)为非局部均值的正则化项；||D‑1x||1为稀疏编码约束的正则化项，λnlm为非局部均值的正则化参数，λs为稀疏编码约束的正则化参数；S2、更新最小二乘项和非局部均值的正则化项；S3、引入中间变量u对稀疏编码约束的正则化项进行优化；S4、引入软阈值对稀疏编码约束的正则化项进行更新；S5、通过迭代进行图像重建；在抑制光声图像中噪声和伪影的同时，很好的保持组织的边缘、纹理等结构信息。
Description
技术领域
本发明涉及数字图像处理技术领域，特别涉及一种基于模型的双约束光声断层图像重建方法。
背景技术
光声断层成像(OAT)是一种新型的混合成像模式，它结合了光学成像高对比度和超声成像组织探测深度的优势。光声断层成像丰富的内源性和外源性造影剂使得该成像技术可用于非侵入性的解剖结构成像、功能成像和分子成像，由于这些独特的优势，光声断层成像被广泛的应用到临床前医学和生物医学研究中。图像重建作为光声成像技术的关键步骤，对光声图像的质量具有很大的影响。光声图像的重建是利用一组实测超声信号恢复物体的初始声压分布，是一个病态的逆问题。重建图像的精度受有限的探测器角度、系统矩阵和信号探测过程中的噪声等因素的影响。因此，为了提高图像重建的质量，可以在重建过程中引入适当的惩罚或先验知识。
通常，光声图像的重建算法可以被分为两类：一类是解析重建算法，如反投影算法；另一类是基于模型的迭代重建算法。对于解析重建算法，王立红研究团队首先提出了一种通用的反投影重建算法(UBP)用于光声图像的重建，该方法适用于平面探测器、柱面探测器和球面探测器三种探测器类型，算法简单，计算速度快，但是由于该方法是基于对声波传播和检测的理想描述，因此难以推广到实际的光声模型当中。此外，该方法易在待重建的图像中引入条形伪影，限制了其在功能成像和分子成像方面的应用。
基于模型的重建方法不是通过求解光声方程的解析解来重建图像，而是通过最小化测量得到的超声信号与成像模型预测到的声信号之间的误差来重建光声图像。基于模型的重建算法通过迭代寻求逆问题的全局最优解，可以抑制传统反投影算法在待重建图像中引入的伪影，如投影数据不完整、低频信息丢失等原因在待重建图像中产生的伪影。基于模型的重建方法的另一个优势是该方法可以将成像场景(包括成像几何和目标对象)的先验知识作为正则化项纳入重建过程，生成更高质量的图像。
虽然基于模型的算法具有上述优点，但是如果模型过于近似，不能保持成像装置固有的物理性质，那么重建的图像就不能反映光学吸收的真实分布。因此，需要运用正则化的方法来求解这个病态逆问题。
针对光声图像重建的不适定问题，研究者提出了多种正则化方法，包括Tikhonov方法、全变分方法和基于稀疏性的方法。在这些方法中，大多数的先验知识都是认为待重建的目标图像应该是平滑的、孤立的、稀疏分布的。然而，这些方法通常会由于原始信号的过平滑而导致细节丢失。
因此，针对现有技术的不足，本发明提出了一种基于模型的双约束光声断层图像重建方法以克服现有技术不足甚为必要。
发明内容
本发明的目的是为了解决现有技术中的缺陷，提供一种基于模型的双约束光声断层图像重建方法，结合了非局部均值约束的去噪特性、边缘保持特性和稀疏编码的稀疏特性，在有效的抑制噪声和伪影的同时，保留了待重建图像的边缘、纹理等结构信息。
为实现上述目的，本发明的技术方案如下：
一种基于模型的双约束光声断层图像重建方法，包括如下步骤：
S1、引入非局部均值的正则化参数和稀疏编码约束的正则化参数对最小化代价函数进行优化得到目标函数，
其中，x为重建后的图像，为最小二乘项，p是大小为的列向量表示投影数据，W的大小为表示系统矩阵，R(x)为非局部均值的正则化项，||D-1x||1为稀疏编码约束的正则化项，λnlm为非局部均值的正则化参数，λs为稀疏编码约束的正则化参数；
S2、基于梯度下降和迭代软阈值化进行图像重建。
为了进一步实现本发明，所述S2具体包括：
S2.1、基于梯度下降法对对目标函数中的最小二乘项和非局部均值的正则化项进行更新得到待重建图像的中间值；
S2.2、引入中间变量u对稀疏编码约束的正则化项进行优化并得到关于中间变量的目标函数；
S2.3、引入软阈值对稀疏编码约束的正则化项进行更新得到更新值；
S2.4、将更新值转移回S2.2中的关于中间变量的目标函数，得到重建后的图像。
为了进一步实现本发明，所述S2.1具体包括有：
S2.1.1、定义子目标函数L(x)：
S2.1.2、对子目标函数L(x)进行迭代，在每次迭代中，将xk+1更新为：
其中，τ是迭代步长；WT中T是矩阵W的转置符号；I表示图像像素的索引；i和j分别表示第i个和j个像素点，Si是以第i个像素为中心的搜索窗的大小。
为了进一步实现本发明，所述S2.2的具体步骤包括引入中间变量u，令x＝Du，可以得到关于的u的另一个目标函数P(u)：
为了进一步实现本发明，所述S2.3的具体步骤包括引入软阈值操作：
ρ>0，将uk+1更新为：
为了进一步实现本发明，所述S2.4的具体步骤包括使用矩阵D进行逆变操作将uk+1转换成图像域，得到更新后的图像：
xk+1＝Duk+1. 公式(VII)
为了进一步实现本发明，所述S2.1.2的具体步骤包括：
S2.1.2.1、定义总迭代次数为K，当前迭代次数为k，且0≤k≤K；
S2.1.2.2、令k＝1；
S2.1.2.3、将采集到的参数代入公式(III)得到待重建图像的中间值x；
S2.1.2.4、判断第k次迭代时重建图像的中间值x的误差是否小于ε或k是否小于对大迭代次数K，若是，则输出k次迭代的待重建图像，否则，则进入S2.1.2.5；
S2.1.2.5、令k＝k+1，返回S2.1.2.3。
为了进一步实现本发明，所述目标函数中λnlm的数值确认如下：
S(1)、定义λs＝0，将目标函数转换为单正则化函数；
S(2)、采用试错法确定λnlm的数值；
所述S(2)中试错法具体步骤如下：
S(2.1)、建立模型λnlm＝β0+μβ1+σ2β2+(μ·σ2)β3，其中，μ为无正则化图像的均值，σ2为无正则化图像的方差，β0、β1、β2、β3均为调整图像重建学习系数；
S(2.2)、采集待调整图像光生投影数据，使用两种方法对采集到的光声投影数据进行图像重构；
一种是经验确定的正则化参数重构，将得到一组由经验确定的λnlm；
另一种是没有正则化重构，得到一组μ和σ2；
S(2.3)、建立一个线性系统λnlm＝xβ，x＝[ones(size(μ)，1)，μ,σ2)；
S(2.4)、通过多元线性回归分析估计线性系统的系数β0、β1、β2、β3，从而得出λnlm；
S(2.5)、将λnlm代入求p和R2的值，其中R为决定系数；
S(2.6)、若S(2.5)中的p小于设定值且R2-1的绝对值小于等于设定值时，则该λnlm正确，否则，进入S(2.2)。
为了进一步实现本发明，所述目标函数中的λs数值确认如下：
S(a)、定义λnlm＝0，将目标函数转换为单正则化函数；
S(b)、采用试错法确定λs的数值。
优选的，上述的基于模型的双约束光声断层图像重建方法，所述S(b)中试错法具体步骤如下：
S(b.1)建立模型λs＝β0+μβ1+σ2β2+(μ·σ2)β3，其中，μ为无正则化图像的均值，σ2为无正则化图像的方差，β0、β1、β2、β3均为调整图像重建学习系数；
S(b.2)获得一系列图像，进行两种重构；
一种是经验确定的正则化参数重构，将得到一组由经验确定的λs；
另一种是没有正则化重构，得到一组μ和σ2；
S(b.3)建立一个线性系统λs＝xβ，x＝[ones(size(μ)，1)，μ,σ2)；
S(b.4)通过多元线性回归分析估计线性系统的系数β0、β1、β2、β3，从而得出λs；
S(b.5)将λs代入求p和R2的值，其中R为决定系数；
S(b.6)若步骤(b.5)中的p小于设定值且R2-1的绝对值小于等于设定值时，则该λs正确，否则，进入步骤(b.2)。
为了进一步实现本发明，所述设定值为0.01或0.001或0.0001。
有益效果
本发明一种基于模型的双约束光声断层图像重建方法采用的双先验约束方法结合了非局部均值先验和稀疏先验的优势，在抑制光声图像中噪声和伪影的同时，很好的保持组织的边缘、纹理等结构信息。该方法形式简洁，不需要对待重建图像进行预处理即可达到很好的去噪效果。对于正则化参数的选择，应用了待重建图像的统计特性，通过建立多元线性回归分析模型，确定正则化参数和待重建图像的均值、方差之间的关系，具有一定的实用性。
附图说明
图1为本发明一种基于模型的双约束光声断层图像重建方法中测试数据集的多元线性回归图；
图2为实施例2仿真实验的重建结果图；
图3为实施例2体膜实验的重建结果图；
图4为实施例2活体实验的重建结果图。
具体实施方式
下面结合附图对本发明作进一步地详细的说明，这些附图均为简化的示意图，仅以示意方式说明本发明的基本结构。
实施例1
光声成像的前向传播模型可以表示为：
p＝Wx. (1)
其中，x是大小为的列向量表示待成像的物体；p是大小为的列向量表示投影数据；W的大小为的系统矩阵。对式(1)的反演可以通过求解以下最小化问题来实现：
其中，表示二范数。这样，图像重建过程就是将实际测量的信号与模型预测的信号之间的差值减到最小。一般情况下，光声图像的重建过程是病态的，即求得的式(2)的解可能是不正确的。因此，需要引入适当的先验或惩罚来改善重建结果。这通常涉及到将一些先验知识以惩罚函数的形式合并到式(2)中，惩罚函数控制未知确定性的某些期望属性或未知随机变量的先验概率分布。
本发明提供一种基于模型的双约束光声断层图像重建方法，该图像重建方法结合非局部均值和基于一范数稀疏编码惩罚的正则化方案，双约束方案的目的是寻找一种既能抑制条纹伪影又能保留组织细节的方法，寻求的最小化的代价函数包含三个部分——数据拟合项、基于非局部均值的惩罚项和基于稀疏性的惩罚项。
具体包括如下步骤：
S1、引入非局部均值的正则化参数和稀疏编码约束的正则化参数对最小化代价函数进行优化得到目标函数，
其中，x为重建后的图像，为最小二乘项，p是大小为的列向量表示投影数据，W的大小为表示系统矩阵，R(x)为非局部均值的正则化项，||D-1x||1为稀疏编码约束的正则化项，λnlm为非局部均值的正则化参数，λs为稀疏编码约束的正则化参数；
S2、基于梯度下降和迭代软阈值化进行图像重建。
S2具体包括：
S2.1、基于梯度下降法对对目标函数中的最小二乘项和非局部均值的正则化项进行更新得到待重建图像的中间值；
S2.1.1、定义子目标函数L(x)：
S2.1.2、对子目标函数L(x)进行迭代，在每次迭代中，将xk+1更新为：
其中，τ是迭代步长；WT中T是矩阵W的转置符号；I表示图像像素的索引；i和j分别表示第i个和j个像素点，Si是以第i个像素为中心的搜索窗的大小；
S2.2、引入中间变量u对稀疏编码约束的正则化项进行优化并得到关于中间变量的目标函数；
S2.3、引入软阈值对稀疏编码约束的正则化项进行更新得到更新值；
S2.4、将更新值转移回S2.2中的关于中间变量的目标函数，得到重建后的图像。
具体的，S2.1具体包括有：
S2.1.1、定义子目标函数L(x)：
S2.1.2、对子目标函数L(x)进行迭代，在每次迭代中，将xk+1更新为：
其中，τ是迭代步长；WT中T是矩阵W的转置符号；I表示图像像素的索引；i和j分别表示第i个和j个像素点，Si是以第i个像素为中心的搜索窗的大小。
具体的，S2.2的具体步骤包括引入中间变量u，令x＝Du，可以得到关于的u的另一个目标函数P(u)：
具体的，S2.3的具体步骤包括引入软阈值操作：
ρ>0，将uk+1更新为：
具体的，S2.4的具体步骤包括使用矩阵D进行逆变操作将uk+1转换成图像域，得到更新后的图像：
xk+1＝Duk+1. 公式(VII)
更进一步的，S2.1.2的具体步骤包括：
S2.1.2.1、定义总迭代次数为K，当前迭代次数为k，且0≤k≤K；
S2.1.2.2、令k＝1；
S2.1.2.3、将采集到的参数代入公式(III)得到待重建图像的中间值x；
S2.1.2.4、判断第k次迭代时重建图像的中间值x的误差是否小于ε或k是否小于对大迭代次数K，若是，则输出k次迭代的待重建图像，否则，则进入S2.1.2.5；
S2.1.2.5、令k＝k+1，返回S2.1.2.3。
ε是一个经验值，通常取非常小的值，在本实施例中ε＝1×10-12，需要说明的是，ε的取值可以根据实际经验进行取值。
目标函数中λnlm的数值确认如下：
S(1)、定义λs＝0，将目标函数转换为单正则化函数；
S(2)、采用试错法确定λnlm的数值；
S(2)中试错法具体步骤如下：
S(2.1)、建立模型λnlm＝β0+μβ1+σ2β2+(μ·σ2)β3，其中，μ为无正则化图像的均值，σ2为无正则化图像的方差，β0、β1、β2、β3均为调整图像重建学习系数；
S(2.2)、采集待调整图像光生投影数据，使用两种方法对采集到的光声投影数据进行图像重构；
一种是经验确定的正则化参数重构，将得到一组由经验确定的λnlm；
另一种是没有正则化重构，得到一组μ和σ2；
S(2.3)、建立一个线性系统λnlm＝xβ，x＝[ones(size(μ)，1)，μ,σ2)；
S(2.4)、通过多元线性回归分析估计线性系统的系数β0、β1、β2、β3，从而得出λnlm；
S(2.5)、将λnlm代入求p和R2的值，其中R为决定系数；
S(2.6)、若S(2.5)中的p小于设定值且R2-1的绝对值小于设定值，则该λnlm正确，否则，进入S(2.2)。
设定值为0.01或0.001或0.0001，这里的R2是一个具体的数值，如果R2＝1，则认为模型是完全正确的，这是理想情况；实际中R2是不可能等于1的，这时，R2越接近1就说明模型越准确。
目标函数中的λs数值确认如下：
S(a)、定义λnlm＝0，将目标函数转换为单正则化函数；
S(b)、采用试错法确定λs的数值。
S(b)中试错法具体步骤如下：
S(b.1)建立模型λs＝β0+μβ1+σ2β2+(μ·σ2)β3，其中，μ为无正则化图像的均值，σ2为无正则化图像的方差，β0、β1、β2、β3均为调整图像重建学习系数；
S(b.2)获得一系列图像，进行两种重构；
一种是经验确定的正则化参数重构，将得到一组由经验确定的λs；
另一种是没有正则化重构，得到一组μ和σ2；
S(b.3)建立一个线性系统λs＝xβ，x＝[ones(size(μ)，1)，μ,σ2)；
S(b.4)通过多元线性回归分析估计线性系统的系数β0、β1、β2、β3，从而得出λs；
S(b.5)将λs代入求p和R2的值，其中R为决定系数；
S(b.6)若步骤(b.5)中的p小于设定值且R2-1的绝对值小于等于设定值，则该λs正确，否则，进入步骤(b.2)。
设定值为0.01或0.001或0.0001，这里的R2是一个具体的数值，如果R2＝1，则认为模型是完全正确的，这是理想情况；实际中R2是不可能等于1的，这时，R2越接近1就说明模型越准确。
为了说明这种方法，我们获取了30张活体裸鼠头部、肝脏和肾脏的横断面图像(每个切片10帧)。执行多元线性回归分析后，我们可以获得最初的系数β，这里，β0＝3.3957×10-5，β1＝-3.2573×10-7，结果β2＝-1.2376×10–7，β3＝1.4102×10-9。这些估计得到的系数的置信区间分别为：[1.9971×10-5,4.7943×10-5],[-4.6644×10-7,-1.8502×10-7],[-1.7008×10-7,-7.7443×10-8],[9.9507×10-10,1.9052×10-9]，因此，我们的线性回归方程可以写成：
λ＝3.3957×10-5-3.2573×10-7μ
-1.2376×10-7σ2+1.4102×10-9μσ2 (11)
回归结果如图1所示。为了评估模型的准确性，我们计算了p值和决定系数R2。p＝3.4762×10-14和R2＝0.9325。由于得到p<0.001并且R2非常接近1，因此式(11)中的模型可以认为是准确的。
本发明中采用的双先验约束方法结合了非局部均值先验和稀疏先验的优势，在抑制光声图像中噪声和伪影的同时，很好的保持组织的边缘、纹理等结构信息。该方法形式简洁，不需要对待重建图像进行预处理即可达到很好的去噪效果。对于正则化参数的选择，应用了待重建图像的统计特性，通过建立多元线性回归分析模型，确定正则化参数和待重建图像的均值、方差之间的关系，具有一定的实用性。
实施例2
如图2所示，在仿真实验和活体动物实验中，比较了采用双先验约束方法和其他四种现有方法在重建光声图像时性能的优劣，在图2中，(a)为原始图像；(b)为采用无先验重建出的图像，(c)为采用吉洪诺夫先验重建出的图像，(d)为采用非局部均值先验重建出的图像，(e)为采用基于小波的稀疏编码先验重建出的图像，(f)为采用本发明所提出的双先验方法重建出的图像，从图2中可以看出双先验约束方法增强了待重建图像的对比度。
计算仿真实验重建结果的均方根误差如表1所示
表1仿真实验重建结果的均方根误差
由表1可知，使用本发明方法(NLM+S先验)得出的均方根误差最小。
为了进一步验证本发明所提出的方法的优越性，图3展示了体膜实验的重建结果，图3中，(a)为原始图像；(b)为采用无先验重建出的图像，(c)为采用吉洪诺夫先验重建出的图像，(d)为采用非局部均值先验重建出的图像，(e)为采用基于小波的稀疏编码先验重建出的图像，将(b)中黑色虚线框标记的区域和白色虚线框标记的区域分别作为计算信噪比和压缩比的目标区域和背景区域。(f)为放大(a)中白色实线框区域的图像，结果表明我们的方法可以更有效地抑制条纹伪影。
图4展示了活体动物实验的重建结果(肾脏部位)，图4中，(a)为原始图像；(b)为采用无先验重建出的图像，(c)为采用吉洪诺夫先验重建出的图像，(d)为采用非局部均值先验重建出的图像，(e)为采用基于小波的稀疏编码先验重建出的图像；(f)放大(a)中标注黑框和白框区域的视图。放大后的区域显示我们的方法在本文所述的所有方法中表现最好，双约束方法重建出的图像更加平滑，细节信息保存完好且图像具有更好的组织轮廓，如图4(e)中的白色箭头所示：肾皮质(箭头1)，脾(箭头2)，小肠(箭头3)，脊髓(箭头4)，鼠标箭头(5)，以及一些细小的血管(箭头6和7)。
以上所述仅为本发明的较佳实施方式，本发明并不局限于上述实施方式，在实施过程中可能存在局部微小的结构改动，如果对本发明的各种改动或变型不脱离本发明的精神和范围，且属于本发明的权利要求和等同技术范围之内，则本发明也意图包含这些改动和变型。
Claims (11)
1.一种基于模型的双约束光声断层图像重建方法，其特征在于，包括如下步骤：
S1、引入非局部均值的正则化参数和稀疏编码约束的正则化参数对最小化代价函数进行优化得到目标函数，
其中，x为重建后的图像，为最小二乘项，p是大小为的列向量表示投影数据，W的大小为表示系统矩阵，R(x)为非局部均值的正则化项，||D-1x||1为稀疏编码约束的正则化项，λnlm为非局部均值的正则化参数，λs为稀疏编码约束的正则化参数；
S2、基于梯度下降和迭代软阈值化进行图像重建。
2.根据权利要求1所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述S2具体包括：
S2.1、基于梯度下降法对目标函数中的最小二乘项和非局部均值的正则化项进行更新得到待重建图像的中间值；
S2.2、引入中间变量u对稀疏编码约束的正则化项进行优化并得到关于中间变量的目标函数；
S2.3、引入软阈值对稀疏编码约束的正则化项进行更新得到更新值；
S2.4、将更新值转移回S2.2中的关于中间变量的目标函数，得到重建后的图像。
3.根据权利要求1所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述S2.1具体包括有：
S2.1.1、定义子目标函数L(x)：
S2.1.2、对子目标函数L(x)进行迭代，在每次迭代中，将xk+1更新为：
其中，τ是迭代步长；WT中T是矩阵W的转置符号；I表示图像像素的索引；i和j分别表示第i个和j个像素点，Si是以第i个像素为中心的搜索窗的大小。
4.根据权利要求3所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述S2.2的具体步骤包括引入中间变量u，令x＝Du，得到关于的u的另一个目标函数P(u)：
5.根据权利要求4所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述S2.3的具体步骤包括引入软阈值操作：
ρ>0，将uk+1更新为：
6.根据权利要求4所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述S2.4的具体步骤包括使用矩阵D进行逆变操作将uk+1转换成图像域，得到更新后的图像：
xk+1＝Duk+1. 公式(VII)
7.根据权利要求5所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述S2.1.2的具体步骤包括：
S2.1.2.1、定义总迭代次数为K，当前迭代次数为k，且0≤k≤K；
S2.1.2.2、令k＝1；
S2.1.2.3、将采集到的参数代入公式(III)得到待重建图像的中间值x；
S2.1.2.4、判断第k次迭代时重建图像的中间值x的误差是否小于ε或k是否小于对大迭代次数K，若是，则输出k次迭代的待重建图像，否则，则进入S2.1.2.5；
S2.1.2.6、令k＝k+1，返回S2.1.2.3。
8.根据权利要求1所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述目标函数中λnlm的数值确认如下：
S(1)、定义λs＝0，将目标函数转换为单正则化函数；
S(2)、采用试错法确定λnlm的数值；
所述S(2)中试错法具体步骤如下：
S(2.1)、建立模型λnlm＝β0+μβ1+σ2β2+(μ·σ2)β3，其中，μ为无正则化图像的均值，σ2为无正则化图像的方差，β0、β1、β2、β3均为调整图像重建学习系数；
S(2.2)、采集待调整图像光生投影数据，使用两种方法对采集到的光声投影数据进行图像重构；
一种是经验确定的正则化参数重构，将得到一组由经验确定的λnlm；
另一种是没有正则化重构，得到一组μ和σ2；
S(2.3)、建立一个线性系统λnlm＝xβ，x＝[ones(size(μ)，1)，μ,σ2)；
S(2.4)、通过多元线性回归分析估计线性系统的系数β0、β1、β2、β3，从而得出λnlm；
S(2.5)、将λnlm代入求p和R2的值，其中R为决定系数；
S(2.6)、若S(2.5)中的p小于设定值且R2-1的绝对值小于等于设定值时，则该λnlm正确，否则，进入S(2.2)。
9.根据权利要求1所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述目标函数中的λs数值确认如下：
S(a)、定义λnlm＝0，将目标函数转换为单正则化函数；
S(b)、采用试错法确定λs的数值。
10.根据权利要求8所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述S(b)中试错法具体步骤如下：
S(b.1)建立模型λs＝β0+μβ1+σ2β2+(μ·σ2)β3，其中，μ为无正则化图像的均值，σ2为无正则化图像的方差，β0、β1、β2、β3均为调整图像重建学习系数；
S(b.2)获得一系列图像，进行两种重构；
一种是经验确定的正则化参数重构，将得到一组由经验确定的λs；
另一种是没有正则化重构，得到一组μ和σ2；
S(b.3)建立一个线性系统λs＝xβ，x＝[ones(size(μ)，1)，μ,σ2)；
S(b.4)通过多元线性回归分析估计线性系统的系数β0、β1、β2、β3，从而得出λs；
S(b.5)将λs代入求p和R2的值，其中R为决定系数；
S(b.6)若步骤(b.5)中的p小于设定值且R2-1的绝对值小于等于设定值时，则该λs正确，否则，进入步骤(b.2)。
11.根据权利要求7或9所述的基于模型的双约束光声断层图像重建方法，其特征在于，所述设定值为0.01或0.001或0.0001。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910479206X | 2019-05-31 | ||
| CN201910479206 | 2019-05-31 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN110610528A true CN110610528A (zh) | 2019-12-24 |
Family
ID=68890489
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910758438.9A Withdrawn CN110610528A (zh) | 2019-05-31 | 2019-08-16 | 基于模型的双约束光声断层图像重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110610528A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111508076A (zh) * | 2020-04-16 | 2020-08-07 | 中国科学院自动化研究所 | 基于自适应参数搜索的弹性网络激发荧光断层重建系统 |
| CN111833412A (zh) * | 2020-07-16 | 2020-10-27 | 中北大学 | 一种基于分数滤波框架的Tikhonov正则化图像重建方法 |
| CN111481172B (zh) * | 2020-04-13 | 2021-08-31 | 南方医科大学 | 一种交错稀疏采样多光谱光声断层成像系统及方法 |
| CN113724351A (zh) * | 2021-08-24 | 2021-11-30 | 南方医科大学 | 一种光声图像衰减校正方法 |
-
2019
- 2019-08-16 CN CN201910758438.9A patent/CN110610528A/zh not_active Withdrawn
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111481172B (zh) * | 2020-04-13 | 2021-08-31 | 南方医科大学 | 一种交错稀疏采样多光谱光声断层成像系统及方法 |
| CN111508076A (zh) * | 2020-04-16 | 2020-08-07 | 中国科学院自动化研究所 | 基于自适应参数搜索的弹性网络激发荧光断层重建系统 |
| CN111508076B (zh) * | 2020-04-16 | 2022-04-19 | 中国科学院自动化研究所 | 基于自适应参数搜索的弹性网络激发荧光断层重建系统 |
| CN111833412A (zh) * | 2020-07-16 | 2020-10-27 | 中北大学 | 一种基于分数滤波框架的Tikhonov正则化图像重建方法 |
| CN111833412B (zh) * | 2020-07-16 | 2023-09-22 | 中北大学 | 一种基于分数滤波框架的Tikhonov正则化图像重建方法 |
| CN113724351A (zh) * | 2021-08-24 | 2021-11-30 | 南方医科大学 | 一种光声图像衰减校正方法 |
| CN113724351B (zh) * | 2021-08-24 | 2023-12-01 | 南方医科大学 | 一种光声图像衰减校正方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110610528A (zh) | 基于模型的双约束光声断层图像重建方法 | |
| CN110097512B (zh) | 基于Wasserstein生成对抗网络的MRI图像去噪模型的构建方法及应用 | |
| Yoon et al. | Motion adaptive patch-based low-rank approach for compressed sensing cardiac cine MRI | |
| Dolui et al. | A new similarity measure for non-local means filtering of MRI images | |
| CN108765313B (zh) | 基于类内低秩结构表示的高光谱图像去噪方法 | |
| Li et al. | A Douglas–Rachford splitting approach to compressed sensing image recovery using low-rank regularization | |
| CN113538260B (zh) | 自监督和有监督联合训练的ldct图像去噪与分类方法 | |
| Wang et al. | A conditional adversarial network for single plane wave beamforming | |
| CN114140442A (zh) | 一种基于频域和图像域退化感知的深度学习稀疏角度ct重建方法 | |
| Yancheng et al. | RED-MAM: A residual encoder-decoder network based on multi-attention fusion for ultrasound image denoising | |
| Jiang et al. | A new nonlocal means based framework for mixed noise removal | |
| Guan et al. | Generative modeling in sinogram domain for sparse-view CT reconstruction | |
| Lee et al. | Speckle reduction via deep content-aware image prior for precise breast tumor segmentation in an ultrasound image | |
| Chan et al. | An attention-based deep convolutional neural network for ultra-sparse-view CT reconstruction | |
| CN111340699B (zh) | 基于非局部先验和稀疏表示的磁共振图像去噪方法及装置 | |
| CN111815527A (zh) | 基于Weibull分布的混合高阶变分超声图像去噪方法 | |
| Zhou et al. | High-resolution hierarchical adversarial learning for OCT speckle noise reduction | |
| Diwakar et al. | Blind noise estimation-based CT image denoising in tetrolet domain | |
| Javed et al. | Developing a bio-inspired multi-gene genetic programming based intelligent estimator to reduce speckle noise from ultrasound images | |
| CN111311531A (zh) | 图像增强方法、装置、控制台设备及医学成像系统 | |
| CN116342444A (zh) | 一种双通道多模态图像融合方法及融合成像终端机 | |
| CN112270725A (zh) | 一种光谱层析成像中的图像重建和编码方法 | |
| Beevi et al. | Denoising transthoracic echocardiographic images in regional wall motion abnormality using deep learning techniques | |
| Lakshmi et al. | A re-constructive algorithm to improve image recovery in compressed sensing. | |
| Dey et al. | Score-based Diffusion Models for Photoacoustic Tomography Image Reconstruction |
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 | ||
| WW01 | Invention patent application withdrawn after publication |
Application publication date: 20191224 |
|
| WW01 | Invention patent application withdrawn after publication |