【图像重建】基于脉冲激发和附加梯度场的磁粒子成像高分辨率方法研究Matlab复现

✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知，求助可私信。

🔥 内容介绍

磁粒子成像（MPI）作为一种新兴的生物医学成像技术，凭借其高灵敏度和无电离辐射的特点，在分子成像、药物递送监测和细胞追踪等领域展现出巨大的应用潜力。MPI的核心原理在于利用外加磁场激发磁性纳米粒子，通过检测其产生的信号来重建图像。传统的MPI成像方法通常依赖于场强梯度产生的无场线（FFL），利用激发磁场激发FFL上的磁性纳米粒子。通过移动和旋转FFL，并结合投影重建算法，可以获得二维和三维的磁粒子分布图像。

然而，传统MPI技术面临的一个主要挑战是图像分辨率受限于FFL的宽度。FFL的宽度直接受到当前可用的最大梯度磁场强度的限制。较宽的FFL会导致模糊的图像，从而降低空间分辨率。为了克服这一限制，本文提出了一种新的MPI成像方法，该方法引入一个与FFL方向相同的附加梯度磁场，并采用脉冲激发策略进行一维空间编码。

空间编码策略：脉冲激发与附加梯度场的协同作用

本文提出的空间编码策略的核心在于利用脉冲激发磁场和附加梯度磁场的协同作用。具体来说，我们在FFL上施加与FFL方向一致的梯度磁场。脉冲激发磁场采用方波形式，包含一个平坦部分。在平坦部分，磁性纳米粒子的弛豫衰减信号被采集用于一维图像重建。通过改变脉冲激发过程中梯度磁场的变化，我们可以对FFL上的磁性纳米粒子进行空间编码。这相当于对FFL上不同位置的磁性纳米粒子施加不同的激发条件，从而获得包含空间信息的弛豫衰减信号。

系统矩阵的构建与一维图像重建

该方法中的一维图像重建是基于脉冲激发方波平坦部分期间的弛豫衰减信号。通过不同梯度激发剖面获得的弛豫衰减信号构建系统矩阵。该系统矩阵包含了不同梯度激发条件下，磁性纳米粒子弛豫信号的空间分布信息。利用该系统矩阵，我们可以从采集到的弛豫衰减信号中反演出磁性纳米粒子在一维空间上的分布情况。

仿真验证：大尺寸磁性粒子的高分辨率成像

通过仿真研究，我们发现，对于大尺寸的磁性粒子，具有更大平坦部分的脉冲激发磁场能够生成具有更高相关性和更高空间分辨率的一维条形幻象图像。这表明我们的方法对于大尺寸磁性粒子的成像具有明显的优势，并且可以通过调整脉冲激发磁场的参数来进一步优化图像质量。

二维成像：平行FFL移动与高分辨率图像重建

在验证了一维成像方法之后，我们进一步将该方法扩展到二维成像。通过平行移动FFL，并利用所提出的空间编码策略，我们可以在每个FFL位置获取一维投影数据。通过对多个FFL位置采集的数据进行处理，并结合二维重建算法，我们可以重建二维图像。我们利用该方法成功重建了人脑大小的Shepp-Logan幻象和临床横断磁共振血管造影（MRA）数据集的高分辨率二维图像。

总结与展望

本文提出了一种新的基于脉冲激发和附加梯度场的磁粒子成像方法。该方法通过引入与FFL方向一致的附加梯度磁场，并结合脉冲激发策略，实现了一维空间编码。利用弛豫衰减信号和系统矩阵，可以重建高分辨率的一维和二维图像。仿真结果表明，该方法尤其适用于大尺寸磁性粒子的成像，并且可以通过调整脉冲激发磁场的参数来优化图像质量。

这项研究的意义在于，它为MPI技术在空间分辨率方面取得了突破。相较于传统的MPI技术，该方法能够利用现有的梯度场强度获得更窄的有效FFL，从而实现更高的图像分辨率。此外，该方法只需要通过简单的脉冲激发和附加梯度场即可实现空间编码，具有较好的实现性和成本效益。未来，我们希望将该方法扩展到三维成像，并进行体外和体内实验验证，以推动MPI技术在生物医学领域的广泛应用。这项研究有望为开发新一代高分辨率MPI成像系统奠定基础，并为未来疾病诊断和治疗提供新的工具。

📣 部分代码

% 变量初始化[N, M] = size(A);   %[1936,1528]x = complex(zeros(N,1)); residual = complex(zeros(M,1));rowIndexCycle = 1:M;% calculate the energy of each frequency component% 计算每个频率分量的能量energy = sqrt(sum(abs(A.*A),1));  % may use a randomized Kaczmarz% 可能使用随机化的Kaczmarzif shuff    rowIndexCycle = randperm(M);end% estimate regularization parameter% 估计正则化参数lambdZero = sum(energy.^2)/N;lambdIter = lambd*lambdZero;for l = 1:iterations        % 迭代次数    for m = 1:M             % 频率分量        k = rowIndexCycle(m);   % 随机索引                if energy(k) > 0        % 对应频率分量能量大于0            tmp = A(:,k).'*x;            beta = (b(k) - tmp - sqrt(lambdIter)*residual(k)) / (energy(k)^2 + lambdIter);            x = x + beta*conj(A(:,k)); %共轭复数            residual(k) = residual(k) + beta*sqrt(lambdIter);        end    end        if enforceReal && ~isreal(x)        x = complex(real(x),0);    end        if enforcePositive        x(real(x) < 0) = 0;    endendend