Efficient Chebyshev Reconstruction for the Anisotropic Equilibrium Model in Magnetic Particle Imaging

TL;DR

This work addresses artifacts in direct Chebyshev MPI reconstruction caused by neglecting nanoparticle anisotropy. It extends the DCR to the EQANIS framework and introduces a fast $p$-rank approximation for the spatially varying deconvolution, achieving $O(N \log N)$ complexity while reducing memory use. Across simulations and six experimental phantoms, DCR-EQANIS consistently improves image fidelity over DCR-EQ and delivers performance comparable to or better than simulated-system-matrix reconstructions, with rank-$p$ results approaching full accuracy as $p$ increases. The approach enables accurate, scalable model-based MPI reconstructions that are practically system-matrix-free, making high-resolution and 3D imaging more feasible in real-world applications.

Abstract

Magnetic Particle Imaging (MPI) is a tomographic imaging modality capable of real-time, high-sensitivity mapping of superparamagnetic iron oxide nanoparticles. Model-based image reconstruction provides an alternative to conventional methods that rely on a measured system matrix, eliminating the need for laborious calibration measurements. Nevertheless, model-based approaches must account for the complexities of the imaging chain to maintain high image quality. A recently proposed direct reconstruction method leverages weighted Chebyshev polynomials in the frequency domain, removing the need for a simulated system matrix. However, the underlying model neglects key physical effects, such as nanoparticle anisotropy, leading to distortions in reconstructed images. To mitigate these artifacts, an adapted direct Chebyshev reconstruction (DCR) method incorporates a spatially variant deconvolution step, significantly improving reconstruction accuracy at the cost of increased computational demands. In this work, we evaluate the adapted DCR on six experimental phantoms, demonstrating enhanced reconstruction quality in real measurements and achieving image fidelity comparable to or exceeding that of simulated system matrix reconstruction. Furthermore, we introduce an efficient approximation for the spatially variable deconvolution, reducing both runtime and memory consumption while maintaining accuracy. This method achieves computational complexity of O(N log N ), making it particularly beneficial for high-resolution and three-dimensional imaging. Our results highlight the potential of the adapted DCR approach for improving model-based MPI reconstruction in practical applications.

Figures (6)

• Figure 1: Example of the anisotropy kernel $\bm K (\bm x, \bm y) = \tfrac{\partial^2}{\partial y_1 \partial y_2}\mathscrbf{E}(\beta \bm{G} \bm y; \mathbb O(\bm x))$ at different positions of the discrete drive-field field-of-view, which is shown in the center of the bottom row. The shape of the convolution kernels changes depending on the position $x$ in the field-of-view. They exhibit a certain symmetry, but the corresponding symmetry axis also changes its orientation. The maximum amplitude also changes depending on the location. The variable $y$ refers to the spatial coordinates within the kernel.
• Figure 2: Behavior of the mean squared error (MSE) of the rank-$p$ approximation of the spatially varying convolution kernel and comparison of the first principal component. (a) rank vs MSE. The error is averaged over all positions within the convolution kernels and all positions of the FOV. (b) Comparison of the Langevin kernel (left) to the first principal component (right) of the rank-$p$ approximation. (c) MSE distribution inside the field-of-view for ranks $p=1,2,...,12$. The error is averaged over all positions within the convolution kernels. The rank is given in the lower left corner of the corresponding MSE heatmap.
• Figure 3: Reconstruction results of the simulation experiments with same computation of $\tilde{\bm c}$ for all methods, followed by different deconvolution methods. 1st row: phantoms used for signal generation. 2nd row: deconvolution results of the DCR-EQ. 3rd - 7th row: deconvolution results for rank-$p$ approximation of DCR-EQANIS with $p=1,3,5,10$ and full rank.
• Figure 4: Reconstruction results of the measurement experiments. 1st row: Photos of the phantoms. 2nd row: Reconstruction using the measured system matrix (SM-MEASURED). 3rd row: Reconstruction using a simulated system matrix following the equilibrium model without anisotropy (SM-EQ). 4th row: Reconstruction using a simulated system matrix following the equilibrium model with anisotropy (SM-EQANIS). 5th row: Reconstruction using the direct Chebyshev reconstruction with the Langevin kernel used for deconvolution (DCR-EQ). 6th-10th row: Reconstruction using the direct Chebyshev reconstruction with the spatially varying anisotropy kernel used for deconvolution (DCR-EQANIS) with rank-$p$ approximation for $p=1,3,5,10$ and full rank.
• Figure 5: Runtime comparison of the DCR-EQANIS with, on the one hand, deconvolution via convolution matrix and, on the other hand, the proposed method via rank-$p$ operator approximation. The average time required for a single FISTA iteration is shown in each case. Left: Behavior of the runtime as a function of the rank for a fixed field-of-view size $N=21^2 = 441$. Right: Runtime as a function of the FOV size.