Geophysical site investigation techniques based on elastic waves have been widely used to characterize rock masses. However, characterizing jointed rock masses by using such techniques remains challenging because of a lack of knowledge about elastic wave propagation in multi-jointed rock masses. In this paper, the roughness of naturally fractured rock joint surfaces is estimated by using a three-dimensional (3D) image-processing technique. The classification of the joint roughness coefficient (JRC) is enhanced by introducing the scan line technique. The peak-to-valley height is selected as a key indicator for JRC classification. Long-wavelength P-wave and torsional S-wave propagation across rock masses containing naturally fractured joints are simulated through the quasi-static resonant column (QSRC) test. In general, as the JRC increases, the S-wave velocity increases within the range of stress levels considered in this paper, whereas the P-wave velocity and the damping ratio of the shear wave decrease. In particular, the two-dimensional joint specimen underestimates the S-wave velocity while overestimating the P-wave velocity. This suggests that 3D joint surfaces should be implicated to obtain the reliable elastic wave velocity in jointed rock masses. The contact characteristic and degree of roughness and waviness of the joint surface are identified as a factor influencing P-wave and S-wave propagation in multi-jointed rock masses. The results indicate a need for a better understanding of the sensitivity of contact area alterations to the elastic wave velocity induced by changes in normal stress. This paper's framework can be a reference for future research on elastic wave propagation in naturally multi-jointed rock masses.