1.本发明属于磁
粒子成像领域,具体涉及一种基于极大极小凹与全变分约束的磁粒子成像图像重建方法、系统、设备。
背景技术:
2.磁粒子成像技术(magnetic particle imaging,mpi)是一种新兴的断层成像技术,它以高时空分辨率对磁性纳米粒子进行在体定量成像。此外,mpi只对磁性纳米粒子进行特异性成像,没有背景
信号干扰,因此,采集到的信号强度与磁性纳米粒子浓度成正比,在提供高灵敏度的同时能够进行磁粒子浓度的量化,这使得mpi在生物医学应用中成为一种非常有前景的成像方式。
3.mpi图像重建主要包括两种方法:基于x-space的重建方法与基于系统矩阵的重建方法。基于x-space的重建方法将mpi成像系统近似为一个线性移不变系统,将系统测量到的时域信号转化为到图像域,再经过网格化得到目标区域内磁粒子的分布,可实现mpi快速成像。但基于x-space的重建理论是在一系列理想假设的基础上推导而来,不能完全等价与真实的系统与磁粒子的物理特性,因此图像适量稍差。在基于系统矩阵的重建方法中,系统矩阵的构建是基于实际测量得到的,包含了复杂的系统特性与磁粒子的物理特性,重建的结果能够更加精准,因此本发明拟采用基于系统矩阵的重建方法进行重建。
4.由于mpi系统矩阵是一个病态矩阵,因此在求解中需要使用合适的正则化来稳定解。目前常使用的两种正则化方法为l2正则化与结合了全变分(total variation,tv)和l1的融合拉索正则化。其中l2正则化结合kaczmarz行迭代方法能够进行快速求解,但是会造成重建结果的过平滑,而结合了融合拉索正则化能够重建出更清晰的边缘,但是会低估图像的幅值,不利于准确定量。因此需要一种合适的正则项,既能重建出清晰的边缘的同时也能保持重建结果的幅值。基于此,本发明提出了一种基于极大极小凹与全变分约束的磁粒子成像图像重建方法。
技术实现要素:
5.为了解决现有技术中的上述问题,即为了解决现有的磁粒子成像图像重建方法重建的mpi图像误差大、精度差、质量低的问题,本发明提出了一种基于极大极小凹与全变分约束的磁粒子成像图像重建方法,该方法包括:
6.s100,获取磁粒子成像系统中待成像重建物体的仿体模型在目标区域内逐像素移动测量后采集的时域电压信号,作为输入信号;
7.s200,将
所述输入信号进行傅里叶变换,并基于傅里叶变换后的输入信号构建系统矩阵;
8.s300,基于所述系统矩阵,根据磁粒子成像原理,构建磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射
方程,作为第一方程;
9.s400,结合各向同性全变分与极大极小凹约束,对所述第一方程进行重构,得到第
二方程;
10.s500,通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,进而进行图像重建。
11.在一些优选的实施方式中,磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程为:
12.sc=u
13.其中,s表示系统矩阵,c表示磁粒子成像的浓度,u表示傅里叶变换后的输入信号。
14.在一些优选的实施方式中,所述第二方程为:
15.min
c≥0
λ
tv
‖c‖
tv
+λ
mc
‖c‖
mc
16.s.t.‖sc-u‖2≤ε
17.其中,‖
·
‖
tv
为各向同性全变分正则项,‖
·
‖
mc
为极大极小凹正则项,λ
mc
与λ
tv
为权重系数,‖sc-u‖2为数据保真项,ε≥0,ε是与噪声大小有关的约束项。
18.在一些优选的实施方式中,通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,其方法为:
19.使用增加混合惩罚函数的交替方向乘子法对所述第二方程进行求解,并引入中间变量z=[z
(0)
,z
(1)
,z
(2)
]h、g=[s
h i i]h;
[0020]
将所述第二方程转化为:
[0021]
min
x f1(c)+f2(z)
[0022]
s.t.pc+qz-s=0
[0023]
其中,p=g,q=-i,s=0,f1(c)=0,f2(z)=ι
e(ε,i,u)
(z
(0)
)+λ
tv
‖c‖
tv
+λ
mc
‖c‖
mc
,[
·
]h为共轭转置,ι
e(ε,i,u)
(z0)表示数据保真项的指示函数,i表示单位矩阵;
[0024]
基于转化后的第二方程,结合拉格朗日乘子构造增广的拉格朗日函数:
[0025][0026]
其中,d为拉格朗日乘子,β为惩罚系数,m为频点数;
[0027]
进而得到交替方向乘子法迭代的步骤如下:
[0028]
a1:
[0029]
a2:
[0030]
a3:d
k+1
=d
k-gc
k+1
+z
k+1
[0031]
其中,k表示迭代次数,zk、dk表示第k次迭代的z、d,c
k+1
、z
k+1
、d
k+1
表示第k+1次迭代的c、z、d;
[0032]
对a1、a2、a3进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度。
[0033]
在一些优选的实施方式中,a1求解的方法为:
[0034]
将a1进行展开:
其中,i∈[0,1,2];
[0035]
通过lu分解法对展开后的a1加速求解。
[0036]
在一些优选的实施方式中,a2求解的方法为:
[0037]
将a2进行展开:
[0038]
b1:b1:
[0039]
b2:
[0040]
a3:
[0041]
对b1问题,使用moreau近端算子,转化为形式,通过下式求解:
[0042]042][0043]
对b2问题,使用chambolle投影算法,转化为为形式,通过下式求解:
[0044][0045][0046]
其中,p为对偶算子,为梯度算子,div为散度算子,τ表示迭代步长,t表示第t次迭代;
[0047]
对b3问题,使用moreau近端算子,转化为转化为形式,通过下式求解:
[0048][0049]
其中,y=c-d,x=c,θ为极大极小凹约束的内部参数。
[0050]
在一些优选的实施方式中,a3求解的方法为:
[0051]
将a3进行展开进行求解,展开后的a3为:展开后的a3为:
[0052]
本发明的第二方面,提出了一种基于极大极小凹与全变分约束的磁粒子成像图像重建系统,该系统包括:信号获取模块、系统矩阵构建模块、第一方程构建模块、第二方程构建模块、图像重建模块;
[0053]
所述信号获取模块,配置为获取磁粒子成像系统中待成像重建物体的仿体模型在目标区域内逐像素移动测量后采集的时域电压信号,作为输入信号;
[0054]
所述系统矩阵构建模块,配置为将所述输入信号进行傅里叶变换,并基于傅里叶变换后的输入信号构建系统矩阵;
[0055]
所述第一方程构建模块,配置为基于所述系统矩阵,根据磁粒子成像原理,构建磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程,作为第一方程;
[0056]
所述第二方程构建模块,配置为结合各向同性全变分与极大极小凹约束,对所述第一方程进行重构,得到第二方程;
[0057]
所述图像重建模块,配置为通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,进而进行图像重建。
[0058]
本发明的第三方面,提出了一种电子设备,包括:至少一个处理器;以及与至少一个所述处理器通信连接的存储器;其中,所述存储器存储有可被所述处理器执行的指令,所述指令用于被所述处理器执行以实现上述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法。
[0059]
本发明的第四方面,提出了一种计算机可读存储介质,其特征在于,所述计算机可读存储介质存储有计算机指令,所述计算机指令用于被计算机执行以实现上述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法。
[0060]
本发明的有益效果:
[0061]
本发明减小了重建的mpi图像的误差,提升了重建的mpi图像的精度以及质量。
[0062]
本发明通过建立磁粒子成像重建的逆问题求解的数学模型,加入各向同性全变分与无偏估计的稀疏非凸约束项进行重建,采用交替方向乘子法进行迭代优化求解,并在利用lu分解、chambolle投影、硬阈值方法对子问题求解进行加速,得到最优的磁粒子浓度分布,进而进行图像重建。各向同性全变分与极大极小凹的混合正则项相比于当前最常使用的l2正则项约束来说,重建图像的边缘信息更加明显;与storath提出的tv于l1正则项相比,mc非凸正则项能够在保持高度稀疏性的情况下降低l1正则项的偏差效应,减少l1正则项对重建幅值的较低估计,达到误差更小、更精确的重建效果。
附图说明
[0063]
通过阅读参照以下附图所做的对非限制性实施例所做的详细描述,本技术的其他特征、目的和优点将会变得更明显。
[0064]
图1是本发明一种实施例的基于极大极小凹与全变分约束的磁粒子成像图像重建
variation,tv)与无偏估计的稀疏非凸约束项进行重建,采用交替方向乘子法(alternatingdirectionmethodofmultipliers,admm)进行迭代优化求解,并在利用lu分解、chambolle投影、硬阈值方法对子问题求解进行加速,最终得到最优结果。tv与极大极小凹(minimaxconcave,mc)的混合正则项相比于当前最常使用的l2正则项约束来说,重建图像的边缘信息更加明显;与storath提出的结合了tv与l1的融合拉索正则项相比,mc非凸正则项能够在保持高度稀疏性的情况下降低l1正则项的偏差效应,减少l1正则项对重建幅值的较低估计,达到更精确的重建效果。具体如下:
[0082]
s100,获取磁粒子成像系统中待成像重建物体的仿体模型在目标区域内逐像素移动测量后采集的时域电压信号,作为输入信号;
[0083]
在本实施例中,利用充满单位浓度磁纳米粒子的仿体模型(即待成像重建物体的仿体模型)在磁粒子成像系统中的目标区域内逐像素移动进行测量,直至测完目标区域内的每一个像素,保存得到的时域电压信号u
t
。
[0084]
s200,将所述输入信号进行傅里叶变换,并基于傅里叶变换后的输入信号构建系统矩阵;
[0085]
在本实施例中,将时域信号u
t
进行傅里叶变换为u,一次排列组成该种类磁纳米粒子在磁粒子成像系统的中的系统矩阵s∈cm×n,其中,m为频点数,n为像素数。
[0086]
s300,基于所述系统矩阵,根据磁粒子成像原理,构建磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程,作为第一方程;
[0087]
在本实施例中,根据磁粒子成像原理,构建磁粒子成像的浓度c∈rn×1与测量信号u∈cm×1之间的映射方程:
[0088]
sc=u(1)
[0089]
其中,u表示傅里叶变换后的输入信号。
[0090]
s400,结合各向同性全变分与极大极小凹约束,对所述第一方程进行重构,得到第二方程;
[0091]
在本实施例中,通过对映射方程引入各向同性全变分与极大极小凹约束,磁粒子成像重建模型(即第二方程)如下:
[0092]
min
c≥0
λ
tv
||c||
tv
+λ
mc
||c||
mc
(2)
[0093]
s.t.||sc-u||2≤ε(3)
[0094]
其中,||
·
||
tv
为各向同性全变分正则项,||
·
||
mc
为极大极小凹正则项,λ
mc
与λ
tv
为权重系数,||sc-u||2为数据保真项,ε≥0,ε是与噪声大小有关的约束项。
[0095]
s500,通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,进而进行图像重建。
[0096]
在本实施例中,根据构造出的第二方程,利用交替方向乘子法将复杂的混合损失函数优化问题分解为较为简单的子问题进行求解。具体如下:
[0097]
使用增加混合惩罚函数的交替方向乘子法对所述第二方程进行求解,并引入中间变量z=[z
(0)
,z
(1)
,z
(2)
]h、g=[shii]h;
[0098]
将所述第二方程转化为:
[0099]
min
x
f1(c)+f2(z)(4)
[0100]
s.t.pc+qz-s=0(5)
[0101]
其中,p=g,q=-i,s=0,f1(c)=0,f2(z)=ι
e(ε,i,u)
(z
(0)
)+λ
tv
||c||
tv
+λ
mc
||c||
mc
,[
·
]h为共轭转置,ι
e(ε,i,u)
(z0)表示数据保真项的指示函数,i表示单位矩阵;
[0102]
基于转化后的第二方程,结合拉格朗日乘子构造增广的拉格朗日函数:
[0103][0104]
其中,d为拉格朗日乘子,β为惩罚系数,m为频点数;
[0105]
进而得到交替方向乘子法迭代的步骤如下:
[0106]
a1:
[0107]
a2:
[0108]
a3:d
k+1
=d
k-gc
k+1
+z
k+1
(9)
[0109]
其中,k表示迭代次数;
[0110]
a1求解的方法为:
[0111]
将a1进行展开:
[0112][0113]
为了避免对(2i+shs)-1
直接求逆,通过lu分解法对展开后的a1加速求解。
[0114]
a2求解的方法为:
[0115]
将a2进行展开:
[0116]
b1:b1:
[0117]
b2:b2:
[0118]
b3:b3:
[0119]
对b1问题,使用moreau近端算子,转化为形式,通过下式求解:
[0120][0120]
[0121]
对b2问题,使用chambolle投影算法,转化为为形式,通过下式求解:
[0122][0123][0124]
其中,p为对偶算子,为梯度算子,div为散度算子,τ表示迭代步长,本发明优选为0.248,t表示第t次迭代;
[0125]
对b3问题,使用moreau近端算子,转化为转化为形式,通过下式求解:
[0126][0127]
其中,y=c-d,x=c,θ表示极大极小凹约束的内部参数,θ的值影响子问题的凸性,为了使子问题保持凸性,这里令θ=2。
[0128]
a3求解的方法为:
[0129]
将a3进行展开进行求解,展开后的a3为::
[0130]
对a1、a2、a3进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度。
[0131]
具体迭代循环过程如下:
[0132]
s50],获取参数:测得电压信号系统矩阵拉格朗日乘子中间变量正则参数λ
tv
,λmc,ε,β≥0,最大迭代次数k
max
;
[0133]
s502,初始化:d=0,z=0,k=0;
[0134]
s503,判断是否达到设定的终止条件或者k>k
max
;如果是,输出否则,循环迭代:
[0135]
[0136][0137][0138][0139][0140][0141][0142]
k=k+1。
[0143]
另外,为了验证本发明的方法,将本方明方法与使用最广泛的基于l2正则的kaczmarz方法、由storath提出的最先进的融合拉索方法(fussed lasso,fl)在不同信噪比的情况下进行重建实验,对比重建效果。在重建实验中,选取如图3所示的多灰度椭圆仿体,大小为51
×
51。重建的结果如图4所示,并选用结构相似性(ssim)、峰值信噪比(psnr)、归一化均方根误差(nrmse)三项指标进行定量评价。
[0144]
从图4可以看出,本发明提出的方法相比于基于l2正则的kaczmarz方法(图4中简称l2正则)来说,能够保持更多的边缘信息并能更好地抑制噪声;相对于融合拉索方法(图4中简称融合拉索)来说,本发明的方法对仿体的幅值重建更加准确,误差更小。由图5、6、7可以看出,从psnr、ssim、nrmse三个方面对重建结果进行定量评价,本发明的方法都能达到更好的重建效果。
[0145]
本发明第二实施例的一种基于极大极小凹与全变分约束的磁粒子成像图像重建的磁粒子成像系统,如图2所示,该系统包括:信号获取模块100、系统矩阵构建模块200、第一方程构建模块300、第二方程构建模块400、图像重建模块500;
[0146]
所述信号获取模块100,配置为获取磁粒子成像系统中待成像重建物体的仿体模型在目标区域内逐像素移动测量后采集的时域电压信号,作为输入信号;
[0147]
所述系统矩阵构建模块200,配置为将所述输入信号进行傅里叶变换,并基于傅里叶变换后的输入信号构建系统矩阵;
[0148]
所述第一方程构建模块300,配置为基于所述系统矩阵,根据磁粒子成像原理,构建磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程,作为第一方程;
[0149]
所述第二方程构建模块400,配置为结合各向同性全变分与极大极小凹约束,对所述第一方程进行重构,得到第二方程;
[0150]
所述图像重建模块500,配置为通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,进而进行图像重建。
[0151]
所述技术领域的技术人员可以清楚的了解到,为描述的方便和简洁,上述描述的系统的具体的工作过程及有关说明,可以参考前述方法实施例中的对应过程,在此不再赘述。
[0152]
需要说明的是,上述实施例提供的基于极大极小凹与全变分约束的磁粒子成像图像重建系统和/或基于极大极小凹与全变分约束的磁粒子成像图像重建的磁粒子成像系统,仅以上述各功能模块的划分进行举例说明,在实际应用中,可以根据需要而将上述功能分配由不同的功能模块来完成,即将本发明实施例中的模块或者步骤再分解或者组合,例如,上述实施例的模块可以合并为一个模块,也可以进一步拆分成多个子模块,以完成以上描述的全部或者部分功能。对于本发明实施例中涉及的模块、步骤的名称,仅仅是为了区分各个模块或者步骤,不视为对本发明的不当限定。
[0153]
本发明第三实施例的一种电子设备,至少一个处理器;以及与至少一个所述处理器通信连接的存储器;其中,所述存储器存储有可被所述处理器执行的指令,所述指令用于被所述处理器执行以实现权利要求上述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法。
[0154]
本发明第四实施例的一种计算机可读存储介质,所述计算机可读存储介质存储有计算机指令,所述计算机指令用于被计算机执行以实现权利要求上述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法。
[0155]
所述技术领域的技术人员可以清楚的了解到,未描述的方便和简洁,上述描述的电子设备、计算机可读存储介质的具体工作过程及有关说明,可以参考前述方法实例中的对应过程,在此不再赘述。
[0156]
下面参考图8,其示出了适于用来实现本技术系统、方法、设备实施例的服务器的计算机系统的结构示意图。图8示出的服务器仅仅是一个示例,不应对本技术实施例的功能和使用范围带来任何限制。
[0157]
如图8所示,计算机系统包括中央处理单元(cpu,central processing unit)801,其可以根据存储在只读存储器(rom,read only memory)802中的程序或者从存储部分808加载到随机访问存储器(ram,random access memory)803中的程序而执行各种适当的动作和处理。在ram803中,还存储有系统操作所需的各种程序和数据。cpu801、rom802以及ram803通过总线804彼此相连。输入/输出(i/o,input/output)接口805也连接至总线804。
[0158]
以下部件连接至i/o接口805:包括键盘、鼠标等的输入部分806;包括诸如阴极射线管、液晶显示器等以及扬声器等的输出部分807;包括硬盘等的存储部分808;以及包括诸如局域网卡、调制解调器等的网络接口卡的通讯部分809。通讯部分809经由诸如因特网的网络执行通讯处理。驱动器810也根据需要连接至i/o接口805。可拆卸介质811,诸如磁盘、光盘、磁光盘、半导体存储器等等,根据需要安装在驱动器810上,以便于从其上读出的计算机程序根据需要被安装入存储部分808。
[0159]
特别地,根据本公开的实施例,上文参考流程图描述的过程可以被实现为计算机软件程序。例如,本公开的实施例包括一种计算机程序产品,其包括承载在计算机可读介质上的计算机程序,该计算机程序包含用于执行流程图所示的方法的程序代码。在这样的实施例中,该计算机程序可以通过通讯部分809从网络上被下载和安装,和/或从可拆卸介质811被安装。在该计算机程序被cpu801执行时,执行本技术的方法中限定的上述功能。需要说明的是,本技术上述的计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是但不限于:电、磁、光、电磁、红外线、或半导体的系统、装置或器件,或者任意以上的组合。计算机可读存储介质的更
具体的例子可以包括但不限于:具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、ram、rom、可擦式可编程只读存储器(eprom或闪存)、光纤、便携式紧凑磁盘只读存储器(cd-rom)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本技术中,计算机可读存储介质可以是任何包含或存储程序的有形介质,该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本技术中,计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号,其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式,包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质,该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输,包括但不限于:无线、电线、光缆等等,或者上述的任意合适的组合。
[0160]
可以以一种或多种程序设计语言或其组合来编写用于执行本技术的操作的计算机程序代码,上述程序设计语言包括面向对象的程序设计语言,如java、smalltalk、c++,还包括常规的过程式程序设计语言,如c语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中,远程计算机可以通过任意种类的网络,包括局域网或广域网连接到用户计算机,或者可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
[0161]
附图中的流程图和框图,图示了按照本技术各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上,流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分,该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意,在有些作为替换的实现中,方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如,两个接连表示的方框实际上可以基本并行地执行,它们有时也可以按相反的顺序执行,这依所涉及的功能而定。也要注意的是,框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合,可以用执行规定的功能或操作的专用的基于硬件的系统来实现,或者可以用专用硬件与计算机指令的组合来实现。
[0162]
术语“第一”、“第二”等是用于区别类似的对象,而不是用于描述或表示特定的顺序或先后次序。
[0163]
术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含,从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素,而且还包括没有明确列出的其它要素,或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
[0164]
至此,已经结合附图所示的优选实施方式描述了本发明的技术方案,但是,本领域技术人员容易理解的是,本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下,本领域技术人员可以对相关技术特征作出等同的更改或替换,这些更改或替换之后的技术方案都将落入本发明的保护范围之内。
技术特征:
1.一种基于极大极小凹与全变分约束的磁粒子成像图像重建方法,其特征在于,该方法包括:s100,获取磁粒子成像系统中待成像重建物体的仿体模型在目标区域内逐像素移动测量后采集的时域电压信号,作为输入信号;s200,将所述输入信号进行傅里叶变换,并基于傅里叶变换后的输入信号构建系统矩阵;s300,基于所述系统矩阵,根据磁粒子成像原理,构建磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程,作为第一方程;s400,结合各向同性全变分与极大极小凹约束,对所述第一方程进行重构,得到第二方程;s500,通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,进而进行图像重建。2.根据权利要求1所述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法,其特征在于,磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程为:sc=u其中,s表示系统矩阵,c表示磁粒子成像的浓度,u表示傅里叶变换后的输入信号。3.根据权利要求2所述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法,其特征在于,所述第二方程为:s.t.||sc-u||2≤ε其中,||
·
||
tv
为各向同性全变分正则项,||
·
||
mc
为极大极小凹正则项,λ
mc
与λ
tv
为权重系数,||sc-u||2为数据保真项,ε≥0,ε是与噪声大小有关的约束项。4.根据权利要求3所述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法,其特征在于,通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,其方法为:使用增加混合惩罚函数的交替方向乘子法对所述第二方程进行求解,并引入中间辅助变量z=[z
(0)
,z
(1)
,z
(2)
]
h
、g=[s
h i i]
h
;将所述第二方程转化为:s.t.pc+qz-s=0其中,p=g,q=-i,s=0,f1(c)=0,f2(z)=ι
e(ε,i,u)
(z
(0)
)+λ
tv
||c||
tv
+λ
mc
||c||
mc
,[
·
]
h
为共轭转置,ι
e(ε,i,u)
(z0)表示数据保真项的指示函数,i表示单位矩阵;基于转化后的第二方程,结合拉格朗日乘子结合拉格朗日乘子构造增广的拉格朗日函数:其中,d为拉格朗日乘子,β为惩罚系数,m为频点数;
进而得到交替方向乘子法迭代的步骤如下:a1:a2:a3:d
k+1
=d
k-gc
k+1
+z
k+1
其中,k表示迭代次数,z
k
、d
k
表示第k次迭代的z、d,c
k+1
、z
k+1
、d
k+1
表示第k+1次迭代的c、z、d;对a1、a2、a3进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度。5.根据权利要求4所述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法,其特征在于,a1求解的方法为:将a1进行展开:将a1进行展开:其中,i∈[0,1,2];通过lu分解法对展开后的a1加速求解。6.根据权利要求4所述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法,其特征在于,a2求解的方法为:将a2进行展开:b1:b2:b3:对b1问题,使用moreau近端算子,转化为形式,通过下式求解:对b2问题,使用chambolle投影算法,转化为转化为形式,通过下式求解:
其中,p为对偶算子,为梯度算子,div为散度算子,τ表示迭代步长,t表示第t次迭代;对b3问题,使用moreau近端算子,转化为转化为形式,通过下式求解:其中,y=c-d,x=c,θ为极大极小凹约束的内部参数。7.根据权利要求4所述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法,其特征在于,a3求解的方法为:将a3进行展开进行求解,展开后的a3为::8.一种基于极大极小凹与全变分约束的磁粒子成像图像重建系统,其特征在于,该系统包括:信号获取模块、系统矩阵构建模块、第一方程构建模块、第二方程构建模块、图像重建模块;所述信号获取模块,配置为获取磁粒子成像系统中待成像重建物体的仿体模型在目标区域内逐像素移动测量后采集的时域电压信号,作为输入信号;所述系统矩阵构建模块,配置为将所述输入信号进行傅里叶变换,并基于傅里叶变换后的输入信号构建系统矩阵;所述第一方程构建模块,配置为基于所述系统矩阵,根据磁粒子成像原理,构建磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程,作为第一方程;所述第二方程构建模块,配置为结合各向同性全变分与极大极小凹约束,对所述第一方程进行重构,得到第二方程;所述图像重建模块,配置为通过增加混合惩罚函数的交替方向乘子法对所述第二方程进行迭代求解,得到所述待成像重建物体对应的磁粒子成像的浓度,进而进行图像重建。9.一种电子设备,其特征在于,包括:至少一个处理器;以及与至少一个所述处理器通信连接的存储器;其中,所述存储器存储有可被所述处理器执行的指令,所述指令用于被所述处理器执行以实现权利要求1-7任一项所述的基于极大极小凹与全变分约束的磁粒子成像图像重建方法。10.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质存储有计算机指令,所述计算机指令用于被计算机执行以实现权利要求1-7任一项所述的基于极大极小凹
与全变分约束的磁粒子成像图像重建方法。
技术总结
本发明属于磁粒子成像领域,具体涉及一种基于极大极小凹与全变分约束的磁粒子成像图像重建方法、系统、设备,旨在解决现有的磁粒子成像图像重建方法重建的MPI图像误差大、精度差、质量低的问题。本发明方法包括:获取时域电压信号,作为输入信号;将输入信号进行傅里叶变换,并构建系统矩阵;基于系统矩阵,构建磁粒子成像的浓度与傅里叶变换后的输入信号之间的映射方程,作为第一方程;结合各向同性全变分与极大极小凹约束,对第一方程进行重构,得到第二方程;通过增加混合惩罚函数的交替方向乘子法对第二方程进行迭代求解,得到磁粒子成像的浓度,进行图像重建。本发明减小了重建的MPI图像的误差,提升了重建的MPI图像的精度及质量。质量。质量。
技术研发人员:
田捷 朱涛 惠辉 杨鑫 卫泽琛
受保护的技术使用者:
中国科学院自动化研究所
技术研发日:
2022.10.10
技术公布日:
2022/12/30