Synopsis

State-of-the-art techniques for denoising functional MRI (fMRI) images consider the problems of spatial and temporal regularization as decoupled tasks. In this work we propose a partial differential equations (PDEs) -based algorithm that acts directly on the 4-D fMRI image. Our approach is based on the idea that large image variations should be preserved as they occur during brain activation, but small variations should be smoothed to remove noise. Starting from this principle, by means of PDEs we were able to smooth the fMRI image with an anisotropic regularization, thus recovering the location of the brain activations in space and their timing and duration.

Introduction

Methods

Let $$$Ω$$$ be a 4-D (3-D space × time) domain and assume Neumann boundary conditions on $$$\partial Ω$$$. We defined a regularization process, based on gradient descent computed with PDEs, such that:

$$\frac{\partial I}{\partial t}=(1- \lambda) \frac{\bar H \star(I_0-H \star I)}{A}+ \lambda \frac{div(D \nabla I)}{B} \quad \text{(1)}$$

where $$$I_0$$$ and $$$I$$$ are the original and the denoised image respectively, $$$A=\|I_0\|$$$ and $$$B= \|div(D \nabla I_0)\|$$$ are normalization factors, $$$H$$$ is the hemodynamic response function⁴ (HRF), $$$\bar H$$$ is the time-reversed HRF, $$$\lambda$$$ is the regularization parameter and $$$D = \frac{\nabla I \nabla I^T}{\| \nabla I\|^2} \star G$$$ is the 4-D structure tensor of $$$I$$$ smoothed by the gaussian kernel $$$G$$$ ($$$\sigma _G=1$$$). In (1) the first term on the right is the data fitting term, which measures the correlation of the residual with $$$H$$$, and the second term minimizes image variations. The convolution with $$$H$$$ and $$$\bar H$$$ were computed only along the time dimension. After computing the operator $$$D$$$, we defined the directions of the image variations by an eigendecomposition of $$$D$$$ such that $$$D = Q \Lambda Q^T$$$, where $$$Q$$$ contains the orthogonal eigenvectors ($$$\theta_{1,2,3,4}$$$) of $$$D$$$ and $$$\Lambda$$$ contains their associated eigenvalues ($$$\lambda_1> \lambda_2> \lambda_3> \lambda_4 $$$). Finally, we recomputed the matrix $$$\tilde D = Q \tilde \Lambda Q^T$$$ as follows: the highest eigenvalue, $$$\lambda_1$$$, was set according to a function of $$$\| \nabla I \|$$$, i.e. $$$f(\|\nabla I\|)$$$, such that if $$$\| \nabla I \| >> 0$$$ the current voxel may be located on a edge and the smoothing is performed just along the other three directions (anisotropic smoothing). Otherwise, if $$$\| \nabla I \| \rightarrow0$$$ the smoothing will be isotropic in all the four directions. $$$\lambda_2$$$,$$$\lambda_3$$$,$$$\lambda_4$$$ were indeed set to 1 to perform an isotropic smoothing. After that, we replaced the operator $$$D$$$ in (1) with $$$\tilde D$$$.

As a proof of concept, similarly to previous works^5,6, we scaled a 3-D activation map computed with the FMRIB Software Library (FSL) in the range [0,3], with a 2-mm isotropic resolution (Fig.1.a). We multiplied it by a piece-wise constant signal of 100s, with one onset of 40s (Fig.1.b). We corrupted the image with model noise, we convolved it with the HRF⁴ and we added gaussian noise thus simulating the fMRI time-courses $$$y(t)$$$. We regularized the image as showed above, and we recovered $$$u^*(t)$$$. Finally, to evaluate the results, we computed the root of the mean square errors (MSE) and standard deviation (STD) and the Pearson correlation and its STD between $$$u(t)$$$ and $$$u^*(t)$$$ averaged among the voxels belonging to the grey matter. We compared our results with those obtained using the Total Activation approach⁶ (TA), implemented in the iCAPs toolbox (https://miplab.epfl.ch/index.php/software/total-activation).

Results

Fig.2 shows examples of regularized spatial maps and time series using our approach ($$$u^*_{PDEs}$$$) and the TA ($$$u^*_{TA}$$$). Fig.3 shows that the MSEs $$$\pm$$$ STDs change for different peak-SNRs (pSNRs) and that they are lower than the ones obtained using TA. Fig.4 shows that the activation recovered with our approach is more correlated with the ground truth, for different pSNRs, whereas the results obtained with TA are more sensitive to noise.

Discussion

Both techniques succeeded at finding the activated regions, but the results obtained using our approach were closer to the ground truth in terms of amplitude and also of correlation between the simulated activations and the recovered ones, for different pSNRs. Our findings shows that our approach enabled us to solve a unique problem, coupling the space and the time dimension thus having to set just one regularization parameter, $$$\lambda$$$, rather than dividing it into two separate problems, i.e. time and space, as it was done by Farouj and colleagues⁶.

Conclusion

Acknowledgements

References

1. Tschumperle, D., & Deriche, R. (2002). Diffusion PDEs on vector-valued images. IEEE Signal Processing Magazine, 19(5), 16-25.
2. Tschumperle, D., & Deriche, R. (2005). Vector-valued image regularization with PDEs: A common framework for different applications. IEEE transactions on pattern analysis and machine intelligence, 27(4), 506-517.
3. Karahanoğlu, F. I., Caballero-Gaudes, C., Lazeyras, F., & Van De Ville, D. (2013). Total activation: fMRI deconvolution through spatio-temporal regularization. Neuroimage, 73, 121-134.
4. Khalidov, I., Fadili, J., Lazeyras, F., Van De Ville, D., & Unser, M. (2011). Activelets: Wavelets for sparse representation of hemodynamic responses. Signal Processing, 91(12), 2810-2821.
5. Costantini, I., Filipiak, P., Maksymenko, K., Deslauriers-Gauthier, S., & Deriche, R. (2018, July). fMRI Deconvolution via Temporal Regularization using a LASSO model and the LARS algorithm. In 40th International Engineering in Medicine and Biology Conference.
6. Farouj, Y., Karahanoğlu, F. I., & Van De Ville, D. (2017, April). Regularized spatiotemporal deconvolution of fMRI data using gray-matter constrained total variation. In Biomedical Imaging (ISBI 2017), 2017 IEEE 14th International Symposium on (pp. 472-475). Ieee.