Efficient data visualization techniques are critical for many scientific applications. Centroidal Voronoi tessellation (CVT) based algorithms offer a convenient vehicle for performing image analysis, segmentation and compression while allowing to optimize retained image quality with respect to a given metric. In experimental science with data counts following Poisson distributions, several CVT-based data tessellation algorithms have been recently developed. Although they surpass their predecessors in robustness and quality of reconstructed data, time consumption remains to be an issue due to heavy utilization of the slowly converging Lloyd iteration. This paper discusses one possible approach to accelerating data visualization algorithms. It relies on a multidimensional generalization of the optimization based multilevel algorithm for the numerical computation of the CVTs introduced in [1], where a rigorous proof of its uniform convergence has been presented in 1-dimensional setting. The multidimensional implementation employs barycentric coordinate based interpolation and maximal independent set coarsening procedures. It is shown that when coupled with bin accretion algorithm accounting for the discrete nature of the data, the algorithm outperforms Lloyd-based schemes and preserves uniform convergence with respect to the problem size. Although numerical demonstrations provided are limited to spectroscopy data analysis, the method has a context-independent setup and can potentially deliver significant speedup to other scientific and engineering applications.