1 Introduction
Functional magnetic resonance imaging (fMRI) offers a unique means of observing the functional brain architecture and its variation during development, aging, or disease. Since its invention in the early 1990s, fMRI has revolutionized neuroscience by characterizing functional connectivity (FC) as the correlation between regional brain activity over time. Despite the insights into network formation and functional growth of the brain, in utero fMRI of living human fetuses, however, remains challenging. Unconstrained and potentially large movements of the fetuses during fMRI acquisition result in artifacts such as inplane image blurring, slice crosstalk and spinhistory artifacts that can considerably affect the image quality and bias any subsequent conclusions on the FC of the developing brain.
Some of the early works in this area established tools from preprocessing of adult fMRI data (eg. MCFLIRT [1]) and censored parts of the data with excessive motion (scrubbing). The reported rejection rate in these studies is very high and many subjects were removed since there was not enough data left after scrubbing. Further attempts to correct for fetal fMRI motionrelated artifacts mainly relied on improving the estimate of the realignment parameters for potentially large movement of the fetuses. Ferrazzi et al. [2] proposed a comprehensive processing framework that incorporates bias field and spinhistory correction, combined with slice to volume registration and scattered data interpolation, which allowed the data to be analyzed within a single consistent space . Furthermore, to avoid that surrounding maternal tissues dominate the similarity metric of registration algorithms, [3]
trained a convolutional neural network for rapid and accurate isolation of the fetal brain and then utilized rigid registration for motion correction. You et al. also in
[4] present an optimized design based on the MCFLIRT registration to identify independently moving objects (i.e., the placenta and fetal brain) and estimate the local motion for each. All these studies, however, used 3D volume interpolation of each time frame independently to reconstruct the whole FC data. More recently, [5]extended super resolution technique (will be discussed in the next section) with separate reconstruction for each 3D volume within the time sequence. Since inutero motion is unconstrained and complex, the rate of scattered points versus the points to infer is less than 1 to 1 in each 3D volume. Therefore, taking the advantage of temporal structure in fMRI imaging can potentially improve the reconstruction of image in this challenging problem. Here we propose a new method that, rather than treating each volume independently, takes both spatial and the temporal domain into account and iteratively reconstruct 4D fetal fMRI time series.
2 Method
Our method directly reconstructs fetal fMRI time series on a 4D regular grid from motion scattered slices. We formulate our problem as a maximum likelihood (ML) estimation in which the conditional probability density function (pdf) of the acquired slices given the current estimate of the reconstruction is maximized:
(1) 
Here, denotes an estimate of the final 4D reconstructed image i.e. , and is the kth acquired slice. is a prior on and is a constant with respect to
. Using the assumption of Gaussian noise with zero mean and standard deviation of
, and independent acquisition of each slice, the above conditional pdf can be written as(2) 
where are the samples from the acquired slices , and are the samples from the estimated slices by considering a slice acquisition model as . is the subsampling matrix, is the blur matrix representing the point spread function (psf), and is the observation noise. The motion matrix
is defined for each slice through slicetovolume registration as a 6DOF 3D rigid transformation (including three rotations and three translations). When the noise residuals are presumed to be drawn from a Gaussian distribution, and using the logarithmic function, the above problem is equivalent to obtaining the Maximum A posteriori Probability (MAP) estimate by minimizing the log likelihood function:
(3) 
The ML formulation is applied in structural fetal T2weighted MRI where a highresolution volume is reconstructed using several orthogonal series of thick slices, called as superresolution technique [6, 7, 8, 9]. The point spread function used in that case is derived directly from the image acquisition process. In fMRI analysis, however, lowfrequency fluctuations of the signal over time is desired and postprocessing commonly includes lowpass filtering to derive temporally smooth signal. In this work we propose, rather than using a separate postprocessing filter, to integrate this step within the reconstruction framework by using a four dimensional point spread function. This also provides further stability of the reconstruction by incorporating information from the neighboring slices in time in case of a large motion and missing data. The 4D psf can be formulated as:
(4) 
where the values represent the point spread in x, y, z, and t, which are the neighborhood of a voxel is defined on both spatial and temporal domain. To create a distance metric in space–time, we scaled the time dimension, as suggested in [2], by a factor of , where res(z) is the throughplane resolution of the image and TR is the repetition time.
The algorithm involves iterative backprojection of the motion scattered slices onto a regular 4D grid using a 4D psf, estimating the slice intensities given the current estimate of the 4D image (forward projection using 4D psf), and computing the residual error between these two. This error is backprojected and added to the regularization term which was computed using a first order Tikhonov function. Optimization terminates when the residual error between estimated and measured data converges.
2.1 Implementation details
We implemented our algorithm of 4D reconstruction based on the publicly available toolkit of Numerical Solver Library (NSol) which is a part of NiftyMIC package [9] for volumetric reconstruction of the structural 3D images. A binary brain mask was manually delineated on the average volume of each fetus and dilated to ensure it covers the fetal brain through all ranges of the motion. A four dimensional estimate of the bias field for spatiotemporal signal nonuniformity correction in fMRI series was obtained using N4ITK algorithm [10] as suggested previously [11]. The motion correction step involved globally coregistering all time points, creating a target volume as proposed in [12] by automatically finding a set of consecutive volumes of fetal quiescence and averaging over them, and finally performing hierarchical slicetovolume registration based on the interleaved factor of acquisition as implemented in NiftyMIC package[9].
3 Experimental Results
We evaluated the proposed algorithm on fetal fMRI scans of 15 subjects between 22 and 39 weeks of gestation (Table.1). None of the cases showed any neurological pathologies. Pregnant women were scanned on a 1.5T clinical scanner (Philips Medical Systems, Best, Netherlands) using singleshot echoplanar imaging (EPI), and a sensitivity encoding (SENSE) cardiac coil with five elements. Image matrix size was 144144, with 1.741.74 inplane resolution, 3 slice thickness, flip angle of 90, and 96 volumes per acquisition. We used a hierarchical approach based on a symmetric block matching registration algorithm to estimate the rigid realignment parameters of each slice. Table.1 summarizes the estimated motional parameters for all subjects, showing especially free excessive rotation in fMRI images of the fetal population. Figure 1 shows the slicewise realignment parameters for a subject with strong rotation (27.9°in leftright direction). When such big motion occurs, gaps might be open up in the reconstructed image between slices that are spatially adjacent. Therefore, interpolation methods using only spatial information, although faster in implementation, may not correct it as there would be no data contributing to the reconstruction of a given regular grid. As fMRI data is temporally correlated, adjacent slices in time can potentially add information to the reconstruction. Here we used a four dimensional Gaussian kernel as psf to take both spatial and temporal structure of the data into account, and iteratively optimized our reconstruction through inverse problem framework.
GA  Max translation (mm)  Max rotation (°)  Sharpness  Standard deviation  

(week+day)  x  y  z  roll  pitch  yaw  Raw  Linear  Ours  Raw  Linear  Ours  
S1  22w6d  4.1  6.1  11.4  25.2  28.2  9.5  16949  11721  15729  1.60  7.49  5.78 
S2  23w4d  5.3  7.4  8  6.8  3.3  3.4  20409  15143  16119  18.87  45.08  48.30 
S3  25w0d  1.6  9.1  3.8  8.2  0.7  1.9  4672  5196  5255  6.37  6.82  7.34 
S4  25w4d  2.9  3.4  5  5.7  7.6  9.1  30539  31353  31010  18.59  35.29  5.06 
S6  28w6d  2.5  1.7  7.6  3.4  6.3  1.5  30855  28746  31298  16.34  25.03  12.22 
S7  29w3d  20.7  6  6.2  5.4  17.1  9.5  32761  31158  32807  19.37  53.90  62.28 
S8  29w4d  2.7  3.8  2.6  16.4  14.7  21.6  32299  30993  33154  43.07  58.17  21.75 
S9  29w5d  4.2  5.8  2.7  5.6  4.5  3.4  33218  32136  33368  46.41  71.01  19.87 
S10  30w2d  4.7  5.6  4.3  14.3  2  21.6  22531  21854  21653  68.15  74.28  21.79 
S11  32w2d  9.2  5.2  4.8  8.4  18.9  12.5  33218  32136  33368  8.16  8.35  4.68 
S12  34w4d  2.3  5  8.3  17.6  12.2  9.8  33198  32230  32744  22.42  81.04  32.02 
S13  34w6d  2.8  1.5  3.8  1.4  4.6  1.7  696  670  1233  9.86  11.79  6.86 
S14  35w6d  11.9  3  3.2  4.2  27.9  7.7  1035  785  1873  16.44  15.10  8.19 
S15  39w2d  1.4  3.7  1.5  3.6  3.8  1.3  803  909  1418  8.89  11.06  8.04 
3.1 Reconstruction accuracy
We compared the proposed 4D reconstruction with standard 3D linear reconstruction. We evaluated the sharpness [13] of the averaged volume after motion correction and reconstruction, and the standard deviation of BOLD signal fluctuations at each voxel over time. Figure 2 illustrates the average volume (top row), and standard deviation of intensity over time (bottom row) of the original data without motion correction, data after single step 3D linear interpolation, and data reconstructed using the proposed algorithm for two subjects with different motion characteristics. Severe blurring is observed in the image obtained from 3D linear interpolation. In contrast, the proposed method yields qualitatively sharper images as the effect of slice motion is appropriately decreased. Quantitative results are reported in Table 1 for all 15 subjects. Figure 3 shows the boxplot for the relative changes of the sharpness (left) and standard deviation (right) in 3D interpolation vs 4D Gaussian reconsruction normalized to the corresponding value in raw data. In 12 out of 14 subjects the proposed method yields higher sharpness than 3D interpolation, with an overall average sharpness higher than for 3D interpolation, in 11 it exhibits lower standard deviation, with an average standard deviation reduction of compared to 3D interpolation. For several subjects with lower gestational age the highest sharpness and lowest standard deviation occurs with uncorrected images. One possible cause is that the advanced gyrification at later gestational ages results in a pronounced sharpness advantage with correction, while for younger fetuses this is not the case and might reflect image noise rather than anatomical structure.
4 Conclusion
In this work we proposed a novel spatiotemporal reconstruction approach for the fetal brain acquired while there is unconstrained motion of the head. Our technique integrates the temporal continuity of fMRI time series with the spatial coherency of the neighboring brain voxels by using a four dimensional Gaussian point spread function for reconstructing motion scattered slices. Qualitative results of the resulting reconstruction of clinical fetal fMRI data shows improvement by the method. Quantitative evaluation confirms this improvement predominantly for higher gestational age. As the developed formulation for fMRI time series reconstruction is general and can be effectively used in adult studies as well, it would be of interest to investigate its performance on the adult population in the future.
5 Acknowledgment
This work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No 765148, and partial funding from the Austrian Science Fund (FWF, P35189). This work was also supported by the Swiss National Science Foundation (project 205321182602). We acknowledge access to the expertise of the CIBM Center for Biomedical Imaging, a Swiss research center of excellence founded and supported by Lausanne University Hospital (CHUV), University of Lausanne (UNIL), Ecole polytechnique fédérale de Lausanne (EPFL), University of Geneva (UNIGE) and Geneva University Hospitals (HUG).
References
 [1] Mark Jenkinson, Peter Bannister, Michael Brady, and Stephen Smith, “Improved optimization for the robust and accurate linear registration and motion correction of brain images,” Neuroimage, vol. 17, no. 2, pp. 825–841, 2002.
 [2] Giulio Ferrazzi, Maria Kuklisova Murgasova, Tomoki Arichi, Christina Malamateniou, Matthew J Fox, Antonios Makropoulos, Joanna Allsop, Mary Rutherford, Shaihan Malik, Paul Aljabar, et al., “Resting state fmri in the moving fetus: a robust framework for motion, bias field and spin history correction,” Neuroimage, vol. 101, pp. 555–568, 2014.
 [3] Saige Rutherford, Pascal Sturmfels, Mike Angstadt, Jasmine Hect, Jenna Wiens, Marion I van den Heuval, Dustin Scheinost, Moriah Thomason, and Chandra Sripada, “Observing the origins of human brain development: Automated processing of fetal fmri,” bioRxiv, p. 525386, 2019.
 [4] Wonsang You, Iordanis E Evangelou, Zungho Zun, Nickie Andescavage, and Catherine Limperopoulos, “Robust preprocessing for stimulusbased functional mri of the moving fetus,” Journal of Medical Imaging, vol. 3, no. 2, pp. 026001, 2016.
 [5] Daniel Sobotka, Roxane Licandro, Michael Ebner, Ernst Schwartz, Tom Vercauteren, Sebastien Ourselin, Gregor Kasprian, Daniela Prayer, and Georg Langs, “Reproducibility of functional connectivity estimates in motion corrected fetal fmri,” in Smart Ultrasound Imaging and Perinatal, Preterm and Paediatric Image Analysis, pp. 123–132. Springer, 2019.
 [6] Ali Gholipour, Judy A Estroff, and Simon K Warfield, “Robust superresolution volume reconstruction from slice acquisitions: application to fetal brain mri,” IEEE transactions on medical imaging, vol. 29, no. 10, pp. 1739–1758, 2010.

[7]
Maria KuklisovaMurgasova, Gerardine Quaghebeur, Mary A Rutherford, Joseph V
Hajnal, and Julia A Schnabel,
“Reconstruction of fetal brain mri with intensity matching and complete outlier removal,”
Medical image analysis, vol. 16, no. 8, pp. 1550–1564, 2012.  [8] Sébastien Tourbier, Xavier Bresson, Patric Hagmann, JeanPhilippe Thiran, Reto Meuli, and Meritxell Bach Cuadra, “An efficient total variation algorithm for superresolution in fetal brain mri with adaptive regularization,” NeuroImage, vol. 118, pp. 584–597, 2015.
 [9] Michael Ebner, Guotai Wang, Wenqi Li, Michael Aertsen, Premal A Patel, Rosalind Aughwane, Andrew Melbourne, Tom Doel, Steven Dymarkowski, Paolo De Coppi, et al., “An automated framework for localization, segmentation and superresolution reconstruction of fetal brain mri,” NeuroImage, vol. 206, pp. 116324, 2020.
 [10] Nicholas J Tustison, Brian B Avants, Philip A Cook, Yuanjie Zheng, Alexander Egan, Paul A Yushkevich, and James C Gee, “N4itk: improved n3 bias correction,” IEEE transactions on medical imaging, vol. 29, no. 6, pp. 1310–1320, 2010.
 [11] Sharmishtaa Seshamani, Xi Cheng, Mads Fogtmann, Moriah E Thomason, and Colin Studholme, “A method for handling intensity inhomogenieties in fmri sequences of moving anatomy of the early developing brain,” Medical image analysis, vol. 18, no. 2, pp. 285–300, 2014.
 [12] Sharmishtaa Seshamani, Mads Fogtmann, Xiaoyin Cheng, M Thomason, Chris Gatenby, and Colin Studholme, “Cascaded slice to volume registration for moving fetal fmri,” in 2013 IEEE 10th International Symposium on Biomedical Imaging. IEEE, 2013, pp. 796–799.

[13]
José Luis PechPacheco, Gabriel Cristóbal, Jesús ChamorroMartinez,
and Joaquín FernándezValdivia,
“Diatom autofocusing in brightfield microscopy: a comparative
study,”
in
Proceedings 15th International Conference on Pattern Recognition. ICPR2000
. IEEE, 2000, vol. 3, pp. 314–317.
Comments
There are no comments yet.