k-space weighted image average (KWIA) for ASL-based dynamic MR angiography and perfusion imaging.

A novel denoising algorithm termed k-space weighted image average (KWIA) was proposed to improve the signal-to-noise ratio (SNR) of dynamic MRI, such as arterial spin labeling (ASL)-based dynamic magnetic resonance angiography (dMRA) and perfusion imaging. KWIA divides the k-space of each time frame into multiple rings, the central ring of the k-space remains intact to preserve the image contrast and temporal resolution, while outer rings are progressively averaged with neighboring time frames to increase SNR. Simulations and in-vivo dMRA and multi-delay ASL studies were performed to evaluate the performance of KWIA under various MRI acquisition conditions. SNR ratios and temporal signal errors between KWIA-processed and the original data were measured. Visualization of dynamic blood flow signals as well as quantitative parametric maps were evaluated for KWIA-processed images as compared to the original images. KWIA achieved a SNR ratio of 1.73 for dMRA and 2.0 for multi-delay ASL respectively, which were in accordance with the theoretical predictions. Improved visualization of dynamic blood flow signals was demonstrated using KWIA in distal small vessels in dMRA and small brain structures in multi-delay ASL. Approximately 5% temporal errors were observed in both KWIA-processed dMRA and ASL signals. Fine anatomical features were revealed in the quantitative parametric maps of dMRA, and the residuals of model fitting were reduced for multi-delay ASL. Compared to other conventional denoising methods, KWIA is a flexible denoising algorithm that improves the SNR of ASL-based dMRA and perfusion MRI by up to 2-fold without compromising spatial and temporal resolution or quantification accuracy.

Introduction

Arterial spin labeling (ASL) uses magnetically labeled arterial blood as an endogenous tracer to noninvasively measure cerebral blood flow (CBF) [1]. In a typical ASL scan, image acquisition is performed after a single post-labeling delay (PLD) that allows the labeled blood to flow into the brain region of interest for quantitative estimation of CBF [2]. The caveat of such single-delay ASL scan is that the assumption of PLD ≥ arterial transit time (ATT) may be violated, and thus affecting the accuracy of CBF quantification. Multi-delay ASL would be ideal that enables the visualization of flow dynamics and quantification of multiple hemodynamic parameters [3], such as CBF and ATT for tissue compartment, and arterial cerebral blood volume (aCBV) and arterial bolus arrival time (aBAT) for macrovascular compartment. Another emerging application of ASL is non-contrast enhanced time-resolved 4D dynamic magnetic resonance angiography (dMRA) [4-6] that provides the depiction of both the architecture of cerebral vasculature with high spatial resolution (mm) and dynamic flow pattern of labeled blood with high temporal resolution (a few hundred ms). Quantitative hemodynamic parameters are also derived from dMRA, such as arterial blood flow (aBF), arterial blood volume (aBV), ATT, and time-to-peak (TTP) [7].

However, the intrinsically low signal-to-noise ratio (SNR) is a critical limitation hampering the clinical translation of ASL-based techniques [8]. Typically, several repetitions of ASL acquisitions and/or prolonged labeling bolus are required to achieve sufficient SNR resulting in increased scan time [2,6]. This will lead to impractically long scan times for multi-delay ASL and/or 4D dMRA scans, and often compromises have to be made between spatial and temporal resolutions. View-sharing techniques [9-11], such as k-space weighted image contrast (KWIC) [9], as well as sparse sampling with constrained reconstruction [12] have been proposed to improve temporal and/or spatial resolution of 4D dMRA and multi-delay ASL without prolonged scan time. These prospective acceleration methods, nevertheless, all require special pulse sequences with a predefined order for k-space data acquisition.

The purpose of this study was to present a novel post-processing technique termed k-space weighted image average (KWIA) to improve the SNR of ASL-based dMRA and multi-delay perfusion imaging, without compromising spatiotemporal resolution or quantification accuracy. KWIA is a denoising technique originally proposed for low dose computed tomography perfusion (CTP) [13], which has the advantages of computational simplicity without altering the image texture or resolution as compared to conventional denoising techniques such as iterative reconstruction. While KWIA was originated from the concept of k-space, it was never applied to MRI. In contrast to conventional view-sharing and sparse sampling techniques that prospectively undersample k-space data, KWIA aims to retrospectively improve the SNR of dynamic MRI acquired with various acquisition conditions and acceleration methods. In this paper, we first presented the theoretical framework of KWIA and then evaluated its performance in terms of SNR, depiction of dynamic flow patterns, temporal blurring, and quantification accuracy using a FORBILD [14] digital phantom and in-vivo data. KWIA was also compared to other conventional denoising methods to demonstrate its advantage in preserving temporal and spatial resolution.

Theory

Without loss of generality, our analysis focuses on MR signal acquired in k-space with Cartesian trajectories. Each sample in k-space contains not only MR signals but also the noise arising from electronics and Brownian motion of molecules in the human body. The noise embedded in the signal can be considered as uniformly distributed additive Gaussian noise with zero-mean for both real and imaginary components. A 2D complex MR image obtained from this k-space using inverse fast Fourier transform (IFFT) should have the same Gaussian noise property [15], and its noise standard deviation (SD) σ is determined by the σ from all k-space samples as shown in Eq. (1): where N and N represent the image size along two dimensions which also serve as the scaling factors of IFFT, σ(k, k) is the SD of complex Gaussian noise in both real and imaginary parts at the k-space coordinate (k, k). According to Eq. (1), manipulating σ across different k-space locations could lead to a global reduction of σ. For a timeseries of images, this can be achieved by frequency and time-weighted averaging of k-space data of neighboring time frames. To better preserve temporal fidelity, more averaging is applied on outer k-space (high frequency component) as the image contrast is primarily determined by the central k-space (low frequency component). Meanwhile, high spatial resolution is maintained since the extent of k-space coverage remains unchanged and signal content of high spatial frequencies preserved during averaging. This is the basic idea of KWIA.

Specifically, KWIA first converts MR images to k-space data by FFT, then the k-space of each image in a time series is divided into multiple rings (other shapes can be used but are less efficient for improving SNR as discussed below). A symmetric moving averaging window with a window size determined by N, the number of rings, slides from the beginning to the end of the time series. When the KWIA window is applied to a selected time frame t′ (i.e., the window center is at t′), the central circle of this k-space remains intact, while outer rings are progressively averaged with neighboring time frames to increase SNR (see an example of 3-ring KWIA in Fig. 1(a)). The radius of the central circle r1 is a user-defined parameter, and the remaining radius of the k-space is equally divided into N − 1 rings: where r is the radius of the i circle, and R is the radius of the full k-space coverage. The time and frequency-dependent weighting function W(r, t, t′) is determined by Eq. 3 for the k-space data at the time frame t, with the KWIA window center at the time frame t′, and radius r. Note that t′ is the discrete index ranging from 1 to the total number of time frames, and t is the discrete index ranging from t′ − N + 1 to t′ + N − 1.

Fig. 1.

(a) Schematic diagram of 3-ring KWIA. When KWIA is applied to the selected time frame, 5 frames within the KWIA window are weighted averaged to form a KWIA-processed k-space. The weighting function and KWIA-processed k-space are shown. (b) Predicted theoretical SNR ratio of KWIA with 2, 3, and 4 rings. SNR ratio increases with a greater number of rings (N) and a smaller ratio (q) of the first ring to the full radius.

For Cartesian k-space sampling, the weighting function can be discretized into the samples on k and k, so that the KWIA-processed k-space data can be expressed by: where F and F represent the k-space data after and before KWIA processing, respectively. Note that for all k-space locations of the selected time frame t′. Since t ranges from 1 to the total number of frames, the KWIA window could be incomplete and asymmetric for the first and last few frames. The equation for expected SNR ratio (ratio of SNR after KWIA versus original SNR) of time frame t′ that has a full KWIA window can be given as below.

Fig. 1(b) shows the predicted SNR ratio as a function of q (r1/R) for 2, 3 and 4 rings respectively. The SNR ratio increases with a higher N, and a smaller q, the ratio of the first ring to the full radius (r1/R). For specific applications, these two parameters can be determined empirically according to the optimal trade-off between SNR improvement and temporal fidelity for that application. The KWIA-processed k-space data will then be converted to denoised images using IFFT. Note that magnitude MRI can be commonly used for KWIA denoising which follows Rician noise distribution that deviates from Gaussian noise for the background and low SNR regions. This may cause deviation from the predicted SNR ratios for practical implementation of KWIA as discussed below.

Methods

Digital phantom simulation

As shown in Fig. 2(a), a 128-by-128 FORBILD digital phantom was simulated with 15 time frames of dynamic signals within three vessels (blue arrow) (16, 4, and 1 pixel) representing large, median, and small arteries in dMRA, respectively. Two tissue region-of-interests (ROIs) (160 pixels) representing gray (yellow arrow) and white matter (red arrow) were simulated for multi-delay ASL signals to evaluate the performance of KWIA. For dMRA, the dynamic signal of the large artery was generated using gamma variate function [16] with shape parameter = 3 and rate parameter = 107 ms which was used as the arterial input function (AIF) based on the average AIF measured in dMRA across five subjects. Subsequently, the dynamic signals of medium and small arteries were simulated by convolving the AIF with an exponential decay residual function [17]. The three dynamic signals are shown in Fig. 2(c). The arterial blood flow (aBF) and arterial transit time (ATT) were set to 360 ml/ml/min (interpreted as ml blood flowing through ml volume per minute) and 300 ms respectively for the medium vessel, and 100 ml/ml/min and 500 ms for the small vessel, respectively, based on the average of measured values in five subjects in the in-vivo dMRA experiment. In addition, gray matter (GM) and white matter (WM) signals (Fig. 2(d)) were simulated at the post labeling delays (PLD) from 500 ms to 2740 ms with an interval of 160 ms based on the single compartment kinetic model of pCASL [3,18] (see online supplementary materials Section 1), where CBF and ATT were set to 60 ml/100 g/min and 1 s for GM, and 20 ml/100 g/min and 1.5 s for WM, respectively.

Fig. 2.

Simulated phantom images (8th time frame) and signal curves. (a) An image without noise. Blue, yellow, and red arrow indicate 3 dMRA vessels (16, 4, 1 pixels), GM and WM, respectively. (b) An image with simulated noise added in k-space. The noise is Gaussian, complex, additive, uniform, zero-mean. (c) dMRA signals in large (blue), medium (red), and small vessel (yellow). (d) ASL signals in GM (blue) and WM (red). (e) Images of 3 lower noise SD conditions (1.41-fold SNR, 1.73-fold SNR, and 2-fold SNR) as well as 6 KWIA conditions applied on the original noisy data (Fig. 2 (b)) that improved SNR by 1.41 (N=2 and q=36/64, and N=3 and q=42/64), 1.73 (N=3 and q=25/64, and N=4 and q=30/64), and 2 times (N=3 and q=11/64, and N=4 and q=20/64). All images processed by KWIA achieve similar SNR level as the corresponding lower noise SD images.

Simulated complex zero-mean Gaussian noise with σ = 32 in both real and imaginary parts (σ = 0.25 in magnitude image according to Eq. (1)) was added in the k-space of the phantom to achieve SNR levels, which are defined as the mean signal divided by noise SD, of 20, 20, 10, 4, and 1.5 for large vessel, medium vessel, small vessel, GM, and WM signal, respectively. By applying IFFT and taking the magnitude of complex numbers, the simulated noisy magnitude image of the phantom can be obtained as shown in Fig. 2(b). Different MRI acquisition conditions could alter the noise properties in k-space, and influence the performance of KWIA. To demonstrate the compatibility of KWIA with some common MRI acquisition conditions, the performance of KWIA was evaluated under the standard condition (square FOV without any filtering or acceleration), rectangular FOV, low-pass filtering, zero filling, partial FT, parallel imaging with SENSE and GRAPPA. Testing of all conditions were performed using the magnitude images as inputs. For rectangular FOV, the width of images was cut to 96 (from 128), and the rings of KWIA were reshaped into elliptical rings according to the ratio of length to width (R was defined as the length of the shortest side). Low-pass filtering (Fermi filtering [19] with Fermi radius = 0.8 and Fermi width = 0.046) was applied on the k-space data along the phase-encoding (PE) direction (top to bottom). Zero filling was implemented by replacing peripheral 50% k-space entries with 0 along PE direction. Partial FT was achieved by filling the last 25% of k-space entries along the PE direction with the conjugate symmetry of the first 25% PE k-space data. Finally, in order to perform parallel imaging, the original phantom image was multiplied with 32 simulated coil sensitivity maps before adding noise. Both SENSE and GRAPPA were applied along PE direction with an accelerating factor of two.

KWIA was implemented in MATLAB (The Mathworks Inc., Natick, MA, USA) on the magnitude images of the digital phantom with the configurations of N = 2 and q = 36/64, and N = 3 and q = 42/64 to achieve a SNR ratio of 1.41, N = 3 and q = 25/64, and N = 4 and q = 30/64 to achieve a SNR ratio of 1.73, and N = 3 and q = 11/64, and N = 4 and q = 20/64 to achieve a SNR ratio of 2, respectively. Detailed configurations are listed in Table 1.

Table 1

KWIA Configurations to achieve 1.41, 1.73 and 2-fold SNR ratio for digital phantom and in-vivo data respectively.

	Digital phantom (R = 64)						In-vivo dMRA (R = 88)	In-vivo ASL (R = 48)
SNR ratio	1.41		1.73		2		1.73	2
	N = 2	N = 3	N = 3	N = 4	N = 3	N = 4	N = 3	N = 4
R1	36	42	25	30	11	20	33	15
R2	64	53	45	41	38	35	61	26
R3		64	64	53	64	49	88	37
R4				64		64		48

To evaluate the performance of KWIA, SNR was measured on the original noisy images and images processed by KWIA under each of the MRI conditions respectively. SNR is determined by the mean signal of the GM ROI divided by the SD (corrected for Rayleigh distribution) of the pixels in the background ROI averaged over time frames that have a full KWIA window (first and last few time frames were excluded). The SNR ratios of KWIA were calculated by dividing the SNR of KWIA-processed images by the SNR of the original images.

The potential temporal error caused by KWIA was measured using normalized root mean square error (NRMSE) between the original time courses without simulated noise and KWIA processed time courses in five ROIs (three vessels and two tissue ROIs) respectively. The NRMSE was defined as the square root of the average of squared errors normalized by the range of the original time courses.

In-vivo data acquisition

Five healthy subjects were recruited for ASL-based dMRA (2 male, age = 24.9±2.0 years) and another five healthy subjects were recruited for multi-delay ASL experiments (4 male, age = 28.2±6.4 years). All subjects signed written informed consent according to protocols approved by the institutional review board of University of Southern California. Intracranial dMRA images were acquired on a Siemens Tim Trio 3 T scanner with a 12-channel head coil using a non-contrast enhanced dMRA protocol with the flow-sensitive alternating inversion recovery (FAIR) labeling scheme and GRE-based Cartesian acquisition (voxel size = 1 × 1 × 1.5 mm, matrix = 220 × 176 × 32, flip angle = 12°, bandwidth = 720 Hz/pixel, TR = 5.6 ms, TE = 3.19 ms, 75% partial FT on both phase and slice directions; in-plane GRAPPA with R = 2, total 6 min). The imaging slab covers the circle of Willis and its branches. A total of 12 time frames were acquired from inversion time of 255 to 1414 ms with an interval of 105 ms. Multi-delay ASL experiment was performed on a Siemens Prisma 3 T scanner with 32-channel head coil using a 14-delay pseudo-continuous ASL (pCASL) with 3D segmented GRASE readout. Two nonselective Hyperbolic Secant (HS) inversion pulses were applied during the PLD and the timing of background suppression was optimized according to [20] for each PLD. Imaging parameters were as follows: voxel size = 2 × 2 × 2 mm3, matrix = 96 × 96, 24 slices covering the basal ganglia and lateral ventricles, 81% partial FT on PE direction, TE = 29 ms, TR = 4060 ms, label duration = 1500 ms, 14 PLDs from 500 ms to 2580 ms with an interval of 160 ms, 8 min 40 s for each scan with an M0 image. A 3D MPRAGE scan (voxel size = 1 × 1 × 1 mm3, matrix = 256 × 240 × 192, TR = 2300 ms, TE = 2.98, scan time = 5 min 12 s) was acquired for each subject.

KWIA processing of dMRA

Control and label images were subtracted to obtain angiography images. KWIA was applied to the magnitude angiography images with the configuration of N = 3 and q = 33/88 to achieve a SNR ratio of 1.73. This KWIA configuration was selected to balance the trade-off between SNR improvement and potential temporal error according to the simulation result in Section 4.1 and the in-vivo result in Section 2 of online supplementary materials. The KWIA-processed dMRA images were compared to the original dMRA images to evaluate the dynamic flow pattern, SNR, temporal error, and quantification accuracy. The middle cerebral artery (MCA), posterior cerebral artery (PCA) and distal small vessel ROIs were drawn manually on all five subjects. The dMRA signal curve extracted from MCA was used as AIF. SNR was measured by dividing the mean of AIF by the SD (corrected for Raleigh distribution) measured in the background region and averaged from 3 to 10 time frame that had a full KWIA window. The temporal error was evaluated using signal curves from the MCA, PCA, and distal small vessel ROIs, and NRMSE was calculated between KWIA-processed signal curves and those of the original images.

Quantification for dMRA was performed using an in-house MATLAB program of Block-circulant singular value decomposition method (cSVD) [21] with 15% threshold for singular value cutoff. A threshold that equals to three times of SD of background noise of each condition was used to differentiate vessel and background regions. The signal extracted from the MCA ROI was used as AIF. Subsequently, the aBF and ATT values in both PCA and distal small vessel ROI were estimated for all five subjects for comparison across all the calculated parametric maps.

KWIA processing of multi-delay ASL

Post-processing was performed offline using custom MATLAB scripts using SPM12 [22]. Motion correction was performed for 14 pairs of label-control images, and image co-registration was applied between MPRAGE and the first M0 volume. GM and WM masks were generated by segmentation of MPRAGE images using SPM12. To demonstrate the accuracy of perfusion measurements in small brain structures, the choroid plexus (CP) and caudate nucleus (CN) were manually segmented from MPRAGE images using ITK-SNAP [23]. Perfusion weighted images (PWI) were derived from the pairwise subtraction between label and control volumes. KWIA was applied to the magnitude images of PWI with the configuration of N = 4 and q = 15/48 to achieve a SNR ratio of 2. As the simulation results in Section 4.1 and the in-vivo result in Section 2 of online supplementary materials demonstrated that ASL is relatively insensitive to the temporal error caused by KWIA, this KWIA configuration was selected to maximize SNR ratio with minimal temporal error. SNR was measured by dividing the mean PWI signal measured in GM by the SD (corrected for Raleigh distribution) measured in the background region and averaged from 4 to 11 time frame that had a full KWIA window. Four PWI signal curves were extracted from GM, WM, CP, and CN, respectively, and the NRMSE was computed between KWIA-processed and the original PWI signal curves to analyze the temporal error caused by KWIA.

CBF, ATT, aCBV, aBAT, and residual maps were derived using an in-house MATLAB program of two-compartment model fitting for 14-delay pCASL [18] (see online supplementary materials section 1). The values of CBF and ATT were extracted from GM and WM, while the values of aCBV and aBAT were measured from an ROI of macrovascular component. Also, the SD values of aCBV and aATT maps were measured from a hand-drawn ROI in WM region without visible vasculatures.

Two-tailed paired t-test was carried out to compare the SNR and all quantitative parameters before and after applying KWIA for dMRA and multi-delay ASL data, respectively.

Comparison with conventional denoising methods

KWIA was compared with conventional denoising methods, including weighted moving average, 2D Gaussian filtering in spatial domain, and BM3D for ASL-based dMRA and multi-delay ASL using the same data described in Section 3.2. KWIA was applied with the same configurations as those described in previous sections which were N = 3 and q = 33/88 for dMRA to achieve 1.73-fold SNR, and N = 4 and q = 15/48 to achieve 2-fold SNR for multi-delay ASL. The weights of weighted moving average were determined according to the expected SNR ratios achieved by KWIA. Besides, the sliding window width of weighted moving average that minimizes the temporal error was determined by evaluating various sliding window widths from 3 to 7 for dMRA and 4 to 9 for multi-delay ASL (see Section 3 in online supplementary materials). Specifically, sliding window widths of 3 and 7 were chosen for dMRA and multi-delay ASL, respectively. The weight for the time frame at the sliding window center was 0.33 for dMRA and 0.48 for multi-delay ASL, while the weight for the neighboring frames was 0.33 for dMRA and 0.09 for multi-delay ASL. The 2D Gaussian filtering was implemented using the MATLAB image processing toolbox with the standard deviation leading to the 1.73-fold ratio for dMRA and 2-fold SNR ratio for multi-delay ASL. For BM3D, the BM3D MATLAB toolbox [24] was used with the required parameter, noise SD, measured from a background region of angiography images for dMRA and PWIs for multi-delay ASL.

Results

Simulation

Fig. 2 shows the digital phantom images processed by KWIA with six sets of configurations applied on simulated noisy magnitude images with the target to improve SNR by 1.41 (N = 2 and q = 36/64, and N = 3 and q = 42/64), 1.73 (N = 3 and q = 25/64, and N = 4 and q = 30/64), and 2 times (N = 3 and q = 11/64, and N = 4 and q = 20/64), respectively. For comparison, the original phantom images with lower added noise SD (1.41-fold SNR, 1.73-fold SNR, and 2-fold SNR) are displayed as references. Compared to the original noisy images, all configurations of KWIA processing were able to suppress the noise to a level visually comparable to the noise level presented in the corresponding reference images with lower noise SD.

The measured SNR ratios (listed in Table 2) in the simulation study under seven different MRI acquisition conditions generally matched the predicted values, except for k-space filtering and zero-filling, where KWIA underperformed by about 10–15%. This is likely due to either reduced noise (low-pass filtering) or no noise in the outer k-space of the original image, which is expected with k-space filtering or zero filling on the KWIA weighting function and is not taken into account in Eq. (5). Specifically, with the designed Fermi window, the predicted SNR ratios reduced from 1.41, 1.73, and 2 to 1.22, 1.46, and 1.75, respectively. With the 50% zero filling applied along the PE direction, the predicted SNR ratios reduced from 1.41, 1.73, and 2 to 1.31, 1.58, and 1.84, respectively. The revised values of predicted SNR ratios were more consistent with the results obtained in the simulations.

Table 2

Measurement of SNR ratios in simulation and in-vivo experiments under various MRI acquisition conditions.

	Expected SNR ratio	1.41		1.73		2
	KWIA config	N = 2	N = 3	N = 3	N = 4	N = 3	N = 4
Simulation	Standard	1.41	1.42	1.74	1.74	2	2
	Rectangular FOV	1.4	1.42	1.73	1.72	1.96	1.96
	Low-pass filtering (Fermi window on PE direction)	1.22	1.19	1.46	1.45	1.75	1.71
	Zero filling (50% on PE direction)	1.31	1.28	1.58	1.54	1.84	1.79
	Partial FT (75% on PE direction)	1.37	1.38	1.65	1.68	1.95	1.96
	SENSE (R = 2 in PE direction)	1.42	1.41	1.74	1.75	1.99	1.99
	GRAPPA (R = 2 in PE direction)	1.41	1.41	1.73	1.71	1.92	1.93
In-vivo dMRA	Partial FT (75% on both PE and slice directions			1.62
	& In-plane GRAPPA (R = 2))			±0.01
In-vivo multi-delay ASL	Standard						1.96
							±0.05

Fig. 3(a) displays the measured NRMSE values under six KWIA conditions for evaluating the temporal error, and Fig. 3(b) shows comparisons of the original and KWIA-processed time curves with targeted SNR ratio of 1.41, 1.73 and 2 respectively. With a SNR ratio of 2, KWIA introduced a minimal temporal error (NRMSE = 0.22%) for GM/WM ROIs and a larger temporal error (NRMSE = 8%) for small vessel ROIs, suggesting that KWIA can be applied to ASL to achieve 2-fold SNR without severe temporal error, while the trade-off between SNR ratio and temporal error should be considered for dMRA using conservative KWIA configurations. For dMRA, with SNR ratios of 1.41 and 1.73, KWIA with smaller N caused lower temporal errors for all vessels ROIs, indicating smaller N is preferred. Specifically, the NRMSE in dMRA caused by KWIA can be reduced to 1.7%, 2.2%, and 4.7% from 2.9%, 3.7%, and 8% for large, medium, and small vessel, respectively, using N = 3 with a SNR ratio of 1.73 compared to N = 4 with a SNR ratio of 2. For ASL, at 2-fold SNR ratio, KWIA with N = 3 led to higher temporal errors for GM and WM compared to KWIA with N = 4. Therefore, 2-fold SNR is achievable for KWIA with N = 4 in ASL applications, while KWIA with N = 3 and 1.73-fold SNR is more desirable for dMRA applications. Similar results were obtained using in vivo data as shown in Section 2 of online supplementary material.

Fig. 3.

(a) Bar plots showing NRMSE between the original and 6 KWIA-processed signals measured in 3 dMRA vessels (first row) and ASL GM/WM regions. (b) Visualization of the temporal error caused by KWIA. 3 KWIA configurations (1.41-fold SNR (N=2, q=36/64), 1.73-fold SNR (N=3, q=30/64), and 2-fold SNR (N=4, q=20/64)) were selected. The temporal errors in GM, WM are negligible, while dMRA vessels seem more susceptible.

In-vivo dMRA

Fig. 4 shows 12 time frames of dMRA maximum intensity projection (MIP) images of a representative subject with and without 1.73-fold-SNR KWIA (N = 3, q = 33/88).). The insets show two enlarged ROIs indicated by the blue and red boxes. With improved SNR, KWIA revealed small distal vessels (red arrows) that were overwhelmed by noise in the original dMRA. This enhancement of distal vessels is a result of improved SNR instead of temporal blurring of neighboring frames because no such clear structures can be visually captured within the KWIA window (five time frames) in the original dMRA images. Better visualization of dynamic flow can be observed in the animation in Section 4 of online supplementary materials. The measured SNR ratio (1.62±0.01) in MCA listed in Table 2 indicates significant SNR improvement was achieved by KWIA (p < 0.0005), and the measured SNR ratio matched the predicted ratio (1.73) although rectangular FOV, partial FT, and GRAPPA were applied during dMRA data acquisition.

Fig. 4.

The entire time series of ASL-based 4D dMRA images with and without 1.73-fold-SNR KWIA (N=3 and q=33/88). Two insets display the enlarged ROI in red and blue boxes. The distal vasculatures (pointed by red arrows) overwhelmed by noise can be recovered by KWIA, suggesting strong enhancement of dynamic visualization by KWIA.

Fig. 5(a) displays dMRA signals of a representative subject extracted from MCA, PCA, and distal small vessel ROIs. The KWIA-processed signal curves closely follow those of the original dMRA with an increasing temporal error from MCA to PCA to distal vessels. On average, as shown in the bar plot in Fig. 5(a), the mean NRMSE values for MCA, PCA, and distal vessels across five subjects are 1.5±0.9%, 2.5±0.9%, and 5.5±3.5%.

Fig. 5.

(a) Comparison between the original and the KWIA-processed dMRA signal curves of MCA, PCA, and distal small vessel from a representative subject, showing that all KWIA signals closely follow the original signals. The bar plot shows that small temporal errors were caused by KWIA. (b) Comparison between the original and the KWIA-processed PWI signals curves in GM, WM, CP, and CN. KWIA signals were almost identical to the original signals with NRMSE about 1%.

Fig. 6 shows quantitative parametric maps of KWIA-processed and original dMRA from two representative subjects, and the following analysis is based on the mean quantitative parametric values across all subjects. For large vessels, no significant change was observed in aBF maps with an 2% increase (p = 0.38), while ATT was slightly higher (7.5% increase, p = 0.14) in KWIA processed maps compared to the original maps. In comparison, KWIA effectively recovered the missing distal vessels in original aBF and ATT maps without severely affecting the value of quantitative flow parameters (8% increase, p = 0.06 for aBF, and 9% increase, p < 0.05 for ATT) for distal small vessels. Quantitative parametric values are listed in online supplementary materials Table S2.

Fig. 6.

Both aBF and ATT maps from KWIA-processed dMRA show better coverage of the distal vasculatures compared to those of the original noisy images. The quantitative flow parameters were not severely affected by KWIA.

In-vivo multi-delay ASL

Fig. 7 shows dynamic perfusion images of a 14-delay pCASL with and without 2-fold-SNR KWIA (N = 4, q=15/48) respectively. The insets show the zoomed images at six selected PLDs to highlight the SNR difference. It can be observed that KWIA processed perfusion images have higher SNR, especially in low perfusion regions such as WM, and a higher contrast between GM and WM compared to the original images. While the perfusion signals in the CP and CN regions (indicated by red and yellow arrows, respectively, in Fig. 7) were severely contaminated by noise in the original images, KWIA enabled the visualization of dynamic perfusion signals in these small brain structures (see the animations in Section 4 of online supplementary materials). In the meantime, no apparent spatial and temporal error or smoothing was introduced. The measured SNR ratio of KWIA was 1.96±0.05 (listed in Table 2), indicating that the SNR of KWIA-processed PWI was significantly greater than that of original PWI (p < 0.0005). This SNR improvement can translate into reduced temporal fluctuation as shown by the results in Section 4 of online supplementary materials. As shown in Fig. 5(b), the PWI signals measured in GM and WM were only slightly affected by KWIA with NRMSE of 0.3±0.1% and 0.9±0.3%, respectively. The NRMSE of the PWI signal measured in CP and CN region was 4.7±1.8% and 2.1±0.9%, respectively.

Fig. 7.

Two-fold SNR ratio was achieved by KWIA (N=4 and q=15/48) resulting in better image quality and dynamic visualization. Strong enhancement of dynamic visualization in CP and CN regions (indicated by yellow and red arrows, respectively) can be observed. No obvious degradation of spatial and temporal resolution was found.

Fig. 8 displays the quantitative parametric maps of multi-delay ASL from two representative subjects, and the mean quantitative parametric values across all subjects are reported in the following analysis. The CBF and ATT maps of KWIA-processed data were visually similar to that of the original data for all subjects. After processing with KWIA, the CBF and ATT values measured in both GM and WM, were slightly lower than the original values with small differences of 0.5% (p = 0.02) for GM CBF, 0.6% (p = 0.02) for WM CBF, 2.2% (p = 0.003) for GM ATT, and 1.4% (p = 0.002) for WM ATT. In aCBV and aBAT maps, the macrovascular structures were not obviously affected by KWIA, while reduced noise (grainy signals) can be observed in the tissue regions without macro-vasculatures. The differences of aCBV and aBAT measured in macro-vasculatures were 5.8% (p = 0.01) and 0.2% (p = 0.4), respectively, and the differences of the SNR values measured in tissue regions were 6.6% for aCBV (p = 0.06) and 19.9% (p = 0.01) for aBAT. Moreover, the residual of model fitting was significantly reduced by 2.7 times (p = 0.0001) using KWIA. Quantitative parametric values and SNR values are listed in online supplementary materials Table S3 and S4.

Fig. 8.

Quantitative parametric maps (CBF, ATT, aCBV, aBAT, and residual) for multi-delay ASL of two representative subjects under original and KWIA (2-fold SNR (N=4, q=15/48)), respectively. Visually comparable CBF, ATT, aCBV, and aBAT maps were generated from the original and KWIA-processed ASL. A significantly suppressed residual of model fitting was observed presumably due to the improved SNR.

Fig. 9 displays the MIP images of dMRA processed by KWIA and that processed by weighted moving average, Gaussian filtering, and BM3D on one representative subject. KWIA achieved a similar SNR level to weighted moving average. However, the latter caused more distal vessels that can be visually captured potentially due to a more servere temporal averaging effect compared to KWIA. A higher NRMSE (8.6±2.1%) was measured in the distal vessels acorss five subjects compared to the NRMSE reported in Section 4.2 for KWIA. Moreover, the evolution of the signal intensity in both PCA and MCA were also more severely smoothed by weighted moving average compared to KWIA with a NRMSE of 7.5±2.2% and 5.8±2.0% measured in PCA and MCA, respectively, which were two to five times higher than the NRMSEs reported for KWIA.

Fig. 9.

Comparison of the denoised MIP images of ASL-based dMRA processed by KWIA with that processed by weighted moving average, Gaussian filtering, and BM3D.

Fig. 10 shows the PWIs of multi-delay ASL processed by KWIA and that processed by weighted moving average. Although KWIA and weighted moving average removed noise from the original PWIs to a similar degree, the signal change in the PWIs processed by weighted moving average did not follow the original signal change as close as the PWIs processed by KWIA did. Especially in the 9 time frame, where an abrupt signal drop happened in the original PWI, the signal intensity in KWIA-processed PWI decreased immediately. However, the PWI processed by weighted moving average failed to reflect such signal change. Quantitatively, NRMSEs of weighted moving average measured in GM and WM were 6.8±1.6% and 10.7±2.0%, respectively. Weighted moving average caused up to twenty times higher NRMSE than the NRMSE reported in Section 4.3 for KWIA.

Fig. 10.

Comparison of the denoised PWIs of multi-delay ASL processed by KWIA with that processed by weighted moving average, Gaussian filtering, and BM3D.

The results of the comparison of KWIA to Gaussian filtering and BM3D are also presented in Fig. 9 and Fig. 10 for dMRA and multi-delay ASL, respectively. Comparing to KWIA, at the same denoising level, Gaussian filtering resulted in more blurring of small vessels in dMRA and a blurred contrast between GM and WM. Despite that a higher SNR can be achieved by BM3D compared to KWIA, BM3D is relatively undesirable for dMRA and multi-delay ASL because of several disadvantages. First, small features and features with low signal intensity could be lost due to the block matching of BM3D. For example, some distal vessels were missing in every time frame of dMRA, and the signal intensity in WM of the PWIs of ASL was significantly reduced by BM3D compared to the original one. Second, the computational time of BM3D is much longer than KWIA. Under the same execution environments (Intel i5-9400F, 16G RAM), it took 175 s for BM3D to process a time series of a dMRA volume data, while only 4 s was needed for KWIA.

Discussion

In this study, we demonstrated that KWIA improved SNR of multi-delay ASL and ASL-based dMRA. There are several potential advantages of KWIA compared to conventional denoising methods: (1) KWIA can be easily implemented in clinical settings since KWIA works in the k-space domain which can be simply converted from MR DICOM images by FFT (i.e., original raw k-space data is not required); (2) KWIA is computationally simple and fast (non-iterative); (3) the spatial and temporal resolution of the original MR images are well-preserved. Compared to conventional view-sharing techniques such as KWIC or sparse sampling techniques such as compressed sensing that prospectively undersample k-space, KWIA can be retrospectively applied to MRI data acquired with various acquisition conditions and acceleration methods. Furthermore, KWIA may be combined with these prospective acceleration methods to increase SNR.

KWIA is a flexible algorithm that allows a variety of configurations, such as different weighting functions and ring shapes. Besides the discrete weighting function demonstrated in this work, a linearly or nonlinearly varying weighting function may be applied in KWIA that continuously decreases from the central to outer k-space. This linear or nonlinear weighting function would allow more efficient SNR gain compared to the discrete version. However, severe temporal blurring could happen as the result of excessive averaging in central k-space. A square ring shape is an alternative to the circular ring. Nevertheless, since the frequency distribution in k-space follows Euclidean distance, circular ring is a more intuitive selection. Moreover, the circular ring more efficiently improves SNR since more high frequency k-space components (4 corners) are averaged by the outmost ring.

The benefit of KWIA in quantification metrics needs to be further explored in more dynamic MRI applications and by using different quantification algorithms. In this study, KWIA either didn’t or only slightly alter the quantification metrics (<10% change) which is within the measurement precision of ASL techniques [25]. For the curve fitting algorithm used in ASL, noise can be removed in the fitted curve by minimizing the mean squared error, which makes the quantification result insensitive to the change of noise level (unless there is a local minimum). For the SVD method used in dMRA, however, the presence of substantial noise often leads to the overestimation of aBF [7]. While KWIA is expected to correct the bias by improving SNR, such effect was not observed potentially due to the cancellation between the benefit of higher SNR and the consequence of temporal blurring. Nevertheless, we observed ~10% SNR increase in aCBV and aBAT maps for multi-delay ASL, and improved visualization of distal small vessels in aBF and ATT maps for dMRA. Overall, our data suggest that KWIA is valuable for denoising and dynamic visualization without adverse effects on quantification.

Several denoising methods have been proposed for ASL, such as image domain filtering BM3D method [26], Non-Local Mean [27], Wavelet-based denoising [28], total generalized variation (TGV) based regularization [29], and deep learning based methods [30,31]. In this work, We compared KWIA with conventional denoising methods including weighted moving average, Gaussian filtering, and BM3D. The results showed that KWIA kept a better temporal fidelity compared to weighted moving average, and was able to preserve temporal and spatial resolution better than Gaussian filtering and BM3D. In addition, a TV-based iterative reconstruction method has been compared with KWIA for CT perfusion, and the results showed a stronger denoising effect but with greater spatial blurring and 100× longer computation time using the former method [13]. Although the above study was performed on CT perfusion, a greater degree of spatial blurring is expected with spatiall regularized denoising methods (e.g., TGV), while a longer computation time is expected with iterative denoising methods (e.g., TGV, NLM, BM3D) compared to KWIA on dynamic MRI. Deep learning-based methods can be executed instantaneously, however the model relies on the available training data. While KWIA was proposed as a post-processing method with magnitude images as the input, it can also be directly applied to the raw k-space data obtained from MRI scanner without the conversion from image to k-space at the beginning of the algorithm. The former is named image-based KWIA, and the latter is called ksp-based KWIA. Compared to ksp-based KWIA, image-based KWIA yields identical performance in denoising but does not correct the bias (1.26× noise SD) in the background region (see online supplementary materials section 8). With the removal of such bias, the performance was similar between image-based and ksp-based KWIA. However, for ASL-based dMRA and multi-delay ASL, complex subtraction between the k-space of label and control images may result in phase errors that requires further investigations. Applying KWIA on the k-space of label and control images separately before subtraction is also prohibited because of the varying intensity across different PLDs caused by background suppression. For other dynamic imaging techniques involving no subtractions, KWIA may be applied to k-space samples of each channel first, and then the reconstruction can be done with parallel imaging that involving multiple coils (see section 9 in online supplementary materials). Besides improving SNR for images acquired with Cartesian trajectories, KWIA may be beneficial for non-Cartesian trajectories. For example, up to 2-fold SNR improvement was achieved for the radial sampled CT perfusion. Although the use of KWIA with radial trajectory was not demonstrated for MRI in this work, CT perfusion has the same concept of radial k-space after Radon transform. Also, KWIA might be compatible with the reconstruction methods that utilize k-space samples from the full extent of the k-space in a time frame, such as iterative reconstruction and compressed sensing. However, KWIA is not suitable for the reconstruction methods with specific accelerated acquisition, such as Keyhole-based dynamic imaging, because of the replicated information in outer k-space. In addition to the presented applications, KWIA is applicable to other dynamic MRI techniques such as dynamic susceptibility contrast (DSC) and dynamic contrast enhanced (DCE) MRI, as well as functional MRI. In the future, the application of KWIA may be generalized to more dynamic MRI techniques or any time series.

Despite the benefit of SNR improvement and simplicity, KWIA has a few limitations. Firstly, KWIA is sensitive to motion since multiple frames are combined during filtering. Images are recommended to be motion-corrected before processed by KWIA. Secondly, KWIA suffers from the potential temporal error caused by averaging between neighboring frames, especially for the fine vessels and/or structures. For ASL perfusion applications, this temporal error is negligible (NRMSE about 1%) for GM and WM, even with the highest (2-fold) SNR ratio because of the low spatial resolution and gradual change in PWI signals. Although there is a stronger temporal blurring effect on the smaller brain structures, such as CN and CP, KWIA enables enhanced visualization of dynamic perfusion signals in those regions, which is difficult to identify in the original images due to noise. However, for dMRA which focuses on finer structures, the temporal error cannot be ignored. Consequently, moderate SNR ratio (1.73-fold) is more suitable for dMRA with NRMSE less than 5%. Therefore, different KWIA configurations can be tuned empirically according to the specific goal of the application. Thirdly, slight texture change was observed in the uniform region of the digital phantom in Fig. 2 and Fig. S3 potentially due to KWIA filtering of the noise across spatial frequencies. Such texture change neither appeared in in-vivo data nor affected the SSIM and PSNR measurements (see online supplementary materials section 8), but further investigation in some extreme cases is required. Finally, the noise property of images processed by KWIA may be altered because of the nonuniform distribution of noise level across k-space and non-Gaussian nature of magnitude images. Despite that no negative impact was observed in this study, the potential effect of KWIA on lesion detection and quantification requires further investigation.

Conclusion

In this study, a novel denoising algorithm termed KWIA is presented for dynamic MRI such as ASL-based dMRA and multi-delay perfusion imaging. Both simulation and in-vivo experiments showed that KWIA can achieve a SNR ratio of 1.73 for dMRA and 2.0 for multi-delay ASL without markedly compromising the spatial and temporal resolution or quantification accuracy. Such significant SNR improvement enables improved visualization of flow dynamics in distal small vessels in dMRA and fine brain structures in multi-delay ASL. With its computational simplicity and ease of implementation in clinical scenarios, KWIA may provide an efficient denoising method for dynamic MRI.