High-Resolution Quantitative MRI Using Contrast-Domain Multiscale Energy (MuSE) Regularization

Introduction

Quantitative MRI estimates tissue parameters such as $T_{1}$ that are more readily comparable across scanners, sessions and subjects. Methods such as MPnRAGE fit a signal model to images acquired during magnetization recovery [1]. In neurodegenerative diseases, early changes affect thin structures, such as the cortical ribbon and the hippocampal subfields, whose dimensions approach 1 mm [2]. However, a 0.5 mm isotropic voxel has one eighth of the volume of a 1 mm voxel; because the signal-to-noise ratio scales with voxel volume and the square root of scan time, preserving it would require 64 times the scan time.

Subspace models are the most widely used approach to accelerate quantitative MRI [3]. They constrain the time profile of each voxel to a few principal components of Bloch-simulated signal evolutions, which leaves a few spatial factors to be estimated. These spatial factors can be further regularized by wavelet sparsity [4] or locally low-rank (LLR) penalties [5]. These methods degrade at high acceleration and low SNR (Sec. 2.3). Diffusion models, which learn multiscale denoisers from images, are now emerging as powerful priors for MRI reconstruction [6]. Most perform posterior sampling, which typically needs many network evaluations; extending them to 3-D quantitative imaging with hundreds of contrast volumes is challenging. In addition, fully sampled training data are impractical to acquire at this resolution.

We introduce, to our knowledge, the first diffusion-model-based reconstruction for quantitative MRI at 0.5 mm resolution. It rests on three design choices. First, we use multiscale energy (MuSE) models [7, 8], diffusion models whose denoiser is the gradient of an explicit energy. The explicit energy poses reconstruction as maximum a posteriori (MAP) estimation, which we solve by majorize–minimize (MM) with quadratic, step-size-free updates. This avoids the posterior sampling used by approaches such as DPS and DAPS [9, 10]. Second, to minimize the risk of overfitting limited, basis-specific training data, we apply the prior to the contrast images synthesized from the subspace model [11] rather than to the spatial factors. The prior is then independent of the basis and the sequence, so a generic 2-D model learned from grayscale face images (CelebA) can be used without MRI training data; such generic priors were recently shown to generalize well across imaging domains [12]. Third, we apply this 2-D prior to slices along all three orientations [13], and each stochastic MM iteration regularizes a rotating subset of contrasts along different axes, which makes hundreds of 3-D contrasts tractable.

We evaluate MuSE on 3-D radial MPnRAGE data at 7 T (8.6 min, and a $2\times$ subset with 4.3 min) against CG-SENSE, L1-wavelet and LLR, on a phantom and in vivo. @references大小规律

Background

MPnRAGE acquisition

We use an MPnRAGE [1] sequence, on a 7 T scanner, which repeats an inversion block (Figure 1, step 1). Each block consists of an inversion pulse, 461 3-D radial spokes (TR 4.94 ms, 2.3 s; the first is discarded) and a 499-ms recovery delay, about 2.8 s in total. The flip angle rises from 4° to 8° along a sigmoid. Spokes with the same index in different blocks share the same time since inversion, so each of the $T = 460$ spoke indices defines one contrast $\tau$, sampled by the spokes of all blocks. The 185 blocks take 8.6 min; the last 92 blocks give the 2× scan (4.3 min).

Subspace model

The contrast images $\mathbf{X} = [\mathbf{x}_1,\ldots,\mathbf{x}_T] \in\mathbb{C}^{M \times T}$ follow $\mathbf{x}_\tau(m) = \rho_m s_\tau(\vartheta_m)$, where $\rho_m$ is the proton density and $s_\tau(\vartheta_m)$ the magnetization at time $\tau$. In this work, $\vartheta_m$ is a function of $T_1$, $B_1^+$ and the inversion efficiency (IE) [1]. The measurements are $\mathbf{b}_\tau= \mathcal{A}_\tau\mathbf{x}_\tau+ \mathbf{n}_\tau$, where $\mathcal{A}_\tau$ includes the coil sensitivities, the trajectory and rigid motion, and $\mathbf{n}_\tau$ is Gaussian noise of variance $\nu^2$. Subspace methods write $\mathbf{X} = \mathbf{U}\mathbf{V}$ [3], where $\mathbf{V} \in\mathbb{R}^{R \times T}$ has orthonormal rows estimated from Bloch simulations and $\mathbf{U} \in\mathbb{C}^{M \times R}$ holds the spatial coefficients, so that $\mathbf{x}_\tau= \mathbf{U}\mathbf{v}_\tau$.

Classical regularizers

L1-wavelet regularization penalizes the $\ell_1$ norm of the wavelet detail coefficients of the columns of $\mathbf{U}$ [4]. LLR soft-thresholds the singular values of the $P \times R$ matrix of coefficients in each block of $P$ voxels [5]. If the block has local rank $k \leq R$, the rank-$k$ approximation retains a fraction of about $k/R$ of the noise energy for $P \gg R$, so with $R = 5$ the denoising is limited unless $k$ is well below $R$.

Multiscale energy models

The prior of 2-D images $\mathbf{u}$ at noise level $\sigma$ is $p_w(\mathbf{u};\sigma) \propto\exp(-\phi_w(\mathbf{u};\sigma)/\sigma^2)$, with $\phi_w(\mathbf{u};\sigma) = \frac{1}{2}\lVert\mathbf{u} - \psi_w(\mathbf{u};\sigma)\rVert_2^2$ and a convolutional network $\psi_w$ [7, 8]. Its gradient defines the denoiser

trained by multiscale denoising score matching [14], which makes $D_\sigma$ approximate the MMSE denoiser. The explicit energy allows classical optimization, such as the majorize–minimize scheme below, instead of posterior sampling. Annealing from large to small further improves convergence: early iterations smooth the energy and fix coarse structure, and later ones recover detail [7, 8].

Proposed method

Multiscale MAP formulation using contrast priors

Contrast prior. Unlike SS-MuSE [11], which learns a prior directly on the spatial factors $\mathbf{U}$, we learn it indirectly through the contrast images $\mathbf{X}_\tau= \mathbf{U}\mathbf{v}_\tau$, $\tau= 1,\ldots,T$, which are linear functions of $\mathbf{U}$. A prior on $\mathbf{U}$ requires training data specific to the basis functions $\mathbf{V}$, which unfortunately is not available in our high-resolution setting. A prior on contrast images can instead be learned from generic images, and such priors were recently shown to generalize well to MRI [12].

Slice-wise energy. A 3-D energy would need much 3-D training data and memory, so we apply the 2-D model to slices. For a volume $\mathbf{x} \in\mathbb{C}^{n \times n \times n}$ and an axis $a \in\{x,y,z\}$, the operator $S_a$ extracts the $n \times n$ slices perpendicular to $a$, each an $n \times n$ image, and the real and imaginary parts are treated as separate images [13, 11]. The energy of the volume along $a$ is the sum of the energies of its slices, $e_\sigma(S_a\mathbf{x}) = \sum_j e_\sigma(\mathbf{x}_{a,j})$, where $e_\sigma= -\log p_w(\cdot;\sigma)$ and $\mathbf{x}_{a,j}$ is the $j$th slice.

MAP estimation. We estimate $\mathbf{U}$ by minimizing

where $\eta$ is a regularization parameter and $\mathcal{P}$ is a pool of $N$ frames spread over the inversion recovery. Since $R \ll T$, a pool much smaller than $T$ suffices. By (1), $\nabla e_\sigma(\mathbf{x}) = (\mathbf{x} - D_\sigma(\mathbf{x}))/\sigma^2$, so $E_\sigma$ is never evaluated explicitly.

Majorization

For an ideal MMSE denoiser, the Hessian of $e_\sigma$ is $\sigma^{-2}(\mathbf{I} - \sigma^{-2}\operatorname{Cov}[\mathbf{u} \mid\tilde{\mathbf{u}}]) \preceq\sigma^{-2}\mathbf{I}$ [15], so at $\mathbf{x}_\tau^k = \mathbf{U}_k\mathbf{v}_\tau$ each term is bounded by a quadratic majorizer whose minimization requires no step size (for the trained network, this bound holds approximately),

Rotating contrast subsets

Evaluating all $3N$ terms at every iteration is expensive, so at iteration $k$ we keep the terms of a subset $\mathcal{K}_k \subseteq\mathcal{P}$ of $K \le N$ frames, each assigned one axis $a_{\tau}$, and scale each kept term by $3N/K$ so that the sum approximates the full sum over $3N$ terms. The subset and the axis assignments rotate between iterations, so that MuSE is applied along different axes and, over the iterations, to all frames of the pool. The resulting majorizer is stochastic: it is an unbiased surrogate of (3) summed over all terms if the frames and axes are drawn uniformly [16], and rotating through interleaved subsets approximates this. Each iteration costs $K$ 2-D passes over the volume, and the descent is not monotone. The selected contrasts are denoised, $\hat{\mathbf{x}}_{\tau}=D_{\sigma}(\mathbf{x}_{\tau}^{k})$ along $a_{\tau}$, and mapped back to the subspace by least squares,

where $\mathbf{V}_{\mathcal{K}_k}$ holds the columns of $\mathbf{V}$ of the selected frames and $\epsilon>0$ is a small constant that stabilizes the inverse. Approximating $\left\|(\mathbf{U}-\widetilde{\mathbf{U}}_k)\mathbf{V}_{\mathcal{K}_k}\right\|_F$ by $\left\|\mathbf{U}-\widetilde{\mathbf{U}}_k\right\|_F$ and absorbing constant scales into $\eta$ gives

with $\lambda_k=\eta^2/\sigma_k^2$, solved by a few warm-started conjugate gradient (CG) iterations.

Continuation and renoising

As $\sigma_k$ is annealed, $\lambda_k$ increases, so the weight on mildly denoised anchors grows. The denoiser is trained on inputs of the form $\mathbf{u}+\sigma\mathbf{z}$, whereas the iterate $\mathbf{x}_{\tau}^{k}$ carries residual aliasing and noise that are not Gaussian at level $\sigma_k$. To reduce this train–test mismatch, as in D-MOD [17], the anchor can be averaged over $L$ renoised inputs, $\hat{\mathbf{x}}_{\tau}=\frac{1}{L}\sum_l D_{\sigma}(\mathbf{x}_{\tau}^{k}+\sigma\mathbf{z}_l)$ with $\mathbf{z}_l\sim\mathcal{N}(0,\mathbf{I})$; $L=0$ means no renoising.

Experiments

Validation strategy and BrainWeb phantom

No ground truth is available for 0.5 mm quantitative MRI in vivo, so we validated the methods on a numerical phantom with known $T_1$ maps and contrasts. The BrainWeb anatomical model [18] provides 1 mm segmentation maps, which we interpolated to 0.5 mm. We assigned representative $T_1$, $B_1^{+}$ inversion efficiency (IE) and proton density (PD) values to gray matter, white matter, CSF and skull, which gives the parameter maps $\vartheta_m$ and $\rho_m$ of the signal model in Sec. 2. The ground-truth contrasts ($T=460$) were generated by Bloch simulation of the MPnRAGE sequence with the blocks, spoke timing and flip-angle schedule of the in vivo scan. The k-space data were then simulated with the forward model $\mathcal{A}_{\tau}$, that is, the in vivo radial trajectory and coil sensitivities, without motion, with in vivo-matched noise added.

In vivo data

One healthy volunteer was scanned on a 7 T GE scanner with a 32-channel head coil (Sec. 2). Images were reconstructed at 0.5 mm isotropic resolution ($512^3$, 256 mm FOV) from 16 compressed virtual coils, with per-block rigid motion included in $\mathcal{A}_{\tau}$.

Training and use of the prior

The network $\psi_w$ is a 3.4 M-parameter SongUNet with EDM preconditioning [19]. It was trained by denoising score matching on 100k grayscale CelebA face images, with no MRI data. To apply it to MRI, each synthesized contrast volume is cut into 2-D slices along the chosen axis. Because the network is trained on real-valued images, the real and imaginary parts of each slice are denoised separately as two grayscale images and then recombined into the complex slice, $\hat{\mathbf{x}}=D_{\sigma}(\operatorname{Re}\mathbf{x})+iD_{\sigma}(\operatorname{Im}\mathbf{x})$.

Compared methods, tuning and metrics

All methods share the basis ($R=5$), the forward model (scaled to $\lambda_{\max}(\mathcal{A}^{H}\mathcal{A})=1$) and the coil maps. CG-SENSE [20] uses 50 CG iterations and initializes the regularized methods. MuSE uses a pool of $N=8$ frames (frames 10–420, TI 69–2095 ms) with MuSE applied along rotating axes (here $K=N$, so the subset is the full pool), no renoising, and the first seven levels of a ten-level EDM schedule from $\sigma_0$ (the noise level of the CG-SENSE images) to $0.08\sigma_0$. Each MuSE iteration uses 20 CG iterations. L1-wavelet uses a six-level orthonormal Daubechies-4 transform and ADMM (10 iterations of 5 CG steps); LLR uses shifted $8^3$ blocks, the threshold $c(\sqrt{P}+\sqrt{R})$ in noise units and ADMM (10 iterations). At $1\times$ and $2\times$, the MuSE weight is $\eta/\sigma_0=0.2$ and 0.15, the wavelet weight is $\lambda=3.2\times10^{-3}$ and $1.6\times10^{-3}$, and the LLR constant is $c=0.25$ and 0.15.

The parameters of all methods were tuned on the phantom to minimize the NRMSE in the brain, averaged over the WM-nulled, MPRAGE-like and CSF-nulled frames (TI 415/1398/855 ms). In vivo, the noise-relative $c$ and $\eta$ were kept, and the wavelet $\lambda$ was scaled by the ratio of noise levels and data gains (0.69 at $1\times$, 0.82 at $2\times$). Figures show contrasts synthesized with $B_1^{+}=\mathrm{IE}=1$. Fine detail was assessed by the high-frequency error norm (HFEN; 3-D LoG, $\sigma=1.5$ voxels) in the brain, and $T_1$ by the median bias and interquartile range (IQR) in WM and GM on three orthogonal slices.

Results and Discussion

Phantom comparison. MuSE gave the lowest whole-brain NRMSE at both scan times (Table 1); at 2×, it was 44% lower than CG-SENSE and 19% lower than L1-wavelet, and even below the 1× L1-wavelet result (0.074 versus 0.080). It also had the lowest HFEN and, at 2×, the narrowest $T_1$ IQR in WM (38 ms, versus 47–117) and GM (159 ms, versus 193–325). L1-wavelet was the closest competitor, and the two methods are nearly equal on the high-SNR MPRAGE-like frame (NRMSE 0.0426 versus 0.0431; PSNR 35.9 versus 35.5 dB). The difference lies in the low-SNR images (Figure 2). In the WM-nulled frame, L1-wavelet leaves a mottled, patchy texture, whereas MuSE is smoother and keeps the thin cerebellar white matter branches (PSNR 34.0 versus 32.2 dB; SSIM 0.918 versus 0.853). The $T_1$ map shows the same pattern: L1-wavelet is speckled, whereas MuSE is homogeneous

	CG-SENSE	LLR	L1-wavelet	MuSE
NRMSE, 1×	0.119	0.105	0.080	0.058
HFEN, 1×	0.197	0.217	0.253	0.146
NRMSE, 2×	0.133	0.110	0.091	0.074
HFEN, 2×	0.317	0.249	0.274	0.244
$T_1$ bias/IQR, WM	-42/117	-6/66	-2/47	-19/38
$T_1$ bias/IQR, GM	-34/325	-34/251	-30/193	-18/159

with sharper folia boundaries (NRMSE 0.159 versus 0.191). MuSE had a WM $T_1$ bias of −19 ms at 2×.

In vivo. Figure 3 shows the in vivo reconstructions of the four methods at 1× and 2×. Since no ground truth exists, we compare them only visually, and a quantitative in vivo evaluation is left for future work.

Limitations

This study is preliminary: it uses one subject and one phantom, the in vivo comparison is visual, memory use is high (about 115 GB on each of two GPUs at 1×), and the face-trained prior is untested in pathology. More elaborate studies will follow.

Conclusion

MuSE regularizes subspace quantitative MRI with a multi-scale energy prior on synthesized contrasts, learned from face images and minimized by a step-size-free stochastic MM algorithm. At 0.5 mm, it outperformed the other methods on a phantom and was demonstrated in vivo at 7 T.

Compliance with Ethical Standards

This study was performed in line with the principles of the Declaration of Helsinki. In vivo data were acquired at the University of Iowa with approval of the Institutional Review Board and written informed consent from the subject.

Acknowledgments

Supported by NIH grants R01-AG067078, R01-EB031169 and R01-EB019961. We thank UVA Research Computing (rc.virginia.edu) for computational resources. Claude (Anthropic) was used to refine the text of this paper. The authors reviewed the content and take full responsibility for it.