In the preceding paper, a rigorous three-dimensional relativistic equation for two-gluon bound states was derived from QCD with massive gluons and represented in the angular momentum representation. In order to apply this equation to calculate the glueball spectrum, in this paper, the equation is recast in an equivalent three-dimensional relativistic equation satisfied by the two-gluon positive energy state amplitude. The interaction Hamiltonian in the equation is exactly derived and expressed as a perturbative series. The first term in the series describes the one-gluon exchange interaction which includes fully the retardation effect in it. This term plus the linear confining potential are chosen to be the interaction Hamiltonian and employed in the practical calculation. With the integrals containing three and four spherical Bessel functions in the QCD vertices being analytically calculated, the interaction Hamiltonian is given an explicit expression in the angular momentum representation. Numerically solving the relativistic equation with taking the contributions arising from the retardation effect and the longitudinal mode of gluon fields into account, a set of masses for the ${0}^{++}{,0}^{\ensuremath{-}+}{,1}^{++}{,1}^{\ensuremath{-}+}{,2}^{++},$ and ${2}^{\ensuremath{-}+}$ glueball states are obtained and are in fairly good agreement with the predictions given by the lattice simulation. In addition, some new glueball states are predicted.