Parallel imaging reconstruction using spatial nulling maps

To develop a robust parallel imaging reconstruction method using spatial nulling maps (SNMs).


INTRODUCTION
Parallel imaging (PI) is widely used in clinical MRI to accelerate data acquisition or correct artifacts through the use of multiple receiving coils where each coil exhibits a unique spatial coil sensitivity map (CSM). [1][2][3][4] PI reconstruction techniques generally fall into three classes: image-domain, 4-8 k-space, [9][10][11][12][13][14][15][16][17][18][19][20] and hybrid-domain category, 21 depending on where the coil sensitivity information is used to perform the reconstruction. The image-domain methods such as SENSE 4 use explicit knowledge of the CSMs. In theory, they can produce optimal least-square reconstruction if accurate CSMs are available. [21][22][23] The image-domain methods also provide the flexibility of incorporating image priors such as image sparsity to suppress noise amplification (e.g., through image-domain regularization). [24][25][26] However, accurate CSMs are sometimes difficult to obtain in practice and their inaccuracies can cause noise amplification and artifacts in the reconstructed images. 23,27 The k-space reconstruction methods (e.g., GRAPPA and SPIRiT) 9,18 use local linear dependencies among k-space samples for reconstruction. They provide robust reconstruction through coil sensitivity autocalibration data, mitigating the problem of spatial mismatch between the coil sensitivity calibration and undersampled data in SENSE reconstruction (e.g., because of motion.) However, the k-space reconstruction methods cannot easily incorporate image priors for further improving reconstruction performance. Among the k-space methods, parallel reconstruction using null operations (PRUNO) is a subspace based reconstruction method in which a k-space nulling system is formed by null-subspace bases of the calibration matrix. 16 PRUNO produces fewer residual artifacts with limited autocalibration signal when compared with GRAPPA. 16 However, as shown in the original ESPIRiT study, 21 PRUNO is computationally prohibitive in practice, and it cannot offer the flexibility of incorporating image priors as in image-domain reconstruction methods. As a hybrid-domain method, ESPIRiT extends the concept of subspace based reconstruction developed in PRUNO. It provides a more efficient and flexible approach, in which the signal-subspace bases are extracted to derive coil sensitivity information for an extended SENSE based reconstruction. 21 By obtaining coil sensitivity information from subspace bases through eigendecomposition, ESPIRiT bridges the gap between image-domain SENSE and k-space GRAPPA while retaining the benefits of both. However, ESPIRiT requires manually thresholding eigenvalues to mask the eigenvector maps to spatially define image background region to minimize the noise propagation during image reconstruction. Therefore, the optimal threshold is image-content dependent and often involves a trial-and-error selection. A sub-optimal threshold can undermine the reconstruction quality (e.g., causing noise increase and residual artifacts.) In addition, ESPIRiT reconstruction is sensitive to the subspace cut-off selection, through which signal-subspace must be accurately separated from null-subspace. Inaccurate subspace division can lead to inaccurate coil sensitivity information and degrade the reconstruction performance. Stein's unbiased risk estimate (SURE) 28 was recently developed to select the ESPIRiT's parameters automatically. It searches for the optimal parameters with minimum mean-square-error reconstruction in a data-driven manner and therefore, avoids manually optimizing eigenvalue thresholding and subspace cut-off. However, this method is computationally intensive although it uses strategies to alleviate the computational burden. Furthermore, because the optimal parameters are image content dependent and need to be estimated for each single slice, SURE-based ESPIRiT reconstruction can become very time-consuming in practice, especially for dynamic imaging reconstruction.
In this study, we combine the concepts of PRUNO and ESPIRiT to present a robust and flexible hybrid-domain based reconstruction method. The proposed method extracts null-subspace bases of the calibration matrix from k-space coil calibration data, and uses them to derive image-domain spatial nulling maps. The subsequent image reconstruction relies on solving the nulling system formed by spatial nulling maps. However, this method does not entail any masking-related procedure and is relatively insensitive to the subspace cut-off selection. It eliminates the need for any intensive parameter fine-tuning, therefore, offering a more efficient alternative to the SURE-based ESPIRiT.

Existing PRUNO and ESPIRiT methods
PRUNO is a k-space PI reconstruction method in which a nulling system is derived using null-subspace bases. 16 A block-wise Hankel matrix called calibration matrix A is constructed from k-space coil calibration data ( Figure 1A). The subspace bases are obtained by applying singular value decomposition (SVD) on the calibration matrix and selecting a cut-off to separate the signal-subspace V || and null-subspace V ⊥ . The columns in V || and V ⊥ are signal-and null-subspace bases, respectively. For each multi-channel k-space local subset d in the calibration data, the k-space nulling system can be formed as ( Figure 1B): (A) Signal-and null-subspace separation. Central fully-sampled multi-channel k-space or coil sensitivity calibration data are vectorized to block-wise Hankel matrix A. Singular value decomposition is then performed on A to calculate singular vector matrix V. By setting a cut-off in V, signal-subspace V || and null-subspace V ⊥ are separated. (B) Parallel reconstruction using null operations (PRUNO) reconstruction. Each null-subspace basis v in V ⊥ is segmented and devectorized into 2D null-subspace convolution kernels f null ij forming the k-space nulling system matrix F. The missing k-space can be estimated by solving the nulling system equation F H d = 0. (C) ESPIRiT reconstruction. f sig ij is obtained from the signal-subspace basis in V || and transformed into an image-domain map s sig ij , which can form a large 2D overdetermined signal system matrix S sig . Pixel-wise eigen-decomposition is then performed on the S sig to calculate eigenvector maps (unmasked ESPIRiT maps) and eigenvalues. Only a single set of eigenvector maps is displayed. After a manual thresholding process, an eigenvalue mask is obtained to mask the eigenvector maps. With masked ESPIRiT maps E, an extended SENSE-based reconstruction can be applied. (D) Reconstruction using spatial nulling maps (SNMs). The null-subspace basis v in V ⊥ is used to obtain f null ij and successively transformed into s null ij , forming a large overdetermined nulling system matrix S null . Instead of solving the overdetermined system, all s null ij are combined to construct multi-channel SNMs denoted by N. Images can be then reconstructed by solving the system established using N and multi-channel images m.
(1) Here, F is the nulling system matrix and F H the Hermitian matrix of F. f null ij represents each 2D null-subspace convolution kernel that is transformed from the th null-subspace basis v (the th column in V ⊥ ) through devectorization of the ith channel segment. n c is the number of channels and n the number of null-subspace bases. The nulling system is essentially composed of multiple 2D nulling convolutions and summations within multi-channel k-space. Missing k-space data can be reconstructed by solving this nulling system across the whole k-space with an iterative conjugate-gradient algorithm. In GRAPPA, as acceleration factor increases, the correlation between target samples and fitting samples becomes weaker, therefore, requiring larger kernels to fit more samples. In contrast, the nulling system in PRUNO always continuously applies the local convolutions at all k-space locations regardless of the reduction factor. 16 With iterative reconstruction, all missing data will eventually be fitted through using all available k-space samples, during which the nulling kernel yields better fitting accuracy because of the larger number of fitted samples and better locality. 16 However, the iteration process in PRUNO reconstruction is computationally demanding and PRUNO does not offer the flexibility of incorporating image priors as in image-and hybrid-domain methods. 21 ESPIRiT extends the subspace concept in PRUNO further and provides a more computationally efficient hybrid-domain approach in which an extended SENSE based reconstruction can be performed. 21 As shown in Figure 1C, the signal-subspace bases are extracted from the calibration matrix to form signal-subspace based shift-invariant convolutional kernels, which differ from the null-subspace based nulling kernels in PRUNO. Each k-space kernel f sig ij is transformed into an image-domain map s sig ij through zero-padding denoted by Z and inverse fast Fourier transform (IFFT) denoted by F −1 : where l is the number of signal-subspace bases. Next, an overdetermined signal system matrix S sig can be formed in image domain: Here, m denotes multi-channel images. Masked ESPIRiT maps E that contain coil sensitivity information are then obtained by (1) performing pixel-wise eigen-decomposition on the system matrix to use its linear relationship to the coil sensitivities and (2) masking the eigenvector maps using a manually chosen eigenvalue threshold to exclude image background region so to minimize noise propagation during image reconstruction. 21 With the masked ESPIRiT maps E, SENSE based reconstruction can be performed to reconstruct the target image by using such estimated coil sensitivity information in E: Here, m 0 is the underlying target image and m a are the aliased multi-channel images reconstructed from undersampled k-space data. F is the fast Fourier transform (FFT) and P an operator that undersamples k-space. The multi-channel images m can then be written as: Although ESPIRiT extends the subspace concept in PRUNO while preserving the image-domain reconstruction flexibility, it requires the manual eigenvalue thresholding for masking coil sensitivity information, which is a cumbersome procedure in practice because such thresholding is image-content dependent. Further, ESPIRiT reconstruction quality can be sensitive to the arbitrary subspace division during signal-subspace extraction process.

Proposed hybrid-domain reconstruction using spatial nulling maps
In this study, we integrate the concepts of the null-subspace reconstruction in PRUNO and the hybrid-domain reconstruction in ESPIRiT and present a new and more robust PI reconstruction method. The proposed method extracts null-subspace bases of k-space calibration matrix and directly derives a set of image-domain spatial nulling maps from these null-subspace bases ( Figure 1D). The spatial nulling maps explicitly represent both coil sensitivity and finite image support information and form an image-domain nulling system for subsequent image reconstruction. The specific reconstruction steps are described below.
First, the central consecutive fully-sampled k-space lines within the multi-channel k-space data are used to construct a block-wise Hankel calibration matrix A, in which the column entries (i.e., vectorized k-space blocks) exhibit strong linear dependencies 9,16,17,21 ( Figure 1A). By performing the SVD, the singular vector matrix V with signal/null-subspace bases of A can be obtained: where V H represents the Hermitian matrix of V. By setting a cut-off, the signal-subspace spanned by V || and null-subspace spanned by V ⊥ are separated. Assuming that the two subspaces are ideally separated, two constraints are satisfied: With extracted V ⊥ , the segment corresponding to the ith channel of each null-subspace basis v is transformed into a 2D null-subspace convolution kernel f null ij through devectorization. According to Equation (8), a convolutional nulling relation between f null ij and multi-channel k-space data can be established: where k i denotes the k-space data from the ith channel. In contrast to PRUNO, each k-space kernel f null ij is transformed into an image-domain map s null ij through zero-padding and IFFT: Using s null ij , an image-domain overdetermined nulling system can be formed: Here, S null is a large 2D overdetermined nulling system matrix. Instead of solving the overdetermined system in Equation (11) for image reconstruction, multiple s null ij are combined to construct multi-channel spatial nulling maps N: where s null ij denotes the conjugation of s null ij . With N, an image-domain nulling system can be built: Here, m can be further represented as m = m a + m m , where the m m denotes the multi-channel images corresponding to the underlying missing k-space data. Therefore, the image-domain nulling system in Equation (13) becomes: The multi-channel images m (i.e., m a + m m ) can then be reconstructed by solving this nulling system. Specifically, with spatial nulling maps N estimated from the central k-space lines, a least-square solution of m m in Equation (14) can be calculated 29 : Here, ‖PF m m ‖ 2 2 is the data consistency term. Further, as in SENSE and ESPIRiT, the reconstruction can readily incorporate regularization terms Ψ to further reduce image noise by forming: where is the regularization weight.

Finite image support information in spatial nulling maps
Both our proposed method and ESPIRiT extend the subspace notion of PRUNO and provide computationally more efficient and flexible hybrid-domain reconstruction. However, unlike ESPIRiT, spatial nulling maps in our proposed method contain both coil sensitivity and finite image support information.
The subspace bases that satisfy the constraints in Equations (7) and (8) essentially represent the underlying low-frequency modulations in MR data, which can be indicated as k-space convolutions with limited kernel size or image-domain multiplications by smooth varying maps. Two dominating low-frequency modulations for multi-channel MR images are coil sensitivity and finite image support because of their nature of slow spatial variation. Note that such finite image support information is related to the finite boundary of the object within the FOV, 30 and it has been used for parallel imaging reconstruction. [31][32][33] In presence of finite image support (i.e., when imaging FOV is larger than the object), a nulling relationship inherently exists between the finite image support p and MR images: Considering Equation (11), the complement of the finite image support (i.e., 1 − p) here is indeed embedded within the null-subspace S null of m in the image domain. Note that, in k-space, such finite image support information, that is, F (1 − p), is represented within the null-subspace convolutional kernels described in Equation (9). Therefore, the spatial nulling maps N computed as in Equation (12) contain both coil sensitivity information and finite image support information. Therefore, our proposed method enables image reconstruction that involves no explicit masking-related procedure (see Figure 1C vs. Figure 1D).

Insensitivity to subspace separation
The proposed method is relatively tolerant of inaccurate subspace separation. In ESPIRiT, the subspace separation determines the number of signal-subspace bases being extracted, which dominates the accuracy of S sig and needs to be fine-tuned ( Figure 1C). Therefore, when performing pixel-wise eigendecomposition on S sig , the coil sensitivity information is affected by the subspace separation. In our proposed method, the null subspace separation determines the number of null-subspace bases and hence influences the accuracy of spatial nulling maps. However, for multi-channel MR data, the number of signal-subspace bases is much smaller than the number of null-subspace bases (i.e., l < n) because of the redundancy and linear dependency within MR data. Therefore, when the subspace separation is inaccurate, that is, extracting relatively excessive or insufficient number of signal/null-subspace bases, the null-subspace spanned by n bases is overall less affected than the signal-subspace spanned by l bases.

Evaluation of the proposed method
The proposed method was evaluated and compared to the ESPIRiT. For ESPIRiT reconstruction, the eigenvector maps that contain coil sensitivity information were termed ESPIRiT maps. Both spatial nulling maps and ESPIRiT maps were estimated from the central consecutive k-space lines. An L 1 -norm wavelet sparsity regularization with the regularization weight = 0.002 was also applied to both the proposed method and ESPIRiT. The kernel size for both the proposed method and ESPIRiT was set to 6 × 6. Two sets of eigenvector maps were applied in ESPIRiT for the extended SENSE based reconstruction. For both ESPIRiT and proposed method, the signal/null-subspace cut-off was set according to a manually determined 2 cut-off relative to the maximum singular value as in previous ESPIRiT. 21 In ESPIRiT reconstruction, the eigenvalue threshold for masking coil sensitivity information was manually optimized. They were 0.996, 0.96, and 0.95 for two brain datasets (with two very different head sizes) and one knee dataset, respectively. The performance with suboptimal threshold selection in ESPIRiT was also evaluated by deliberately setting the threshold smaller or larger than the optimal value.
The final reconstructed images were generated by combining all coil images using the square root sum-of-squares method. In addition to comparing the images reconstructed from the proposed method and ESPIRiT, the residual error maps were calculated as the square root sum-of-squares of the coil-by-coil difference between the reconstructed images and reference images reconstructed from the fully sampled data. Normalized root mean squared error (NRMSE) and structural similarity index measure (SSIM) 39 were calculated within the brain or knee region. L 1 -norm of sparsifying wavelet transform with different regularization weights ( = 0.001, 0.002, 0.003, and 0.004) was also applied for evaluation. In the implementation, the objective function minimized the joint-sparsity 18,40-42 of multiple channels, which differed from the channel-combined image sparsity in the ESPIRiT. 21 The proposed method was further evaluated using one 8-channel dataset with a FOV of 200 × 150 mm 2 that was smaller than the head size. A total of 24 central consecutive k-space lines of 256 lines were preserved.
Our proposed method and its evaluation were implemented using MATLAB R2020b (The MathWorks). All codes and test data can be obtained online (https://github. com/jiahao919/SNMs) or from the authors on request. Figure 2 displays the coil sensitivity and finite image support information in the spatial nulling maps (SNMs) from one 6-channel T1W brain dataset, together with the

F I G U R E 2
Coil sensitivity and finite image support information in spatial nulling maps (SNMs) from one 6-channel T 1 -weighted gradient echo dataset. Representative maps with clear spatial boundary are displayed for illustrating the property of SNMs. The corresponding unmasked ESPIRiT maps and masked ones via manual eigenvalue thresholding are also shown. Both SNMs and ESPIRiT maps exhibited coil sensitivity information, whereas the SNMs also provided additional information on the finite image support. For the proposed method, the complement of finite image support could be demonstrated by combined SNM. In unmasked ESPIRiT maps, the image object region could not be distinguished from the background (as shown by red arrows) and the combined unmasked ESPIRiT map was an identity or uniform map as expected. Only in masked ESPIRiT maps as well as their combined map, the finite support of the object was delineated. Note that the matrix size for SNMs and ESPIRiT maps was 256 × 218. They were computed from 24 central consecutive k-space lines.
corresponding unmasked ESPIRiT maps as well as masked ESPIRiT maps through manually optimized eigenvalue thresholding. It is evident that unmasked ESPIRiT maps and SNMs both provided coil sensitivity information. However, SNMs provided additional information on the finite image support. Note that, sum-of-squares combination of SNMs clearly delineated the image support boundary, indicating that finite image support information was indeed embedded in SNMs. In contrast, combination of unmasked ESPIRiT maps did not contain such information. Figure 3 compares the reconstruction using SNMs and manually optimized ESPIRiT at various acceleration factors (R = 2, 3, and 4). Our proposed method provided reconstruction results with noise level and residual artifacts that were highly comparable to ESPIRiT results at

F I G U R E 3
Reconstruction results of the 6-channel T 1 -weighted gradient echo data from two subjects with different head sizes (i.e., varying finite image support) at R = 2, 3, and 4. ESPIRiT reconstruction was manually optimized for the eigenvalue threshold selection and subspace separation to ensure its best performance. The proposed spatial nulling maps (SNMs) method provided reconstruction results that were highly comparable to ESPIRiT at low acceleration factors (R = 2 and 3). At R = 4, reconstruction using SNMs outperformed ESPIRiT visually, which was also evident from the normalized root mean squared error in percentage and structural similarity index measure values measured within the brain region. Error maps are displayed with enhanced brightness (×5).

F I G U R E 4
Comparison between our proposed method and the ESPIRiT reconstruction with different eigenvalue masking thresholds, ranging from 0.8 to 0.999, at R = 3 and 4. The results with manually optimized thresholds are shown with red boxes. Either larger or smaller threshold values strongly degraded ESPIRiT reconstruction performance. In contrast, our proposed spatial nulling maps method was free from such phenomenon, yet its performance was always highly comparable to the best ESPIRiT results. low acceleration factors (R = 2 and 3). As the acceleration factor increased (R = 4), reconstruction using SNMs slightly outperformed the ESPIRiT in terms of noise and artifact levels within the brain region as indicated by NRMSE and SSIM values. Note that, ESPIRiT reconstruction was conducted by carefully optimizing the eigenvalue threshold and subspace separation to ensure its optimal reconstruction performance. In contrast, our proposed method required no masking-related procedure because the finite image support information was already embedded in the SNMs. Figure 4 shows the reconstruction results of our proposed method versus the ESPIRiT reconstruction results that used different eigenvalue thresholds (0.8 ∼ 0.999) for

F I G U R E 5
Reconstruction results of our proposed method and ESPIRiT for 6-channel proton density-weighted knee data at R = 2, 3, and 4. ESPIRiT reconstruction was performed using different eigenvalue thresholds, ranging from 0.5 to 0.99. The proposed method led to less noise amplification or/and fewer residual artifacts than ESPIRiT in terms of normalized root mean squared error and structural similarity index measure values, especially at high acceleration factors. With threshold larger or smaller than the optimal value, ESPIRiT reconstruction exhibited severe artifacts and noise, especially at R = 3 and 4. masking coil sensitivity information. Such threshold selection clearly affected the ESPIRiT reconstruction performance in terms of noise or artifact levels. When the threshold was smaller than the optimal one, larger eigenvalue mask occurred. Consequently, the reconstructed images exhibited significant noise increase in both brain region and background as expected. Conversely, when the threshold was larger than the optimal one, mask size decreased, leading to severe artifacts and noise. Note that the optimal threshold was highly dependent of the image content here, mainly object size and shape. Therefore, it must be manually determined in a trial-and-error manner by judging the reconstruction quality. In contrast, the reconstruction using SNMs completely circumvented such

F I G U R E 6
Reconstruction results of our proposed method and ESPIRiT for 6-channel SSFP cardiac data at R = 2, 3, and 4. In ESPIRiT reconstruction, with a threshold larger than the optimal one, the eigenvalue masking became problematic. The signal void and edge regions were masked out, leading to severe residual artifacts and/or noise. Large noise amplification was also observed when using a threshold smaller than the optimal one. In contrast, the proposed method did not entail any masking procedure and robustly reconstructed images that were comparable to the best ESPIRiT results. masking-related procedure, demonstrating a more robust reconstruction than ESPIRiT.
The reconstruction results of 6-channel knee data using the proposed method and ESPIRiT with different thresholds are shown in Figure 5. The thresholds in ESPIRiT reconstruction ranged from 0.5 to 0.99. Eigenvalue masks in ESPIRiT were also displayed to show the influences of different thresholds. Threshold values smaller than the optimal value yielded larger masks, causing more noise amplification. Larger threshold values led to smaller masks and significantly increased artifacts. Meanwhile, the proposed method achieved overall low noise and artifact level when compared with ESPIRiT, especially at high acceleration. Figure 6 shows the reconstruction of the 6-channel cardiac data using the proposed method and ESPIRiT with different thresholds. Compared to the head and knee datasets, large signal voids and susceptibility artifacts were present within the abdominal region. As shown in these results, ESPIRiT reconstruction performance could be affected by threshold selection. When using a threshold larger than the optimal one, the low intensity region and the edge region were masked out, causing severe artifacts and/or noise. Large noise amplification was observed when using a threshold smaller than the optimal one. In contrast, the proposed method consistently achieved quality reconstruction results that were comparable to the optimized ESPIRiT results at R = 2 and 3, and slightly better than the optimized ESPIRiT result at R = 4.
To examine the effect of different subspace divisions, Figure 7 compares the reconstruction results of the proposed method and the ESPIRiT with varying σ 2 cut-off for one brain dataset at R = 3. The number of bases in the signal-subspace l and null-subspace n corresponding to

F I G U R E 7
Reconstruction results of the proposed method versus ESPIRiT results (using manually optimized eigenvalue threshold) with different 2 cut-off for subspace division at R = 3 for two 6-channel brain datasets. The proposed spatial nulling maps (SNMs) method clearly offered overall better reconstruction quality than ESPIRiT for each cut-off selection, demonstrating its relative insensitivity to inaccurate 2 cut-off selection. different 2 cut-off is shown in Table S1. Both the proposed SNMs and the ESPIRiT performance degraded at suboptimal 2 cut-off . However, the proposed method was overall less sensitive to the 2 cut-off selection. Specifically, the proposed method only showed a slightly increased noise level for extremely small 2 cut-off or residual artifacts for extremely large 2 cut-off and it could tolerate approximately 30% deviation from optimal 2 cut-off . In contrast, ESPIRiT with suboptimal 2 cut-off exhibited severe noise (for smaller 2 cut-off ) or residual artifacts (for large

F I G U R E 8
Reconstruction results of the proposed method with different L 1 -norm sparsifying regularization weights ( = 0, 0.001, 0.002, 0.003, and 0.004) at R = 3 and 4. These results demonstrate the spatial nulling maps (SNMs) method as a hybrid-domain method that was flexible and could incorporate L 1 -norm regularization to reduce the noise at the expense of certain degree of image blurring.  Figure S1. The corresponding number of bases is shown in Table S2.
To demonstrate the incorporation of regularization into our proposed method, Figure 8 displays reconstruction results with different L 1 -norm sparsifying regularization weights at R = 3 and 4. Results with = 0.001, 0.002, and 0.003 are shown. With the regularization, the noise was suppressed in the reconstructed images as expected. Figure S2 plots the NRMSE and SSIM values that corresponded to the results shown in Figures 4-6. They again indicated that optimal ESPIRiT reconstruction required the optimization of the eigenvalue threshold values for masking, whereas our proposed method involved no such optimization ( Figure S2A). In general, the proposed method yielded slightly better performance than ESPIRiT in terms of NRMSE and SSIM values within the object regions, especially at high acceleration. The reconstruction improvement of our proposed method over ESPIRiT became more apparent for knee results. As shown in Figure S2B, the proposed method achieved overall better reconstruction than ESPIRiT. More importantly, our method was more tolerant of inaccurate 2 cut-off selection for subspace separation as evident from the NRMSE and SSIM values.
Reconstruction from twofold undersampled data with a FOV smaller than the head size is shown in Figure S3. SENSE using a single set of smooth coil sensitivity maps could not recover the correct coil sensitivity information, leading to severe artifacts. In contrast, GRAPPA, ESPIRiT (with 2 sets of ESPIRiT maps), and SNMs method were able to reconstruct the data. SNMs method provided quality results comparable to ESPIRiT in terms of the noise and artifact levels as indicated by NRMSE and SSIM values.
In addition to the reconstruction using 24 central consecutive k-space lines (Figure 3), the proposed method was further evaluated by reducing the number of central consecutive k-space lines to 16 ( Figure S4). Again, the proposed SNMs method exhibited its robustness and provided reconstruction results that were highly comparable to ESPIRiT at low acceleration factors (R = 2 and 3) and slightly better than ESPIRiT at high acceleration factor (R = 4).

DISCUSSION
This study presents a novel hybrid-domain reconstruction method. It extracts null-subspace bases of calibration matrix constructed from coil calibration data (i.e., autocalibrating or additional calibration scan data) and calculates image-domain spatial nulling maps (SNMs) that contain both coil sensitivity and finite image support information. This null-subspace based method completely eliminates the need for any masking-related procedure. It is also relatively insensitive to subspace cut-off selection, offering a robust reconstruction in practice.

Coil sensitivity and finite image support information in SNMs
The ESPIRiT extracts signal-subspace bases and uses eigendecomposition for calculating ESPIRiT maps, whereas the proposed method extracts null-subspace bases and transforms them into SNMs. They are theoretically different approaches in terms of both the subspace extraction and the method of calculating maps, but share a similar idea of estimating information from subspace bases. The subspace bases calculated through SVD of the constructed block-wise Hankel matrix theoretically can represent all the underlying low-frequency modulations (k-space convolutions with a limited kernel size or image-domain multiplications by smooth maps) in MR data. Coil sensitivity and finite image support are two dominant low-frequency modulations within the MR data. 43,44 The coil sensitivity information is successfully estimated in ESPIRiT and enables signal-subspace reconstruction. In our proposed method, besides estimating coil sensitivity information, we mathematically incorporate the finite image support information by extracting null-subspace bases and combining their image-domain transformations as in Equation (12). This finite image support utilization is based on our analysis that the multiplication between the complement of the finite image support and MR images formulates an image-domain nulling equation as described by Equation (17), indicating the existence of the finite image support information within the null-subspace. The finite support information embedded in SNMs is also evidence from the experimental results shown in Figure 2, which illustrates that the combined SNM delineated the complement of the finite image support of the object to be imaged (i.e., 1 − p). Therefore, our null-subspace based approach and resulting SNMs enable a masking-free hybrid-domain reconstruction method as demonstrated in  In our experiments, the SNMs method directly solved multi-channel images and enhanced joint-sparsity, 18,[40][41][42] whereas ESPIRiT reconstructed and enhanced sparsity on channel-combined image. 21 When the imaging FOV is smaller than the object, the signal cannot be represented correctly using a single set of sensitivity maps within the restricted FOV. Therefore, ESPIRiT uses two sets of ESPIRiT maps to represent the fold-over images and their corresponding sensitivity information. Transformed from multiple null-subspace kernels, the spatial nulling maps are capable of representing such fold-over sensitivity information and can be used to reconstruct the small FOV data. When this occurs, the proposed method achieved quality reconstruction comparable to ESPIRiT ( Figure S3).

Subspace separation
The signal-or null-subspace separation is required in both ESPIRiT and our proposed method. The subspace bases represent the low-frequency modulations within MR data. ESPIRiT obtains coil sensitivity information from signal-subspace bases while our method obtains both coil sensitivity and finite image support information from null-subspace bases. In ESPIRiT, a precise subspace separation is required to ensure the extracted bases satisfy the signal-subspace constraint in Equation (7). If insufficient or excessive bases are extracted, the constraint formed by those bases will be inaccurate. This leads to an inaccurate estimation of coil sensitivity information, causing noise or residual artifacts in the reconstructed results ( Figure 6). In our proposed method, null-subspace bases satisfying the null-subspace constraint in Equation (8) need to be extracted. The excessively extracted bases beyond the null-subspace (i.e., including certain signal-subspace data) cannot satisfy the null-subspace constraint, which may cause incorrect s null ij . Similarly, insufficient null-subspace bases can generate insufficient s null ij . Both scenarios can lead to the inaccuracy of SNMs. However, our proposed method is less sensitive to such subspace separation compared with ESPIRiT, as demonstrated in Figure 6 and Figure S1. Such tolerance to inaccurate subspace separation primarily arises from the fact that the number of signal-subspace bases corresponding to dominant singular values is much smaller than the null-subspace bases because of the high redundancy and high linear dependency within MR data. Another potential factor contributing to such relative insensitivity is that, compared with our proposed method that directly combines multiple s null ij to obtain SNMs, ESPIRiT requires one more eigendecomposition step and introduces manual thresholding before the reconstruction. This procedural difference may cause stronger error propagation in the reconstruction when the subspace separation is inaccurate.

Non-Cartesian imaging
The proposed SNMs method can be extended from Cartesian to non-Cartesian imaging. For non-Cartesian imaging, 45 nonuniform fast Fourier transform (NUFFT) 46,47 can be used to establish the connection between the acquired non-Cartesian data and estimated Cartesian data. Specifically, images can be calculated by inverse NUFFT and sequentially transformed into Cartesian estimations. Null-subspace kernels can then be constructed from the Cartesian estimations and transformed into spatial nulling maps. The images can be reconstructed iteratively as illustrated in Figure S5. In each iteration, the spatial nulling maps are updated using the reconstructed images from the previous iteration, and used for reconstructing updated images. Data consistency is enforced both during the optimization procedure of Equations (15) or (16) and between every iteration. Results of reconstructing spiral data using the SNMs method are shown in Figure S6.

Computational time
In our implementation, we used a personal desktop computer with 4-core i5-6500 CPU. The reconstruction codes were written using MATLAB R2020b. Our proposed method is computationally efficient. In this study, the calculation of 6-channel SNMs for a single 2D 256 × 218 slice at R = 4 and the reconstruction took ∼2 and ∼6 s without and with L 1 -norm sparsifying regularization, respectively. The computational times for the 2D single-slice reconstruction using the proposed method and ESPIRiT from the first subject in Figure 3 are summarized in Table 1.
The computational efficiency of our proposed method was similar to ESPIRiT, yet much higher than PRUNO that requires many iterations to converge and each iteration can take ∼1 to 3 s for 2D single-slice reconstruction (e.g., it took ∼3 mins for the 6-channel 2D single-slice reconstruction at R = 4.)

Comparison to k-space PRUNO
As the original null-subspace method, PRUNO solves the inverse problem in the k-space and the proposed method converts the problem and solves it in the image domain. They share the same idea of using null-subspace and produce similar quality results ( Figure S7). However, the proposed method offers advantages in terms of computational efficiency and explicit usage of maps containing coil sensitivity information. First, our hybrid-domain SNMs method is computationally more efficient than k-space PRUNO.

T A B L E 1
Computational time corresponding to the single-slice 2D reconstruction from the first subject shown in Figure 3.
No regularization When solving the inverse problem by least square minimization, the image domain minimization in the proposed method can exhibit different converging behavior (e.g., converging speed and numerical stability 48 during conjugate gradient descent) compared to the k-space PRUNO. This is evident from that PRUNO reconstructs data very slowly and needs the extra step to perform GRAPPA as an initial guess for increasing convergence speed and/or more accurate estimates. 16 Even so, PRUNO still entails a long reconstruction time. Compared to PRUNO, the proposed SNMs method is much faster and can reconstruct one slice within a few seconds (Table 1). More importantly, the SNMs method does not require any initial guess, and still converges fast and robustly. Therefore, the proposed method is efficient and this feature is important in practice. Furthermore, with the fast convergency, the proposed method can be extended to potentially real-time 3D or dynamic imaging in the future. Second, the proposed method derives image-domain spatial nulling maps with explicit coil sensitivity information. Such image-domain coil maps offer more choices in using prior information related to coil sensitivity. For example, in recent studies, 49,50 the coil sensitivity information can be directly estimated by a deep learning neural network from the undersampled imaging data, where prior knowledge such as subject-coil geometry information is used. They offer parallel imaging approaches that require a few number of or no ACS lines. Such deep learning strategy can also be combined with the proposed SNMs method, again illustrating the flexibility of the proposed hybrid-domain method when compared to k-space PRUNO.

CONCLUSION
We provide a flexible and efficient hybrid-domain parallel imaging reconstruction method that extracts null-subspace bases of calibration matrix to calculate image-domain SNMs. Multi-channel images are reconstructed by solving a nulling system formed by SNMs that explicitly contain coil sensitivity and finite image support information. The proposed hybrid-domain method shows quality reconstruction that is highly comparable to ESPIRiT with optimal manual masking. More importantly, it eliminates the need for coil sensitivity masking and is relatively insensitive to the subspace division, therefore, offering a robust parallel imaging reconstruction procedure in practice.

SUPPORTING INFORMATION
Additional supporting information may be found in the online version of the article at the publisher's website. Figure S1. Reconstruction results of the proposed method versus ESPIRiT results (using manually optimized eigenvalue threshold) with extremely small or large 2 cut-off for subspace division at R = 3 for two 6-channel brain datasets. The datasets were the same as those in Figure 7. Even with 2 cut-off that was far from the optimal selection, the proposed method achieved fewer artifacts or less noise than ESPIRiT, demonstrating that the proposed method was less sensitive to subspace cut-off selection than ESPIRiT. Figure S2. (A) NRMSE and SSIM values corresponding to the brain image reconstruction results from two subjects shown in Figure 4 and knee image results from Figure 5 with different eigenvalue thresholds. (B) NRMSE and SSIM values with different 2 cut-off for subspace division, corresponding to the brain image reconstruction results from two subjects shown in Figure 6. Figure S3. Reconstruction from twofold undersampled data with a FOV smaller than the head size. The dataset was from ESPIRiT open sources 21 and acquired using a 2D spin-echo sequence (TR/TE = 550/14 ms, FA = 90 • , BW = 19 kHz, matrix size = 320 × 256, slices = 6, slice thickness = 3 mm). 24 central consecutive k-space lines were preserved out of 256 lines. SENSE could not reconstruct the image correctly. The SENSE reconstructed image suffered from severe artifacts. GRAPPA, ESPIRiT (with 2 sets of ESPIRiT maps), and SNMs method were able to reconstruct the image. SNMs method provided quality results comparable to ESPIRiT using two sets of maps, in terms of the noise and artifact levels as indicated by NRMSE and SSIM values. Figure S4. Reconstruction results of the 6-channel T1W GRE brain data by preserving only 16 central consecutive k-space lines. Figure S5. Flowchart of reconstructing non-Cartesian data. Nonuniform fast Fourier transform (NUFFT) establishes the connection between the acquired non-Cartesian data and estimated Cartesian data. The images can be reconstructed by iteratively solving the nulling system equation, during which the spatial nulling maps are updated and data consistency is enforced. Figure S6. Reconstruction results of 2D 8-channel T2W SE spiral data. Fully sampled data were acquired according to an 8-shot regular spiral trajectory, with acquisition window = 21 ms, TR/TE = 2700/54 ms, FOV = 220 × 220 mm 2 , slice number = 9, and spectral presaturation with inversion recovery used for fat suppression. The spiral interleaves spread evenly along the azimuthal direction and had the same rotation angle among adjacent slices. Undersampling was performed by discarding the spiral shots in an interleaved way. Images were reconstructed from 2 shots using the SNMs method. The non-Cartesian reconstruction further demonstrates the flexibility of the SNMs method. Figure S7. Reconstruction results of the 6-channel T1W GRE brain data using SNMs method, ESPIRiT with manually optimized parameters, and PRUNO. A total of 24 central consecutive k-space lines were preserved. SNMs method and PRUNO produced similar reconstruction results in terms of noise and artifacts. Table S1. The number of signal-and null-subspace bases corresponding to different 2 cut-off . Table S2. The number of signal-and null-subspace bases for a wider range of 2 cut-off corresponding to Figure S1.