Abstract
Slice-to-volume reconstruction (SVR) method can deal well with motion artifacts and provide high-quality 3D image data for fetal brain MRI. However, the problem of sparse sampling is not well addressed in the SVR method. In this paper, we mainly focus on the sparse volume reconstruction of fetal brain MRI from multiple stacks corrupted with motion artifacts. Based on the SVR framework, our approach includes the slice-to-volume 2D/3D registration, the point spread function- (PSF-) based volume update, and the adaptive kernel regression-based volume update. The adaptive kernel regression can deal well with the sparse sampling data and enhance the detailed preservation by capturing the local structure through covariance matrix. Experimental results performed on clinical data show that kernel regression results in statistical improvement of image quality for sparse sampling data with the parameter setting of the structure sensitivity 0.4, the steering kernel size of and steering smoothing bandwidth of 0.5. The computational performance of the proposed GPU-based method can be over 90 times faster than that on CPU.
1. Introduction
Magnetic resonance imaging (MRI) is an ideal diagnostic technique for researchers to investigate the development of the fetal brain [1]. Its advantages are the absence of ionizing radiation, the availability of different contrast options (T1-weighted, T2-weighted, and diffusion-weighted imaging), and the superior contrast of soft tissue compared with ultrasonography, and MRI is also a safe and noninvasive procedure for patients and fetuses [2–4]. For these reasons, MRI has been widely used to investigate the developing fetal brain in vivo [5]. For fetal brain MRI, the high-quality volume representation of 3D acquisition has significant clinical meaning [6]. By the observation of the reconstructed volume data, researchers can study the mechanism of brain development and maturation [7] and identify the fetal brain abnormality or potential injury [8, 9], such as brain tumors, vascular malformations, and posterior fossa abnormalities. Fetal brain MRI can provide abundant information about aid clinical management, prognostication, and counseling [10].
The duration of an examination is typically 45 to 60 minutes for fetal brain MRI [1]. One major problem of fetal brain MRI is motion artifacts caused by fetal and maternal motion, because of the long acquisition times of 3D MRI scanning. Maternal motion may be avoided by some measures, but fetal motion is usually fast and unpredictable, especially for the younger fetus. Thus, it is still challenging to reconstruct high-fidelity image for fetal brain MRI due to the presence of fetal motion. For fetal motion, different strategies can be adopted to reduce the motion artifacts on MRI [11]. The first strategy tries to prevent the motion occurring during the examination, such as maternal sedation. The second one tries to quicken the data sampling speed. The faster the acquisition techniques for fetal brain MRI are, the lower the motion occurs. For example, the single-shot fast spin echo (SSFSE) T2-weighted imaging can acquire a slice at 1-second speed [12]. On the other hand, sparse data sampling technique can be applied to shorten the time of data acquisition. The last strategy tries to reconstruct high-quality image through advanced postprocessing motion detection and correction algorithms, such as the SVR method [13].
For the SVR framework [14], it includes the following steps to reduce the fetal motion and reconstruct the high-quality 3D result: motion identification and exclusion step, registration step, reconstruction step, and regularization step. For the motion identification and exclusion step, we should estimate the amount of motion and exclude the slices with large amount of motion corruption. Early reconstruction approaches need to manually exclude the motion corrupted slices. The intersection-based motion correction approach can automatically detect and reject motion corrupted and incorrect registration slices by the abnormal level of their mean squared intensity difference with respect to all other intersecting slices [15]. In [16], Kainz et al. have proposed an approach to automatically estimate the amount of motion based on the low-rank decomposition for linearly correlated image slices [17]. Using this approach, we can reject stacks with large motion and choose the stack with the least motion as the template to prepare for the registration step. Registration step can be utilized to correct the motion between slice and the reconstructed volume. Rousseau et al. [18] combined the 2D/3D registration with the PSF to achieve the 3D reconstruction. PSF [14] is a mathematical function to model the actual appearance of data points in physical space. By PSF, we can physically correct estimation of the image acquisition process. Subsequently, the SVR method was modified to improve the robustness of the 2D/3D registration [19]. For the reconstruction step, superresolution methods [20, 21] are utilized to reconstruct the 3D volume. In [22], Gholipour and Warfield combined the superresolution method with slice-to-volume registration to reduce the burring effect. Because the motion identification and exclusion steps can exclude the slice of which the motion amount is greater than the threshold, the amounts of the slightly corrupted slices are still preserved for reconstruction. Using the robust superresolution volume reconstruction method [23], the weight of slightly corrupted and misaligned slices would be reduced to minimize the effect of motion. During the process of superresolution reconstruction, maximum likelihood estimation (MLE) is treated as an optimum solution to estimate the point’s value [24]. To get better results, we should minimize the difference between the estimated slices and the acquired slices. Since the minimization only depends on the acquired samples, the estimation in the MLE framework is ill-posed and inaccurate when the samples are sparse [23]. The regularization step is used to solve the overfitting problem, and it can reduce image noise and registration errors. In [25], Charbonnier et al. proposed a deterministic edge-preserving regularization method to deal with image. However, this method makes it difficult to avoid the smoothing of edges. Adaptive regularization techniques can be employed to reduce the smoothing effects of regularization [26]. In [27], Rousseau et al. took advantage of total variation regulation to extend the superresolution reconstruction method.
The general SVR framework with the superresolution reconstruction method has been developed in [28]. One important way to alleviate fetal motion is to quicken the data acquisition time by the sparse data sampling technique. However, the traditional SVR method could not deal well with the sparse sampling problem and cannot provide high-quality image. In this paper, we utilize the SVR method with adaptive kernel regression to cope with the sparse volume reconstruction with minimum motion artifacts under the condition of sparse data acquisition. The key improvements compared to previous works are as follows: firstly, we make use of the sparse samples to get faster speed of data acquisition in fetal brain MRI. Next, the adaptive kernel regression-based reconstruction method [29] with robust statistics calculation [24] can reconstruct high-quality volume under the condition of sparse sampling. In general, our comprehensive reconstruction method for fetal brain MRI mainly includes slice-to-volume registration, the robust statistics calculation, the PSF-based volume update, and adaptive kernel regression-based volume update.
The rest of the paper is organized as follows. The detailed methodology is discussed in Section 2. We design the actual implementation of the algorithm in Section 3. Section 4 involves the experiment results and compares with those of superresolution methods. In this section, we also discuss how to determine the optimal values of related parameters using GPU-based fast reconstruction. Finally, we make a brief conclusion in Section 5.
2. Methods
2.1. Model of Data Acquisition and Motion Estimation
During data acquisition of fetal brain MRI, we collected several stacks of 2D slices in different orientations. Because of the fetal motion, the movement could be observed between these slices. Assume that the acquired misaligned 2D slices are , , and the corresponding sparse 2D slices are , . During the slice acquisitions of MRI, the inhomogeneity of the magnetic field , , affects the intensities of the slices and the scaling factor , , is potentially different for each acquired slices. In [30], the logarithmic transformation was chosen to make the bias additive. However, field in-homogeneities are known to be multiplicative. Differently, we use the multiplicative bias field to form the multiplicative exponential model which replaces the logarithmic model. So the scaled and bias corrected slice can be modeled as where is the sparsely sampled slice coming from the sparse operator , is the vectorization operator that transforms a -pixel () image into a vector of intensity values . The corresponding -aligned 2D ground-truth slices are , . The relationship between corrected slices and the ground-truth slices can be denoted as follows: where is the motion error, and denotes the unknown motion transformation parameter of slice . Then, we can define the following data matrix: where , , , and denote the observed data matrix, reconstructed data matrix, motion error matrix, and the rigid transformation matrix. Given these definitions, the observed data matrix can be described as . The motion error matrix is mainly caused by misaligned slices. The misaligned slices can cause the inaccurate reconstructed volume, and we want to exclude the stack which has many misaligned slices. However, we cannot directly calculate the amount of stack motions for the observed data matrix , but a low-rank approximation as surrogate estimate can be used to evaluate the stack motion indirectly [16]. It has been shown that provides the best approximation to [31]. The difference value between and measures the motion error . The smaller difference value indicates that the stack has fewer motions. To provide the low-rank approximation, the singular value decomposition is used to decompose the data matrix as . The singular value decomposition of produces three matrices , , and . and are both orthogonal matrices, and is the diagonal matrix containing the singular values on the diagonal. And the singular value decomposition of is the first singular values of the original matrix , i.e., , . and are the first columns of and , and is the top left submatrix of . The relative error based on the Frobenius norm is used to measure the approximation between and , i.e. . For the different values of , we can find the minimal rank for each stack that satisfies the given threshold , i.e., . Combining and , the surrogate estimate for the amount of motion is given by .
Based on the low-rank decomposition method, we can choose one stack with minimal motion as the target template and first perform the 3D rigid volumetric registration between the target template and the other stacks (stack to template registration). During the first registration, we can get the corresponding rigid global transformation matrix . Then, second, the 3D rigid volumetric registration between the reconstructed volume and all slices (slice to reconstructed volume registration) can produce local transformation matrix . The prerequisite for two registrations is that all stacks and reconstructed volume should be mapped to the world coordinates. Thus, we need to define two transformations to map each pixel in the 2D slice and each voxel in the reconstructed volume to a continuous location in the world coordinates. The first one is world transformation that transforms the discrete coordinates of a pixel in the acquired slice to the continuous local world coordinates. The second one is world transformation that transforms the discrete coordinates of a voxel in the reconstructed volume to the continuous local world coordinates. Meanwhile, the mapping and registrations can be combined and formulated as Equation (4). Thus, Figure 1 illustrates the whole transformation process from the pixels in the sparse slice to voxels in the 3D reconstructed volume.
2.2. PSF-Based Volume Update
To model the actual appearance of sampling data points in physical space, the point spread functions (PSFs) are used to make the exact estimation for every voxel value in the reconstructed target volume. For the MRI ssFSE sequence in this paper, the exact shape of the PSF has been measured using a phantom and rotating imaging encoding gradient in [14]. The resulting shapes of the PSF in in-plane and in through-slice are given by a sinc function and the slice profile, respectively. Since the ideal rectangle profile has the very dense and inefficient spatial sampling, Kuklisova-Murgasova et al. [28] have proposed to use the 3D Gaussian function with the full width at half maximum (FWHM) equal to the slice thickness as an approximation for the sinc function. The PSF function based on 3D Gaussian profile is used to approximately model the SSFSE sequence and is expressed as follows: where , , and are the offsets from the center of a reconstructed voxel, and are the full width at half maximum (FWHM) in the in-plane - and -directions, and the equals to the slice thickness in the through-plane direction. For each pixel in the sampled slice, the is applied to obtain the corresponding PSF coefficient matrix. Since every sampling pixel (i.e., ) does not perfectly align itself with the reconstructed voxel (i.e., ), one contributes to more than one . To model this, every voxel is sampled around its local surrounding neighbor in the reconstructed volume to make sure that it has at least one corresponding pixel in the acquired slices. Then, the PSF coefficients are used to weigh the pixel’s contribution during theth iteration. where is the operation that finds the nearest voxel center in the space of the reconstructed volume. The reconstructed volume is updated iteratively through the PSF-based data sampling model, and every voxel of is filled at an arbitrarily chosen voxel size.
2.3. Robust Outlier Removal
Once the target volume is updated based on the Gaussian PSF, the simulated slices can be generated from the updated reconstructed volume. Then, the misaligned error between the corrected acquired sparse slices and simulated slices can be computed as
In [28], an EM model-based robust statistics approach was proposed to classify each slice pixel into two classes: inliers and outliers. Specially, the probability density function (PDF) for the inlier class is modeled as a zero-mean Gaussian distribution with variance : , and the PDF for the outlier class is modeled as a uniform distribution with constant density, which is a reciprocal of the range :. Then, the likelihood of the observing error can be expressed as where is a mixing proportion of inliers representing the correctly matched voxels. Then, the posterior probability of a voxel being an inlier can be computed using the expectation step as
The variables and are updated by the following maximization step: where is the number of the pixels in the slice. By constantly iterating, we can get the best parameters and c. The inlier probability can be used to weigh the PSF-based volume update. By the same way, each slice is classified into inlier and outlier as well using the EM algorithm. The probability of an inlier slice is defined as . The slices inferred to be an outlier are excluded from the PSF-based volume update to remove artifacts of motion corruption and misregistration.
The purpose of the outlier removal is to make the framework more robust by rejecting the outlier slices. The outlier removal module is adopted directly from the cited previous work [16], where the accuracy of the motion recognition and outlier removal has been evaluated in detail by simulating the slice motion at a variety of amplitudes and comparing the known motion amplitude to the surrogate measure provided through rank approximation. They have shown that there was strong correlation between the amplitude of the known motion and the values of derived from the stack data matrices.
2.4. Steering Kernel Regression-Based Volume Update
For sparse reconstruction, it is experimentally found that the reconstructed volume still remains unallocated or inaccurate voxels after PSF-based volume update and the reconstructed result is noise as shown in Figure 2.
In [32], the kernel regression can make better nonparametric estimation for the empty pixels. In this paper, the steering kernel regression approach [29] is introduced to update the voxels for the previous sparse volume data. The model for the kernel regression is expressed as where is the function of kernel regression, is the 3D coordinate of the voxel, is a zero-mean Gaussian noise with variance as , and is the voxel after PSF-based Gaussian volume update.
Assuming that the voxel is close to the known voxel in the reconstructed volume, we have the following approximation for using the -term-order Taylor series: where and are, respectively, the gradient () and Hessian operators; , which is the voxel value of interest; and the vectors and are defined as is the half-vectorization operator that transforms the upper triangular portion of a symmetric matrix into a column-stacked vector, i.e.,
Based on the least-squares formula, we can optimize Equation (12) as where is the number of known voxels within the neighborhood window, is the distance-weighted kernel function which penalizes distance away from the local position, and is the smoothing parameter that controls the strength of the penalty. The kernel function is chosen as the exponential function, Gaussian function, or other functions which satisfy the following conditions:
For the computation simplicity, the Gaussian-based kernel function is chosen in the steering kernel regression [33]. The steering kernel adapts locally to image structures (e.g., edges, flat, and texture areas), which are captured by the kernel footprint. For example, the kernel footprint is large in the flat areas, elongated in edge areas, and compact in texture areas. The 3D steering kernel function takes from where is the norm and is the symmetric covariance matrix. Since the local image structure is highly related to the gradient covariance, we can make the data-dependent covariance matrix estimation utilizing the local edge gradients: where is a local analysis window and , , and are the gradients along the -, -, and -directions.
Equation (16) can be expressed in the matrix form as where is the vector set of all known voxels, is the vector set of all estimated parameters, is the diagonal matrix whose elements on the diagonal are the value of , and the other elements are zero. According to the least-squares method, we have the following solution: where , is the voxel value estimated by the steering kernel regression, is applied for computing the symmetric covariance matrix iteratively, and is a coordinate matrix expressed as follows:
Once the reconstructed volume is updated based on the steering kernel regression, we update the simulated slices and the misaligned error according to Equation (7). To remove artifacts caused by motion corruption and misregistration and enhance image edges, we further update the reconstructed volume using the following equation:
3. Implementation
The experiment computer is equipped with Intel Core i5 2.6 GHz CPU, and the operating system is Windows 7 64 bit. We have implemented the proposed algorithm using the Microsoft Visual Studio 2012 and Image Registration Toolkit (IRTK) software package which includes many useful methods to do registration, transformation, and other image processing. In this section, we discuss the key implementation details. The diagram of the total algorithm is expressed in Figure 3.
The first step is to evaluate the stack motion according to the method of low-rank decomposition. We estimate the amount of the stack motion by the surrogate and choose the stack with the minimum amount of stack motion as the template. The second step is to perform the global registration, which calculates the matrix of global transformation from the other stacks to the template. The third step is the iterative registration-based volume reconstruction, which consists of the outer registration step and the inner reconstruction step. The outer loop step includes the PSF-based volume update, robust outlier removal, steering kernel regression-based volume update, and slice to volume registration. The PSF-based volume update step makes the initial estimation of the reconstructed volume based on Equation (6). Then, the simulated slices are created and used for the robust misaligned error calculation between the simulated slices and the acquired slices as described in Section 2.3. The robust statistic calculation achieves the classification of outlier slices and inlier slices. The outlier slices are excluded to remove artifacts of motion corruption and misregistration. The slice to volume registration is to calculate the local transformation from slices to reconstructed volume. The whole transformation process is described by Equation (4). The volume update based on the adaptive steering kernel regression is aimed at reconstructing the accurate volume iteratively as shown in Figure 4. The initial gradients are estimated by the classical kernel regression. Then, the gradient information is used to calculate the covariance smoothing matrix (i.e., Equation (19)). We use smoothing matrix to update the voxel value and its corresponding gradients according to Equation (21), respectively. To obtain a more reliable voxel estimation, the process is iterated three times in our experiment.
4. Experimental Results and Evaluation
4.1. Evaluation of Image Quality
In the experimental evaluation, we used the datasets from the fetal MRI datasets [16], which were acquired by a Philips Achieva 3 T MR scanner. During the experiment, the volunteers were lying at a tilt on the left side to avoid the pressure on the inferior vena cava. The volunteer’s womb was scanned with single-shot fast spin echo (SSFSE) T2-weighted sequence. Three stacks of images from axial, coronal, and sagittal orientation are used to construct the final high-resolution volume. To obtain the sparse stacks, we randomly remove different proportions of pixels of the stack once every 10% proportion ranging from 10% to 90%. The different removal proportions control the removal number of pixels. The typical 30th slice of the collected stack and its corresponding simulated spare slices are illustrated in Figure 5.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
For different data removal ratios, the sparse stacks are used to reconstruct the high-resolution 3D fetal brain MRI volume with the method of Kainz et al. [16] (SVR with superresolution) and our proposed method. Figure 6 shows the reconstructed results by Kainz et al.’s method and the proposed method for the sparsely sampled dataset with once every 10% data removal ratio ranging from 0% to 90%, respectively. In Figure 6, we can observe that as the removal ratio increases, the reconstructed results by Kainz et al. method have much more noise for the sparse sampled dataset compared with our proposed method. On the other hand, the proposed method is capable of reconstructing high-resolution images without obvious artifacts even for the 90% data removal ratio.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(q)
(r)
(s)
(t)
For the sake of quantitative evaluation, the image quality assessment index of root mean square error (RMSE) [9] and mean structure similarity (MSSIM) [34] is introduced to quantitatively assess the algorithms under different removal ratios. The RMSE score can be computed by the following equation: where is the reconstructed result, is the ground-truth volume, and is the number of voxels. A good reconstruction method is capable of estimating the removal data very close to the original data. Given and, a low RMSE value represents that the estimated result is satisfying while a high RMSE means that the interpolation accuracy is poor.
The structure similarity (SSIM) index explores the structural information for image quality assessment based on the main idea that the pixels have strong interdependency when they are spatially close. The SSIM metric is calculated based on the intensity, contrast, and structure and is computed as where , , , , and denote the mean, variance, and covariance on square window, which moves pixel by pixel in images and , respectively. The two variables and are used to stabilize the division with weak denominator. Here, is the dynamic range of pixel value (e.g., 255 for 8-bit grayscale image), with and by default. Since the SSIM metric is calculated on various windows of a volume image, the mean SSIM (MSSIM) index is used in this experiment to assess the overall image quality: where is the number of local windows in the image. ; the higher MSSIM indicates better structural similarity between two images.
For the clinical datasets, it is impractical to obtain the ground-truth volume in advance. For the sake of fair comparison among different methods, the quantitative evaluation is performed based on an average reconstructed volume. We first use the original stacks without data removal to reconstruct a complete volume by Kainz et al.’s method (2015) and our method (e.g., Figures 6(a) and 6(b)), respectively. Both volumes are adopted to create an average volume as the ground truth. Table 1 shows the quantitative results of the RMSE and MSSIM values with different data removal ratios for each method. As can be seen, the results of Kainz et al.’s method produce the highest RMSE scores and lowest scores for all sampling rates. Both the high RMSE value and low MSSIM value for Kainz et al.’s method indicate poor image quality because of the artifacts and noise. For all levels of sampling rate, the proposed method performs better than the Kainz et al.’s method. More importantly, both of the difference of the RMSE and MSSIM index between Kainz et al.’s method and our method increase while the data removal ratio increases, indicating that our method outperforms much more compared with the Kainz et al.’s method when the data removal ratio increases.
4.2. Evaluation of Computational Efficiency
Our approach is capable of reconstructing the accurate volume from the highly sparse sampling dataset, but it requires largely computational burden as well due to the iterative kernel regression estimation. To reduce the long processing time of the adaptive kernel regression, the proposed method is accelerated by the GPU-based parallel implementation based on the NVIDIA GeForce GTX 1080 and CUDA 8.0 libraries. In the experiment, we make the evaluation of the computational efficiency of the adaptive kernel regression method, including the computation of the gradient information, the covariance smoothing matrix, and the steering kernel regression. The computational efficiency of the other modules (i.e., motion estimation, stack-to-template registration, PSF-based volume update, robust outlier removal, and slice-to-volume registration) has been evaluated in detail in [16]. The comparisons are based on the single-threaded CPU, multithreaded CPU, and GPU for the dataset of 80% data removal ratio under the parameter setting as the kernel size and the smoothing parameter in the initial gradient estimation step based on the classical kernel regression, the steering kernel size and the steering smoothing parameter in the steering kernel regression step, and the window size , the regularization parameter , and the structure sensitivity . The practical running time for the proposed method is shown in Table 2. For single-threaded CPU, the running time of the adaptive kernel regression method is 1865.769 s, which includes the computation of the gradient information in 416.96 s, the covariance smoothing matrix in 46.082 s, and the steering kernel regression in 1402.727 s. For the multithreaded CPU, we use 4 threads to run the adaptive kernel regression method and its computational time is 1045.293 s, indicating less improvement compared with the single-threaded CPUs. The running time of the GPU implementation is 20.718 s in total. From Table 2, it can be observed that the GPU-based processing time has significantly decreased by 98.89% and 98.02%, compared with the single-threaded CPU and the multithreaded CPU, respectively.
4.3. The Choice of the Adaptive Kernel Regression Parameters
There are seven parameters which can be adjusted to affect the reconstructed image quality for the proposed method. These parameters include the kernel size and the smoothing parameter (i.e., the kernel bandwidth) in the initial gradient estimation step based on the classical kernel regression, the steering kernel size and the steering smoothing parameter in the steering kernel regression step, and the window size , the regularization parameter , and the structure sensitivity () in the covariance matrix estimation step. In our method, and are related with the initial calculation of gradient information and have a negligible effect in the experiment. For the adaptive sparse reconstruction, covariance matrix estimation and steering kernel estimation are the two of the important steps and their parameters (i.e., , , , , and ) play an important role in the volume reconstruction and deserve much more investigation.
With the help of GPU-based fast implementation, we firstly adjust the parameters (i.e., , , and ) of the covariance matrix estimation one by one. The window size decides how many neighbor points in the gradient matrix are taken for the estimation of the covariance matrix. Table 3 shows the RMSE and MSSIM values and the running time for different window sizes. Both of the RMSE and MSSIM values differ slightly, indicating that the window size has a negligible influence on the reconstructed image quality, as shown in Figure 7. However, the running time increases with the increase of window size. It can be observed that the window size of is chosen because of its faster implementation and lower RMSE value.
(a)
(b)
(c)
(d)
Table 4 shows the RMSE and MSSIM values influenced by the structure sensitivity parameter , and the lowest RMSE value and highest MSSIM value are obtained for the structure sensitivity indicating the best performance of the algorithm. Figure 8 shows the corresponding reconstructed images for different values. As can be seen, the result with large structure sensitivity (e.g.,) results in oversmoothing image, while small structure sensitivity (e.g.,) overemphasizes the image edges. The experiment shows that the structure sensitivity has a significant influence on the reconstructed volume.
(a)
(b)
(c)
(d)
Under different regularization parameter settings, the RMSE and MSSIM measurements of the reconstructed results are calculated and shown in Table 5. The illustrative results are further shown in Figure 9. The regularization parameter is used to suppress the noise. However, the regularization parameter has negligible influence on the reconstructed image quality in the experiments.
(a)
(b)
(c)
(d)
The next group parameters (i.e., and ) come from the steering kernel regression for the adaptive voxel value estimation. The kernel window size has a great impact on the processing time for the kernel regression-based algorithm under different data removal proportions. When the kernel window increases, the estimation of each voxel involves more nearby pixels and leads to more computation [32]. The smaller the kernel window size is, the faster our algorithm runs. On the other hand, if the size of the kernel window is too small, we could obtain the fault result, because there are not enough samples to make the current voxel estimation, especially for large data removal proportion. The larger the data removal proportion is, the sparser the sampled data will be. The RMSE and MSSIM index and processing time measurement of the reconstructed results under different kernel window sizes are shown in Table 6. With the increase of the kernel window size, the running time of steering kernel regression function is becoming longer. The corresponding images of different steering kernel sizes are shown in Figure 10. The proper kernel window size (i.e., ) produces a trade-off between the processing time and the reconstruction accuracy under different removal proportions.
(a)
(b)
(c)
(d)
Table 7 shows the RMSE and MSSIM values and running time with different steering smoothing parameters . As can be seen, the results with the steering smoothing parameter (i.e., ) achieve the lowest RMSE value and highest MSSIM value among these settings. The reconstructed results produced by different steering smoothing parameters are shown in Figure 11. In [33], it has been given that the steering smoothing parameter indicates the “footprint” of the kernel function. The large footprint of the kernel function could reduce the noise but at the cost of oversmoothing details, while small footprints are desirable to preserve the edges. In the experiment, the footprint setting is chosen for reaching a trade-off between the noise reduction and edge preservation. Finally, all parameters of the adaptive steering kernel regression algorithm are determined as follows: the window size , the regularization parameter , the structure sensitivity , the steering kernel size , and the steering smoothing parameter . Under such parameter setting, the RMSE value decreases from 126.47 to 78.89, indicating the quality improvement by 37.62%.
(a)
(b)
(c)
(d)
5. Conclusion
In this paper, we proposed an adaptive reconstruction method to deal with the sparse sampling dataset for fetal brain MRI. Our method combines the latest SVR framework, including the stack motion evaluation, PSF-based volume update, robust outlier removal, slice-to-volume registration, and the proposed adaptive kernel regression-based volume update. Compared with the existing SVR framework, our method has advantages of sparse volume reconstruction and is capable of reconstructing superresolution image even for 80%~90% data removal. With the capability of sparse reconstruction, the data sampling time can be greatly shortened and thus, the motion artifacts can be reduced indirectly. To accelerate the voxel estimation, we use the CUDA to implement the steering kernel regression approach. For the proposed method, the running times of GPU-based implementation are speeded up to 90x than that of the CPU. The GPU-based parallel implementation of the proposed method is more practical to meet the requirements of fetal brain MRI. Meanwhile, we make the detailed investigation on the choice of parameters for the adaptive kernel regression-based volume reconstruction with the help of GPU-based fast implementation. To summarize, the structure sensitivity and the steering kernel window size are two of the important parameters on sparse kernel regression volume reconstruction. Meanwhile, the kernel window size has a strong relationship with the running time. Larger window size requires longer processing time. Overall, our approach is used to reconstruct superresolution image from the highly sparse sampled dataset of fetal brain MRI corrupted with motion artifacts. One of its potential applications includes other motion organ MRI reconstruction, such as the heart MRI with the heart beating motion artifacts.
Data Availability
The test data was downloaded from the publicly available dataset on GitHub (https://github.com/bkainz/fetalReconstruction.git).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Qian Ni and Yi Zhang contributed equally to this work and should be considered co-first authors.
Acknowledgments
This study was financed partially by the National Key R&D Program of China (No. 2018YFA0704102), the National Natural Science Foundation of China (Nos. 81827805, 61401451), and the Shenzhen Key R&D Program (No. JCYJ20200109114812361).