MPI重建知识整理

​ 最近学习了基于X空间的磁粒子成像理论。正好凑时间复习/整理一下最近学习的磁粒子成像的基本理论框架，包含两部分：

​ 1.基于系统矩阵的MPI重建知识整理

​ 2.基于x空间的MPI重建知识整理

​ 梳理内容包括：仿真常用的公式，成像（仿真）流程。刚刚入门MPI如有错误还望指正！

1. 基于系统矩阵的MPI重建知识整理

1.1 基于测量的系统矩阵方法：

正向信号采集过程：

​ 首先，系统矩阵作为设备的“默认属性“，需要被先获取到，才能进行后续的成像工作；因此，应首先采样FOV内所有位置处的系统函数，然后将所有位置处的系统函数组成系统矩阵（System Matrix, SM）$S$，作为该设备的已知“属性”。获取到系统矩阵后，即可进行后续的扫描成像任务。

​ 对于任意形状的仿体或临床中的任意”目标“，被MPI设备扫描后，通过接收线圈获取到对应的响应信号$\hat{u}$，那么一个完整的成像过程可以被描述为：
$$Sc=\hat{u}(1)$$
​ 其中：
$$\begin{gathered} \hat{u}:=({\tilde{u}}_k^P{)}_{k=0}^{K-1}\in {\mathbb{C}}^K \\ c:=(c(r_n{)}_{n=0}^{N-1})\in {\mathbb{R}}^N \\ S=(\Delta V{\hat{s}}_k(r_n){)}_{k=\mathrm{0...},K-1;n=0,...,N-1}\in {\mathbb{C}}^{K\times N}(2) \end{gathered}$$
​ $k$:表示粒子响应信号的$k$次谐波；$K$：表示考虑的总的谐波数量；$r_n$表示FOV中的第$n$个位置点；$N$表示总的位置数；$c$表示粒子浓度；$\Delta V$：表示位置$n$处样本的体积。

特点：需要采集FOV内所有位置处的系统函数。在FOV中特定的特定点处，放置已知浓度的微小样本进行采样，获取该点处样本的响应信号，作为此处粒子样本的存在标识。然后，移动样本至下一个位置，进行相同的过程，直到采样完FOV中所有的位置点。

重建过程：重建过程，可以被简单的描述为，求解粒子在仿体/目标中的空间分布（浓度）。通过上面的过程描述可知，系统矩阵$S$事先被获取到，信号$\hat{u}$是仿体或任意目标中粒子的响应，这两个变量都可作为已知量，因此可以通过求解上式（1）获得重建图像。常用的重建算法有Kaczmarz、ADMM等。关于Kaczmarz算法的原理可以参考：Kaczmarz迭代算法在磁粒子成像（MPI）中的应用。

优、缺点：暂不具体描述

1.2 基于模型的系统矩阵方法：

​ 在介绍基于模型的系统矩阵之前，首先阐述一下为什么引入基于模型的系统矩阵方法？通过1.2小结的描述可知，基于测量的SM需要逐个扫描FOV中的每个位置以获取对应位置处的系统函数，这一过程非常耗费时间。此外，SM严重依赖设备参数，当MPI设备（磁场、粒子、线圈等）发生变化时，往往需要重新测量SM，进一步增加时间成本。因此，引入基于模型的SM方法，希望在保证重建质量的情况下，降低SM的获取成本。

​ 什么是基于模型的SM方法呢？基于模型中的”模型“是指的描述粒子特性的模型，通常使用郎之万（Langevin）模型和福克普朗克（Fokker-planck）模型。因此，基于模型的SM方法，通常是基于这两种模型（之一）和设备特性推导出理想的SM: ${\hat{s}}_K^{model}$,然后实测少量位置处的系统函数：${\hat{s}}_k^{meas}$，使用实测的${\hat{s}}_k^{meas}$估计真实的SM。

​ 接下来又有一个新的问题：通过实测${\hat{s}}_k^{meas}$和模型${\hat{s}}_K^{model}$估计出真实SM？一般的，可以通过求解优化问题，得到${\hat{s}}_k^{meas}$和${\hat{s}}_K^{model}$之间的放缩系数${\hat{a}}_k$，其中$k$表示频率分量索引, 整个过程用公式表示为：
$$\chi ({\hat{a}}_k)=\sum _{n=0}^{N^{\prime }-1}|{\hat{s}}_k^{meas}(r_n)-{\hat{a}}_k{\tilde{s}}_k^{model}(r_n){|}^2(3)$$

​ ${\hat{a}}_k$可以通过下式求解：
$${\hat{a}}_k=\frac{{\sum }_{n=0}^{N^{\prime }-1}{\hat{s}}_k^{meas}(r_n)({\tilde{s}}_k^{model}(r_n){)}^{\ast }}{{\sum }_{n=0}^{N^{\prime }-1}|{\tilde{s}}_k^{model}(r_n){|}^2}(4)$$
​ 理解：${\hat{a}}_k$其实就是一个缩放因子，用来拟合基于模型计算的系统矩阵和对应位置测量得到的系统矩阵。通过有限个位置的测量系统矩阵，最小化拟合误差求解出一个合适的${\hat{a}}_k$。

2. 基于X空间的MPI重建知识整理

2.1 X空间方法介绍

​ 基于X空间的方法将正向信号产生的过程描述为，粒子分布与点扩散函数（Point Spread Function， PSF）的卷积。此处，将郎之万函数的导数作为系统的PSF。公式表示为：X-Space论文参考
$$s(t)=B_1\frac{d\Phi }{dt}=B_{1}m\rho (x)\ast \dot{L}[kGx]{|}_{x=x_s(t)}kG{\dot{x}}_s(t)(5)$$
​ 怎么由信号得到图像呢？论文中描述粒子密度与PSF的卷积即为重建图像。表示为：
$$\begin{array}{rl}IMG(x_s(t))& =\frac{s(t)}{B_{1}mkG{\dot{x}}_s(t)}\\ & =\rho (x)\ast \dot{L}[kGx]{|}_{x=x_{s}t}\end{array}$$
​

​ 实际计算的时候，通常是将信号与位置对应后，对速度归一化，即为某一位置粒子浓度。基于X空间的方法，理论上较简单，但是目前个人理解的还不深入，仅供参考。

2.2 X空间为什么被叫做X空间

​ 一个问题：X空间方法为什么被叫做X空间？（以下为个人理解）

​ 回答这个问题，需要回到第一篇X空间的文章（上面那篇），文章在Background中提到了一句话我觉得是理解X空间命名的关键，即“In this paper, we will show that MPI can be understood as a x-domain scanning process and, as such, reconstructed quickly and simply”。MPI可以被理解为一个X域的扫描过程。

​ 什么是X域?

​ 这篇文章是第一篇提出X空间的文章，并且只在一维空间验证X空间理论，一维空间用什么表示呢？位置x。文章中的理论证明都是基于位置x推导，如位置x处的总场强，位置x处的粒子磁化强度，位置x处的信号。所以综上，X域就是表示空间内的某一方向。那么X空间（X-Space）就是位置空间（成像空间，或者FOV？），这点在多维X空间的文献中也可以被证实，文章中的理论推导部分使用x表示真实中间中的位置，例$x=[x,y,z{]}^T$。多维X空间文献。

3. 仿真常用公式整理

默认参数：

参数名	符号表示
粒子直径	D
粒子饱和磁化强度	M_S
四氧化三铁密度	d_Fe3O4
四氧化三铁摩尔质量	m_Fe3O4
溶液铁浓度	c_Fe
溶液体积	V_s
真空磁导率	u_0
玻尔兹曼常数	k_B
温度	T_p
接收灵敏度	S

磁场参数：

参数名	符号表示
激励场幅值	A_E
激励场频率	f_E
梯度场梯度	G

公式：

粒子体积：
$$V=\frac{1}{6}\ast \pi \ast D^3$$
磁矩：
$$m=V\ast M_S$$
粒子核质量：
$$m_c=d_{Fe3O4}\ast V$$
粒子核中铁摩尔量：
$$m_{Fe}=\frac{3\ast m_c}{m_{F}e3O4}$$

激励场：
$$H_E=A_E\ast cos(2\ast \pi \ast f_E\ast t)$$
郎之万函数中$\beta$计算：
$$\beta :=\frac{u_0\ast m}{k_B\ast T_p}$$
溶液中四氧化三铁浓度：
$$c=\frac{c_{Fe}}{56\ast m_{Fe}}$$
灵敏度：
$$S=\frac{S}{u_0}$$
FOV范围：
$$FOV=\frac{A_E\ast 2}{G}$$
郎之万函数：
$$L(\xi )=\{\begin{array}{r}(coth(\xi ))-\frac{1}{\xi }\xi \ne 0\\ 0\xi =0\end{array}$$
粒子磁化特性：(此处的H为总磁场强度，可以理解为某位置处粒子随磁场变化的响应)
$$M(H)=cmL(\beta H)$$
粒子磁化特性的导数：
$$M^{\prime }(H)=cm\beta L^{\prime }(\beta H)$$
郎之万函数的导数：
$$L^{\prime }(\xi )=\{\begin{array}{r}(\frac{1}{{\xi }^2}-\frac{1}{sin{h}^2(\xi )})\xi \ne 0\\ \frac{1}{3}\xi =0\end{array}$$
粒子响应信号：
$$u^P(t)=-u_0{\int }_{Object}S(r)\frac{\partial M(r,t)}{\partial t}{d}^{3}r$$