This paper presents a novel spatial Bayesian method for simultaneous activation detection and hemodynamic response function (HRF) estimation of functional magnetic resonance imaging (fMRI) data. A Bayesian variable selection approach is used to induce shrinkage and sparsity, with a spatial prior on latent variables representing activated hemodynamic response components. Then, the activation map is generated from the full spectrum of posterior inference constructed through a Markov chain Monte Carlo scheme, and HRFs at different voxels are estimated non-parametrically with information pooling from neighboring voxels. By integrating functional activation detection and HRFs estimation in a unified framework, our method is more robust to noise and less sensitive to model mis-specification.