一种加速多步贪婪扩展稀疏Kaczmarz方法

一种加速多步贪婪扩展稀疏kaczmarz方法
技术领域
1.本发明属于磁粒子成像技术领域，具体涉及一种加速多步贪婪扩展稀疏kaczmarz方法、设备。


背景技术：

2.磁粒子成像(mpi，magnetic particle imaging)是一种临床前成像技术，能够可视化磁性纳米粒子的时空分布。这种快速动态过程的图像重建依赖于有效解决不适定逆问题。目前从其测量中重建示踪剂浓度的方法要么适应mpi的图像特征，但计算复杂度较高且收敛较慢，要么速度快但重建图像的图像质量不足。
3.kaczmarz方法基于系统矩阵行进行迭代，适合于在线实时图像重建，目前已经提出了很多改进的kaczmarz方法，如扩展kaczmarz方法和贪婪随机kaczmarz方法。另外，已有工作基于mpi图像的稀疏先验特征和kaczmarz迭代方法，提出稀疏kaczmarz方法，希望在改善mpi图像重建质量的同时取得相对快的重建速度。然而，针对大孔径mpi系统，其系统矩阵规模增大，信号噪声强，经典稀疏kaczmarz的图像重建质量和速度都明显下降。
4.综上，本发明利用mpi图像的稀疏先验特征求得稀疏解；利用扩展kaczmarz的正交投影算法有效消除信号的高斯噪声；利用方程残差信息，每次迭代过程仅选择若干行进行迭代，有效提高收敛速度和重建质量，提出一种用于大孔径磁粒子成像系统实时图像重建的加速多步贪婪扩展稀疏kaczmarz方法。本发明能够实现高质量、快速的大孔径mpi系统图像重建，有望解决人体尺寸mpi实时成像难题。


技术实现要素：

5.为了解决现有技术中的上述问题，即为了解决现有的mpi图像重建方法不能同时实现高质量、快速的大孔径mpi系统图像重建的问题，本发明第一方面，提出了一种加速多步贪婪扩展稀疏kaczmarz方法，应用于大孔径磁粒子成像系统，用于获取磁粒子浓度，进而进行图像重建；该方法包括：
6.s100，获取大孔径磁粒子成像系统的系统矩阵、接收信号、噪声标准差、成像视野内粒子浓度总和以及设定的迭代参数，作为输入数据；所述迭代参数包括迭代过程选用的列数和行数、lasso项系数、步长参数、总浓度约束项系数；
7.s200，基于所述输入数据，构建大孔径磁粒子成像系统图像重建的无约束模型，作为第一模型；对所述第一模型中保真项和总浓度约束项进行合并，进而得到简化后的第一模型，作为第二模型；所述无约束模型包括保真项、总浓度约束项、lasso项、非负实数示性函数项；
8.s300，通过预构建的加速多步贪婪扩展稀疏kaczmarz方法对所述第二模型进行迭代求解，得到所述大孔径磁粒子成像系统中待进行图像重建的目标物体对应的磁粒子浓度，进而进行图像重建。
9.在一些优选的实施方式中，简化后的大孔径磁粒子成像系统图像重建的无约束模
型为：
[0010][0011]
其中，为保真项与总浓度约束项合成后的项，a表示系统矩阵，γ表示总浓度约束项系数，b表示接收信号，c表示成像视野内粒子浓度总和，t表示转置，λ‖φx‖1为lasso范数项，λ表示lasso项系数，φ表示可逆线性变换，x表示磁粒子浓度，i
+
(x)为非负实数集合示性函数项。
[0012]
在一些优选的实施方式中，在通过预构建的加速多步贪婪扩展稀疏kaczmarz方法对所述第二模型进行迭代求解中，终止迭代的条件为：
[0013]
且k≥n
[0014]
其中，k表示迭代次数，δ表示噪声标准差，n表示总迭代次数。
[0015]
在一些优选的实施方式中，在通过预构建的加速多步贪婪扩展稀疏kaczmarz方法对所述第二模型进行迭代求解中，列数的选取方法为：
[0016]
在每次迭代中，计算的绝对值
[0017]
对绝对值进行降序排序，排序后，选取前m1个绝对值对应的列，作为选取的列数；
[0018]
或者对绝对值进行归一化处理，处理后，依据概率不放还随机选取m1列，作为选取的列数；
[0019]
其中，h表示共轭转置，z表示到列空间的正交投影变量。
[0020]
在一些优选的实施方式中，在通过预构建的加速多步贪婪扩展稀疏kaczmarz方法对所述第二模型进行迭代求解中，行数的选取方法为：
[0021]
在每次迭代中，计算的绝对值
[0022]
对绝对值进行降序排序，排序后，选取前m2个绝对值对应的行，作为选取的行数；
[0023]
或者对绝对值进行归一化处理，处理后，依据概率不放还随机选取m2行，作为选取的行数。
[0024]
在一些优选的实施方式中，通过预构建的加速多步贪婪扩展稀疏kaczmarz方法对所述第二模型进行迭代求解，得到所述大孔径磁粒子成像系统中待进行图像重建的目标物体对应的磁粒子浓度，其方法为：
[0025]
s301，初始化参数：k＝1，μ0＝0；vk表示中间变量；
[0026]
s302，当且k＜n，依据绝对值选取列数nc；
[0027]
将nc作为迭代次数，依次迭代求解：
[0028][0029]
判断迭代结束，若结束，则得到判断迭代结束，若结束，则得到并跳转s303，否则迭代求解z
k-1
；其中，m表示大孔径磁粒子成像系统的系统矩阵的行数，t、β、μ、α表示调整yk、xk、vk迭代步长的参数
[0030]
s303，依据绝对值选取行数nr；
[0031]
将nr作为迭代次数，依次迭代求解：
[0032]yk-1
＝αv
k-1
+(1-α)x
k-1
[0033][0034][0035][0036]
判断迭代结束，若结束，则计算稀疏解，并令zk＝z
k-1
，k＝k+1，跳转s302，否则迭代求解x
k-1
。
[0037]
在一些优选的实施方式中，所述稀疏解的获取方法为：
[0038]
xk＝φ-1
shrink(φx
k-1
，λ)
[0039]
x
k-1
＝max(real(x
k-1
)，0)
[0040]
其中，max表示max函数，shrink表示广义收缩算子。
[0041]
本发明的第二方面，提出了一种电子设备，包括：至少一个处理器；以及与至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的加速多步贪婪扩展稀疏kaczmarz方法。
[0042]
本发明的第三方面，提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的加速多步贪婪扩展稀疏kaczmarz方法。
[0043]
本发明的有益效果：
[0044]
本发明能够实现高质量、快速的大孔径mpi系统图像重建。
[0045]
本发明利用mpi图像的稀疏先验特征求得稀疏解；利用扩展kaczmarz的正交投影算法有效消除信号的高斯噪声；利用方程残差信息，每次迭代过程仅选择若干行进行迭代，有效提高收敛速度和重建质量，提出了一种用于大孔径磁粒子成像系统实时图像重建的加速多步贪婪扩展稀疏kaczmarz方法，能够实现高质量、快速的大孔径mpi系统图像重建，有望解决人体尺寸mpi实时成像难题。
附图说明
[0046]
通过阅读参照以下附图所做的对非限制性实施例所做的详细描述，本技术的其他特征、目的和优点将会变得更明显。
[0047]
图1是本发明一种实施例的加速多步贪婪扩展稀疏kaczmarz方法的简略流程示意图；
[0048]
图2是本发明一种实施例的加速多步贪婪扩展稀疏kaczmarz方法的详细流程示意图；
[0049]
图3是本发明一种实施例的加速多步贪婪扩展稀疏kaczmarz方法与其他重建方法重建的结果的对比示意图；
[0050]
图4是本发明一种实施例的适于用来实现本技术实施例的电子设备的计算机系统的结构示意图。
具体实施方式
[0051]
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0052]
下面结合附图和实施例对本技术作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
[0053]
需要说明的是，在不冲突的情况下，本技术中的实施例及实施例中的特征可以相互组合。
[0054]
本发明的加速多步贪婪扩展稀疏kaczmarz方法，如图1所示，包括以下步骤：
[0055]
s100，获取大孔径磁粒子成像系统的系统矩阵、接收信号、噪声标准差、成像视野内粒子浓度总和以及设定的迭代参数，作为输入数据；所述迭代参数包括迭代过程选用的列数和行数、lasso项系数、步长参数、总浓度约束项系数；
[0056]
s200，基于所述输入数据，构建大孔径磁粒子成像系统图像重建的无约束模型，作为第一模型；对所述第一模型中保真项和总浓度约束项进行合并，进而得到简化后的第一模型，作为第二模型；所述无约束模型包括保真项、总浓度约束项、lasso项、非负实数示性函数项；
[0057]
s300，通过预构建的加速多步贪婪扩展稀疏kaczmarz方法对所述第二模型进行迭代求解，得到所述大孔径磁粒子成像系统中待进行图像重建的目标物体对应的磁粒子浓度，进而进行图像重建。
[0058]
为了更清晰地对本发明加速多步贪婪扩展稀疏kaczmarz方法进行说明，下面结合附图1、2对本发明方法一种实施例中各步骤进行展开详述。
[0059]
s100，获取大孔径磁粒子成像系统的系统矩阵、接收信号、噪声标准差、成像视野内粒子浓度总和以及设定的迭代参数，作为输入数据；所述迭代参数包括迭代过程选用的列数和行数、lasso项系数、步长参数、总浓度约束项系数；
[0060]
在本实施例中，先获取输入参数，即获取大孔径磁粒子成像系统的系统矩阵am×n、
接收信号bm×1、噪声标准差δ、成像视野内粒子浓度总和c以及设定的迭代参数，作为输入数据；所述迭代参数包括迭代过程选用的列数nc和行数nr、lasso项系数λ、步长参数η、总浓度约束项系数γ。
[0061]
s200，基于所述输入数据，构建大孔径磁粒子成像系统图像重建的无约束模型，作为第一模型；对所述第一模型中保真项和总浓度约束项进行合并，进而得到简化后的第一模型，作为第二模型；所述无约束模型包括保真项、总浓度约束项、lasso项、非负实数示性函数项；
[0062]
在本实施例中，先构建大孔径磁粒子成像系统图像重建的无约束模型，无约束模型由保真项总浓度约束项(总浓度约束项表示成像视野内的磁粒子浓度总和不变，通过调整系数γ可以增强此约束项作用，γ＝0表示此约束项不起作用)、lasso范数项λ||φx||1(lasso范数项表示x在某种可逆线性变换域内具有稀疏性，当φ为单位矩阵时，表示mpi图像本身具有稀疏性)，非负实数集合示性函数项i
+
(x)构成。
[0063]
然后，通过扩展系统矩阵a和接收信号b，对所述第一模型中保真项和总浓度约束项进行合并，进而得到简化后的第一模型，作为第二模型；即令保真项与总浓度约束项合成从而使得第一模型可以简化为：
[0064][0065]
s300，通过预构建的加速多步贪婪扩展稀疏kaczmarz方法对所述第二模型进行迭代求解，得到所述大孔径磁粒子成像系统中待进行图像重建的目标物体对应的磁粒子浓度，进而进行图像重建。
[0066]
在本实施例中，构建加速多步贪婪扩展稀疏kaczmarz方法。对第二模型进行求解，具体步骤如下：
[0067]
s301，初始化参数：k＝1，μ0＝0，vk表示中间变量；
[0068]
s302，当且k＜n，引入到列空间的正交投影变量，计算绝对值然后依据绝对值选取列数nc；其中，绝对值越大，被选中的概率越大，如此可以加快zk收敛速度。列数选取的具体过程为：在每次迭代中，计算的绝对值对绝对值进行降序排序，排序后，选取前m1个绝对值对应的列，作为选取的列数；或者对绝对值进行归一化处理，处理后，依据概率不放还随机选取m1列，作为选取的列数；其中，h表示共轭转置，z表示到列空间的正交投影变量。
[0069]
将nc作为迭代次数，依次迭代求解：
[0070]
[0071]
判断迭代结束，若结束，则得到判断迭代结束，若结束，则得到并跳转s303，否则迭代求解z
k-1
；
[0072]
s303，用代替经典kaczmarz算法中的计算绝对值减小中噪声部分的范数，提高算法收敛精度，并依据绝对值选取行数nr；其中，绝对值；其中，绝对值越大，被选中的概率越大，如此可以加快xk收敛速度。行数的选取过程具体为：在每次迭代中，计算的绝对值对绝对值进行降序排序，排序后，选取前m2个绝对值对应的行，作为选取的行数；或者对绝对值进行归一化处理，处理后，依据概率不放还随机选取m2行，作为选取的行数。
[0073]
将nr作为迭代次数，依次迭代求解：
[0074]yk-1
＝αv
k-1
+(1-α)x
k-1
ꢀꢀꢀ
(3)
[0075][0076][0077][0078]
在上述步骤中，类比于nesterov加速算法，这里添加了yk、vk，t、β、μ、、α等参数，yk、vk都是加速迭代步骤中的中间变量，t、β、μ、、α为一些调整yk、xk、vk迭代步长的参数，具体可见文献：“xiang x，liu x，tan w，et al.an accelerated randomized extended kaczmarz algorithm[c]//journal of physics：conference series.iop publishing，2017，814(1)：012017”。
[0079]
判断迭代是否结束，若结束，则计算稀疏解，并令zk＝z
k-1
，k＝k+1，跳转s302，否则迭代求解x
k-1
。
[0080]
其中，稀疏解的获取方法为：
[0081]
x
k-1
＝max(real(x
k-1
)，0)
ꢀꢀꢀ
(7)
[0082]
xk＝φ-1
shrink(φx
k-1
，λ)
ꢀꢀꢀ
(8)
[0083]
其中，max表示max函数，shrink表示广义收缩算子。
[0084]
当式成立时，终止迭代过程，得到大孔径磁粒子成像系统中待进行图像重建的目标物体对应的磁粒子浓度，进而进行图像重建。
[0085]
在迭代的过程中，选用列数nc和行数nr时，均不超过总行数的20％。
[0086]
另外，加速多步贪婪扩展稀疏kaczmarz方法与其他重建方法重建的结果的对比示意图，如图3所示，其中，在图3中，各子图从左到右、从上到下依次分别为：图3中的(a)为仿体(即模体)、图3中的(b)为经典的kaczmarz方法对应的重建结果、图3中的(c)为经典的稀疏kaczmarz方法对应的重建结果、图3中的(d)为基于多步最大残差的稀疏kaczmarz方法对应的重建结果、图3中的(e)为基于多步最大残差的扩展稀疏kaczmarz方法对应的重建结
果、图3中的(f)为基于多步最大残差的正则化稀疏kaczmarz方法对应的重建结果、图3中的(g)为基于多步最大残差的加速扩展kaczmarz方法对应的重建结果、图3中的(h)为基于多步最大残差的正则扩展稀疏kaczmarz方法对应的重建结果、图3中的(i)为基于多步最大残差的加速正则扩展稀疏kaczmarz方法对应的重建结果。
[0087]
本发明第二实施例的一种电子设备，包括：至少一个处理器；以及与至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的加速多步贪婪扩展稀疏kaczmarz方法。
[0088]
本发明第三实施例的一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的加速多步贪婪扩展稀疏kaczmarz方法。
[0089]
所述技术领域的技术人员可以清楚的了解到，为描述的方便和简洁，上述描述的电子设备、计算机可读存储介质的具体工作过程及有关说明，可以参考前述方法实例中的对应过程，在此不再赘述。
[0090]
下面参考图4，其示出了适于用来实现本技术方法、系统、装置实施例的服务器的计算机系统的结构示意图。图4示出的服务器仅仅是一个示例，不应对本技术实施例的功能和使用范围带来任何限制。
[0091]
如图4所示，计算机系统包括中央处理单元(cpu，central processing unit)401，其可以根据存储在只读存储器(rom，read only memory)402中的程序或者从存储部分408加载到随机访问存储器(ram，random access memory)403中的程序而执行各种适当的动作和处理。在ram 403中，还存储有系统操作所需的各种程序和数据。cpu 401、rom 402以及ram 403通过总线404彼此相连。输入/输出(i/o，input/output)接口405也连接至总线404。
[0092]
以下部件连接至i/o接口405：包括键盘、鼠标等的输入部分406；包括诸如阴极射线管(crt，cathode ray tube)、液晶显示器(lcd，liquid crystal display)等以及扬声器等的输出部分407；包括硬盘等的存储部分408；以及包括诸如lan(局域网，local area network)卡、调制解调器等的网络接口卡的通讯部分409。通讯部分409经由诸如因特网的网络执行通信处理。驱动器410也根据需要连接至i/o接口405。可拆卸介质411，诸如磁盘、光盘、磁光盘、半导体存储器等等，根据需要安装在驱动器410上，以便于从其上读出的计算机程序根据需要被安装入存储部分408。
[0093]
特别地，根据本公开的实施例，上文参考流程图描述的过程可以被实现为计算机软件程序。例如，本公开的实施例包括一种计算机程序产品，其包括承载在计算机可读介质上的计算机程序，该计算机程序包含用于执行流程图所示的方法的程序代码。在这样的实施例中，该计算机程序可以通过通讯部分409从网络上被下载和安装，和/或从可拆卸介质411被安装。在该计算机程序被中央处理单元(cpu)401执行时，执行本技术的方法中限定的上述功能。需要说明的是，本技术上述的计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是——但不限于——电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子可以包括但不限于：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机访问存储器(ram)、只读存储器(rom)、可擦式可编程只读存储器(eprom或闪存)、光纤、便携式紧凑磁盘只读存储器(cd-rom)、光存储器件、磁存储器
件、或者上述的任意合适的组合。在本技术中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本技术中，计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括但不限于：无线、电线、光缆、rf等等，或者上述的任意合适的组合。
[0094]
可以以一种或多种程序设计语言或其组合来编写用于执行本技术的操作的计算机程序代码，上述程序设计语言包括面向对象的程序设计语言—诸如java、smalltalk、c++，还包括常规的过程式程序设计语言—诸如“c”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(lan)或广域网(wan)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
[0095]
附图中的流程图和框图，图示了按照本技术各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分，该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
[0096]
术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
[0097]
术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
[0098]
至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。