In the fields of 3D modeling, analysis of discontinuities and engineering calculation, surface extraction is of great importance. The rapid development of photogrammetry and Light Detection and Ranging (LiDAR) technology facilitates the study of surface extraction. Automatic extraction of rock surfaces from 3D rock-mass point clouds also becomes the basis of 3D modeling and engineering calculation of rock mass. This paper presents an automated and effective method for extracting rock surfaces from unorganized rock-mass point clouds. This method consists of three stages: (i) clustering based on voxels; (ii) estimating major orientations based on Gaussian Kernel and (iii) rock surface extraction. Firstly, the two-level spatial grid is used for fast voxelization and segmenting the point cloud into three types of voxels, including coplanar, non-coplanar and sparse voxels. Secondly, the coplanar voxels, rather than the scattered points, are employed to estimate major orientations by using a bivariate Gaussian Kernel. Finally, the seed voxels are selected on the basis of major orientations and the region growing method based on voxels is applied to extract rock surfaces, resulting in sets of surface clusters. The sub-surfaces of each cluster are coplanar or parallel. In this paper, artificial icosahedron point cloud and natural rock-mass point clouds are used for testing the proposed method, respectively. The experimental results show that, the proposed method can effectively and accurately extract rock surfaces in unorganized rock-mass point clouds.