Supremum norm loss is intuitively more meaningful to quantify function estimation error in statistics. In the context of multivariate nonparametric regression with unknown error, we propose a Bayesian procedure based on spike-and-slab prior and wavelet projections to estimate the regression function and all its mixed partial derivatives. We show that their posterior distributions contract to the truth optimally and adaptively under supremum-norm loss. The master theorem through tests with exponential errors used in Bayesian nonparametrics was not adequate to deal with this problem, and we developed a new idea such that posterior under the regression model is systematically reduced to a posterior arising from some quasi-white noise model, where the latter model greatly simplifies our rate calculations. Hence, this paper takes the first step in showing explicitly how one can translate results from white noise to regression model in a Bayesian setting.