We have calculated the direct correlation function, c(1, 2), via a high-order virial expansion, for systems of hard spheres and spheroids, in both the isotropic and nematic phases. For hard spheres, we find that truncation at sixth order in density gives good agreement with simulation data. Close to freezing, the virial series still appears to converge to the simulation results, but there are significant discrepancies, particularly at very small separations. In the non-overlap region, the virial theory begins to capture the features found from simulation. We also calculate the pair distribution function, g(1, 2), from our estimates of c(1, 2), and find good agreement with simulation data at all densities up to freezing. For hard, prolate spheroids of aspect ratio 3 : 1, we calculate c(1, 2) and g(1, 2) in the isotropic phase, again finding good agreement with simulation data at moderate densities. Finally, we present the result of our calculations on c(1, 2) for 3 : 1 spheroids in the nematic phase.