Abstract: TR Centroidal Voronoi tessellation (CVT) is a fundamental geometric structure that finds many applications in computational science and engineering, including computer graphics. The prevailing method for computing CVT is Lloyd's method, which has linear convergence and is inefficient in practice. Our goal is to develop efficient methods for CVT computation, justify the fast convergence of these methods theoretically and demonstrate their superiority with experimental examples in various cases. Specifically, it is shown that the CVT energy function has C2 smoothness in convex domains and in most other commonly encountered domains with smooth density, correcting the view in the literature that this function is non-smooth (that is, merely C0 but not C1 ). Due to its C2 smoothness, it is therefore possible to minimize the CVT energy functions using Newton-like optimization methods and expect fast convergence. We apply quasi-Newton methods to computing CVT and demonstrate their faster convergence than Lloyd's method and their better robustness than the Lloyd-Newton method, a previous attempt at CVT computation acceleration. The application of these results to surface remeshing in computer graphics is also studied.