Given the spline representation of the boundary of a three dimensional domain, constructing a volumetric spline parameterization of the domain (i.e., a map from a unit cube to the domain) with the given boundary is a fundamental problem in isogeometric analysis. A good domain parameterization should satisfy the following criteria: (1) the parameterization is a bijective map; and (2) the map has lowest possible distortion. However, none of the state-of-the-art volumetric parameterization methods has fully addressed the above issues. In this paper, we propose a three-stage approach for constructing volumetric parameterization satisfying the above criteria. Firstly, a harmonic map is computed between a unit cube and the computational domain. Then a bijective map modeled by a max-min optimization problem is computed in a coarse-to-fine way, and an algorithm based on divide and conquer strategy is proposed to solve the optimization problem efficiently. Finally, to ensure high quality of the parameterization, the MIPS (Most Isometric Parameterizations) method is adopted to reduce the conformal distortion of the bijective map. We provide several examples to demonstrate the feasibility of our approach and to compare our approach with some state-of-the-art methods. The results show that our algorithm produces bijective parameterization with high quality even for complex domains.