Fourier Magnetic Imaging

• Fourier magnetic imaging is a method that encodes spin and magnetic field data into reciprocal (k) space using pulsed gradients and fast Fourier transforms for high-resolution image reconstruction.
• NV-center implementations employ pulsed optical initialization and dynamical decoupling along with compressed sensing to achieve resolutions from a few nanometers down to sub-nanometer scales.
• Alternative Fourier-domain techniques, including spread-spectrum MRI, holography, and lensless imaging, address inversion challenges with noise regularization and incorporation of physical priors.

Fourier magnetic imaging denotes a class of measurement and reconstruction strategies in which magnetic, spin, or current information is encoded, sampled, or inverted in Fourier space rather than obtained solely by direct real-space rastering. In the nitrogen-vacancy (NV) implementation introduced by Arai et al., pulsed magnetic field gradients phase-encode spatial information on NV electronic spins in wavenumber or k-space followed by a fast Fourier transform to yield real-space images with nanoscale resolution, wide field-of-view (FOV), and compressed sensing speed-up (Arai et al., 2014). Closely related Fourier-domain constructions appear in spread-spectrum MRI, nanoscale nuclear-spin MRI, magnetic current imaging from stray-field maps, and lensless magnetic imaging schemes based on holography, ghost imaging, and Fourier-plane optical processing (Puy et al., 2011, Nichol et al., 2013, Senthilnath et al., 16 Jul 2025, Chen et al., 2018, Martínez et al., 2022, Sandoval et al., 2022).

1. Fourier-domain encoding and image formation

At its most basic, Fourier magnetic imaging follows the standard magnetic resonance relation between an object and its reciprocal-space signal. In Fourier MRI, the measured signal at k-space location $(k_x,k_y)$ is

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$

with $\rho(x,y)$ the transverse magnetization distribution in the slice (Tang, 23 Jun 2025). In nanoscale pulsed MRI, the same structure appears as

$$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$

where the encoding wave-vector is accumulated by the applied time-varying gradient $G(t)$ (Nichol et al., 2013).

For NV-center Fourier magnetic imaging, the phase acquired by a spin at position $\mathbf r$ is

$$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$

and the measured fluorescence at each sampled point in k-space is

$$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$

with $\theta$ carrying any additional phase, for example from a local AC field (Arai et al., 2014). This formal analogy with MRI is exact at the level of reciprocal-space encoding.

The standard reciprocal-space relations between sampling density, FOV, and spatial resolution also carry over. For NV Fourier magnetic imaging,

$$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$

so the k-space step $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$0 sets the FOV and the maximum extent $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$1 sets the nominal spatial resolution (Arai et al., 2014). In Fourier MRI formulated distributionally, rectangular sampling on a grid with interval $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$2 generates periodic replications of the image with period $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$3 in each direction; to avoid aliasing one chooses $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$4 (Tang, 23 Jun 2025). Truncating k-space to $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$5 corresponds to multiplying the signal by a rectangular window, hence to a sinc point-spread function in image space with width $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$6; tapered or Gaussian windowing reduces ringing at the expense of resolution (Tang, 23 Jun 2025).

These relations clarify that Fourier magnetic imaging is not defined by a single sensor platform. What unifies the methods is the use of reciprocal-space sampling, reciprocal-space inference, or Fourier-plane optical processing as the primary route from measurement to magnetic structure.

2. NV-center implementations

The original NV implementation uses optically detected magnetic resonance with pulsed field gradients generated by on-chip microcoils. A 3–5 $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$7s green laser pulse at $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$8 nm initializes the NV into $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$9, a dynamical-decoupling microwave sequence such as spin-echo or $\rho(x,y)$0-pulse CPMG/XY protects coherence during a total free-precession duration $\rho(x,y)$1, pulsed magnetic field gradients are applied during each free-precession interval, and a final $\rho(x,y)$2 pulse maps the accumulated phase into fluorescence read out in the 640–800 nm band (Arai et al., 2014). The gradient coils are anti-Helmholtz pairs separated by $\rho(x,y)$3 $\rho(x,y)$4m; a 1 A pulse generates a $\rho(x,y)$5 G/$\rho(x,y)$6m gradient at the center, uniform within 1% over $\rho(x,y)$7 (Arai et al., 2014).

Uniform k-space grids were demonstrated in one and two dimensions. In a one-dimensional example, $\rho(x,y)$8 was sampled in 58 steps of $\rho(x,y)$9, giving FOV $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$0 $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$1m and $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$2 nm from the standard Fourier-imaging relations (Arai et al., 2014). In two dimensions, grids such as $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$3 or $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$4 points were acquired. After tapered-cosine windowing with 10% taper, a discrete Fourier transform yields a complex real-space image whose magnitude reveals NV positions and whose phase encodes local AC field information (Arai et al., 2014).

Representative performance metrics establish the scale of the method. Arai et al. reported 3.5 nm spatial resolution for 1D single-NV imaging with SNR $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$5, 30 nm pixel resolution for 2D multiple-NV imaging resolving a 121 nm separation, field-of-view up to $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$6 in hybrid real + k imaging with dynamic range $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$7, magnetic sensitivity of $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$8 $$S(k)=\int \rho(r)\,e^{-i\,k\cdot r}\,dr,\qquad k(t)=\gamma\int_0^t G(t')\,dt',$$9T Hz$G(t)$0 per NV, and gradient sensitivity of $G(t)$1 nT/nm Hz$G(t)$2 (Arai et al., 2014). The same work implemented compressed sensing by randomly sampling only $G(t)$3 out of $G(t)$4 k-points, using $G(t)$5 minimization with a partial DFT matrix. For $G(t)$6, the acceleration factor was $G(t)$7, reconstructed NV separations agreed within error bars with full-sample values, and phase-difference field estimates of $G(t)$8 nT were obtained versus $G(t)$9 nT for full sampling (Arai et al., 2014).

Later work pushed the approach to the sub-nanometer regime. A compact ambient platform with thermal drift compensation generated a pulsed magnetic field gradient of up to 13.5 G/$\mathbf r$0m and achieved localization of a single NV center with a spatial resolution of $\mathbf r$1 nm and a magnetic field measurement deviation of 9 nT (Lei et al., 24 Mar 2026). In that implementation, a spin-echo sequence with longest gradient-pulse evolution $\mathbf r$2 $\mathbf r$3s and maximum current $\mathbf r$4 mA produced $\mathbf r$5 G/$\mathbf r$6m at the NV location, corresponding to $\mathbf r$7 nm$\mathbf r$8 and a theoretical $\mathbf r$9 nm, while the observed full-width at half-maximum was $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$0 nm (Lei et al., 24 Mar 2026). The same platform reported $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$1 ms and magnetic sensitivity $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$2 $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$3T/$$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$4 (Lei et al., 24 Mar 2026).

3. MRI, spread spectrum, and optimal Fourier sampling

Fourier magnetic imaging inherits much of its reconstruction theory from MRI. In standard under-sampled MRI, with $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$5 the vectorized image and $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$6 an $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$7-sparse expansion in an orthonormal sparsity basis $$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$8, the measurement model is

$$\phi(\mathbf r)=\gamma\int_0^T \mathbf G(t)\cdot\mathbf r\,dt \approx 2\pi\,\mathbf k\cdot\mathbf r,$$9

where $$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$0 is the row-restricted Fourier operator (Puy et al., 2011). In the spread-spectrum variant, a diagonal unit-modulus modulation $$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$1 or $$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$2 is applied before Fourier sampling:

$$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$3

Digital random pre-modulation uses Rademacher or Steinhaus sequences; the analog version uses a quadratic-phase chirp $$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$4 (Puy et al., 2011).

The central effect of pre-modulation is coherence reduction. Without modulation, the mutual coherence between Fourier sensing rows and sparsity atoms can be high, forcing $$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$5. With random modulation, Puy et al. show that, with probability at least $$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$6,

$$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$7

where the modulus-coherence is

$$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$8

If $$s(\mathbf k)\propto \cos\bigl(2\pi\,\mathbf k\cdot\mathbf r+\theta\bigr),$$9 is the unit-magnitude Fourier basis, then for any $\theta$0 one has $\theta$1, so the modulated coherence is essentially $\theta$2 (Puy et al., 2011).

This leads to the universality claim. Standard compressed sensing gives recovery of every $\theta$3-sparse $\theta$4 by $\theta$5 minimization provided

$\theta$6

After random pre-modulation this becomes

$\theta$7

and for a universal sensing basis such as Fourier or Hadamard the number of measurements scales as

$\theta$8

independent of the choice of sparsity basis $\theta$9 (Puy et al., 2011). Numerical phase-transition experiments at $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$0 showed that, with random pre-modulation, Dirac, Haar, and Fourier sparsity cases collapse onto the same Donoho–Tanner optimal curve (Puy et al., 2011). In numerical MRI examples using $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$1 random phase-encoded lines, chirp modulation with $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$2 improved reconstructions from PSNR $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$3 dB and SSIM $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$4 without chirp to PSNR $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$5 dB and SSIM $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$6 with chirp, while per-iteration complexity remained $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$7 (Puy et al., 2011).

Two later developments address reconstruction beyond standard $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$8 recovery. First, uncertainty quantification for Fourier MRI can be sharpened by reweighting without-replacement sampling so that the virtual Gram matrix satisfies $$\mathrm{FOV}\approx \frac{2\pi}{\Delta k}, \qquad \Delta r\approx \frac{\pi}{k_{\max}},$$9. In a Shepp–Logan phantom experiment with $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$00, $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$01, and complex Gaussian noise with $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$02, reweighted debiasing at $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$03 reduced $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$04 from $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$05 to $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$06, reduced $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$07 from $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$08 to $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$09, and achieved confidence-interval coverage of $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$10 overall and $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$11 on support (Hoppe et al., 2024). Second, random anisotropic sampling combined with dualizable shearlet frames yields asymptotic optimality for cartoon-like functions, with

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$12

and numerical experiments on $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$13 images reported 5–7 dB PSNR gain over wavelet-based schemes (Kutyniok et al., 2015).

4. Fourier inversion of magnetic field maps into current density

A distinct branch of Fourier magnetic imaging reconstructs current density from measured magnetic fields. For a two-dimensional current density $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$14 in the plane $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$15, the out-of-plane field at height $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$16 is given by the Biot–Savart integral

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$17

in the SPIM formulation (Senthilnath et al., 16 Jul 2025). Its two-dimensional Fourier transform satisfies

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$18

together with the divergence-free constraint

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$19

(Senthilnath et al., 16 Jul 2025). The current-to-field operator is therefore diagonal in Fourier space, and the real-space convolution is replaced by multiplication by an inversion kernel.

The SPIM framework adds spatial pre-processing before the FFT-based inversion. Stage I rotates lock-in detector channels $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$20 and $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$21 by a global angle $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$22 chosen to maximize the total variation of the new in-phase image $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$23. In the reported 3D-spiral SQUID scan, $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$24, which sharpens $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$25 by 0.3% and suppresses $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$26 by 25% (Senthilnath et al., 16 Jul 2025). Stage II applies affine alignment with small rotation $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$27 or skew $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$28; values of $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$29 and $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$30 eliminated visible misalignment (Senthilnath et al., 16 Jul 2025). Stage III converts $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$31 into magnetic field using the SQUID calibration constant $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$32 $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$33T/V, computes $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$34 via FFT, forms $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$35 with a hard cut-off $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$36, and uses inverse FFT to recover $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$37 and $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$38. The best results were found for $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$39 with $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$40 $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$41m (Senthilnath et al., 16 Jul 2025).

The classical Fourier route is effective but intrinsically ill-conditioned because the kernel decays exponentially at large $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$42. In the wavelet-based analysis of two-dimensional magnetic current imaging, the kernel is

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$43

so naïve deconvolution amplifies noise strongly at high spatial frequencies (Miller et al., 2024). Fourier methods therefore use low-pass filters such as cosine taper, Gaussian cut-off, or sharp cut-off, but this forces a trade-off between noise and resolution. The reported consequences are explicit: at $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$44 and $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$45, Fourier-filtered error can exceed 30%, whereas an $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$46-curl, divergence-free wavelet method keeps error below 15%; more generally, across all tested noise levels and standoffs, the new method reduces relative $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$47 error by roughly a factor of two compared to the Fourier approach (Miller et al., 2024). The Fourier formulation remains the baseline against which such regularized inversions are measured.

5. Lensless and Fourier-plane magnetic imaging

Related methods recover magnetic structure from Fourier-plane intensity, holographic interference, or higher-order correlations rather than from gradient phase encoding. In thermal-neutron Fourier-transform ghost imaging, a spatially incoherent polarized neutron beam is split into a sample arm with a bucket detector and a reference arm with a position-sensitive detector. The key observable is the fourth-order correlation function

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$48

and, using fermionic Wick factorization and Pauli anti-symmetry, the covariance of intensity fluctuations becomes negative:

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$49

With the choice $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$50, this yields a coincidence signal directly proportional to $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$51, the modulus-squared of the lateral Fourier transform of the projected scattering potential, including both atomic and magnetic terms (Chen et al., 2018). The method is lensless, and the stated resolution limit is the neutron de Broglie wavelength (Chen et al., 2018).

Fourier-transform holography (FTH) provides a different lensless route. The measured hologram is

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$52

so a known reference isolates a phase-locked cross term and overcomes the phase problem in one single step of calculation (Martínez et al., 2022). In magnetic tomography, the phase contrast obeys

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$53

making the reconstructed phase directly proportional to the projection of magnetization along the beam (Martínez et al., 2022). This approach produced a 3D full-vectorial image of a 800 nm-thick extended Fe/Gd multilayer in a 5 $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$54m-diameter circular field of view with a resolution of approximately 80 nm (Martínez et al., 2022).

Fourier-plane optical processing can also be used to disambiguate magnetic spectra. In NV-ensemble vector magnetometry, Fourier-plane amplitude masks and a linear polarizer generate four measurements $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$55, where $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$56 is a $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$57 circulant weight matrix; inversion gives the isolated ODMR lineshapes for each NV orientation even when the eight resonances overlap (Backlund et al., 2017). This enables vector magnetic imaging at arbitrarily low fields and extends the field-dynamic range by a factor $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$58 without increasing $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$59 (Backlund et al., 2017). The measured crosstalk is $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$60, while the relative SNR is approximately 0.47 for NA 1.49/oil and approximately 0.25 for NA 0.75/air (Backlund et al., 2017).

A camera-based Fourier-space MOKE setup uses a high-NA objective to map a broad angular spectrum of incident and reflected wave-vectors onto detector pixels and fits the resulting intensity maps to first-order Kerr expressions. With left and right circular input polarization and no analyzing optics, the method retrieves the three magnetization components together with the optical and magneto-optical constants (Sandoval et al., 2022). In simulations with realistic camera noise, a single shot gives magnetization orientation error $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$61, ten averages give $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$62, and 250 averages give $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$63; errors in optical and Voigt constants fall correspondingly from the 20–30% range to $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$64–2% depending on parameter and averaging (Sandoval et al., 2022).

6. Performance envelope, limitations, and evolving inverse models

Across the surveyed literature, Fourier magnetic imaging spans several experimentally distinct regimes. The following representative values are explicitly reported.

Implementation	Fourier quantity	Representative performance
NV Fourier magnetic imaging	Pulsed-gradient k-space of NV electronic spins	3.5 nm in 1D; 30 nm in 2D; FOV up to $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$65 (Arai et al., 2014)
Sub-nanometer NV localization	Spin-echo k-space encoding with pulsed gradients	$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$66 nm spatial resolution; 9 nT magnetic field deviation (Lei et al., 24 Mar 2026)
Nanoscale Fourier-transform MRI	Pulsed nuclear-spin Fourier encoding	Two-dimensional projection with approximately 10-nm resolution (Nichol et al., 2013)
3D magnetic FTH tomography	Reciprocal-space hologram and inverse FT	800 nm-thick Fe/Gd, 5 $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$67m FOV, approximately 80 nm resolution (Martínez et al., 2022)

Several limitations recur. In gradient-based Fourier imaging, achievable $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$68 is bounded by gradient strength, coherence time, and thermal or timing stability; in the NV case, stronger gradients and longer $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$69 are identified as the route below 1 nm, while in nanoscale nuclear-spin MRI spin relaxation during long encoding reduces high-$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$70 contrast (Arai et al., 2014, Nichol et al., 2013). In Fourier inversion of magnetic field maps, the forward operator is a low-pass filter, so any formal inverse must control noise amplification with cut-offs or regularization, producing the familiar blur-versus-noise trade-off (Miller et al., 2024, Senthilnath et al., 16 Jul 2025). In Fourier-plane optical decomposition, mask throughput and half-pupil blocking reduce sensitivity and introduce anisotropic PSF elongation of about $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$71 along the blocked axis (Backlund et al., 2017). In lensless reciprocal-space methods, phase retrieval is a central issue unless a reference architecture such as FTH is used (Martínez et al., 2022).

Recent work indicates a shift from purely kinematic Fourier inversion toward physics-informed Fourier models. In a Fourier-space approach to magnetization reconstruction from NV stray-field measurements, the forward model uses FFT-based stray-field calculations and Fourier-space upward continuation inside a variational functional

$$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$72

where $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$73 includes exchange, demagnetizing, anisotropy, and Zeeman energies (Setescak et al., 19 Feb 2026). On synthetic data with Néel-type walls at $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$74 nm and 3% Gaussian noise, the recovered height converges within 2 nm of 80 nm over a wide $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$75-plateau; for $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$76 noise the height error remains below 5 nm (Setescak et al., 19 Feb 2026). On $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$77 data with FOV $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$78 nm and $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$79 pixels, the optimal reconstruction gives $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$80–80 nm, reproduces the measured stray-field map to within $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$81 RMS error, typically converges in 500–1 000 iterations, and yields final stray-field RMS residual $$S(k_x,k_y)=\iint_{\mathbb R^2}\rho(x,y)\,e^{-2\pi i\,(x\,k_x+y\,k_y)}\,dx\,dy =\bigl(\mathcal F\{\rho\}\bigr)(k_x,k_y),$$82 MA/m (Setescak et al., 19 Feb 2026). This suggests a convergence between explicit Fourier encoders, FFT-based forward operators, and micromagnetically constrained inverse problems.

The resulting picture is that Fourier magnetic imaging is not a single apparatus but a reciprocal-space paradigm. In one branch, pulsed gradients encode position directly into spin phase and the image is obtained by inverse Fourier transformation. In another, the magnetic observable is reconstructed from Fourier-domain field kernels, holograms, or Fourier-plane optical signatures. In both branches, the decisive quantities are the accessible bandwidth in reciprocal space, the conditioning of the inverse map, and the extent to which physical priors can be incorporated without erasing high-spatial-frequency magnetic information.