Traditional statistical wavelet analysis carries out modeling and inference under a given, predetermined wavelet transform. This approach can quickly lose efficiency for multi-dimensional data (e.g., observations measured on a multi-dimensional grid), because a predetermined wavelet transform does not adaptively exploit the information or energy distribution in a problem-specific manner. This work aims to overcome this challenge within the multivariate wavelet analysis framework by incorporating adaptivity into the wavelet transform itself in a principled manner. By exploiting a connection between wavelet transforms and permutations on the index space of multi-dimensional functions, we show that the desired adaptive wavelet transform can be achieved by adopting a layer of Bayesian hierarchical modeling on the space of such permutations. In particular, when combined with the Haar basis, exact Bayesian inference under the model can be achieved analytically through a recursive message passing algorithm with an efficient computational complexity that is linear with sample size. We also provide recipe for incorporating block shrinkage and general wavelet bases into the framework, all while maintaining such adaptivity. We demonstrate via extensive numerical experiments that with our framework even simple 1D Haar wavelets can achieve excellent performance in the context of 2D and 3D image reconstruction, outperforming state-of-the-art wavelet and non-wavelet methods especially in noisy, low signal-to-noise ratio settings at a fraction of the computational cost. Furthermore, we investigate the source of the gain by quantitatively comparing the efficacy of energy concentration under our adaptive wavelet transform with that of classical fixed wavelet transforms.