The Monte Carlo code MCML (Monte Carlo modeling of light transport in multi-layered tissue) has been the gold standard for simulations of light transport in multi-layer tissue, but it is ineffective in the presence of three-dimensional (3D) heterogeneity. New techniques have been attempted to resolve this problem, such as MCLS, which is derived from MCML, and tMCimg, which draws upon image datasets. Nevertheless, these approaches are insufficient because of their low precision or simplistic modeling. We report on the development of a novel model for photon migration in voxelized media (MCVM) with 3D heterogeneity. Voxel crossing detection and refractive-index-unmatched boundaries were considered to improve the precision and eliminate dependence on refractive-index-matched tissue. Using a semi-infinite homogeneous medium, steady-state and time-resolved simulations of MCVM agreed well with MCML, with high precision (～100%) for the total diffuse reflectance and total fractional absorption compared to those of tMCimg (<70%). Based on a refractive-index-matched heterogeneous skin model, the results of MCVM were found to coincide with those of MCLS. Finally, MCVM was applied to a two-layered sphere with multi-inclusions, which is an example of a 3D heterogeneous media with refractive-index-unmatched boundaries. MCVM provided a reliable model for simulation of photon migration in voxelized 3D heterogeneous media, and it was developed to be a flexible and simple software tool that delivers high-precision results.