  Having accurate completeness functions is crucial to the determination of the
rest-frame ultraviolet luminosity functions (UVLFs) all the way back to the
epoch of reionization. Most studies use injection-recovery simulations to
determine completeness functions. Although conceptually similar, published
approaches have subtle but important differences in their definition of the
completeness function. As a result, they implement different methods to
determine the UVLFs. We discuss the advantages and limitations of existing
methods using a set of mock observations, and then compare the methods when
applied to the same set of Hubble Legacy Field (HLF) images. We find that the
most robust method under all our mock observations is the one that defines
completeness as a function of both input and output magnitude. Other methods
considering completeness only as a function of either input or output magnitude
may suffer limitations in a presence of photometric scatter and/or steep
luminosity functions. In particular, when the flux scatter is >0.2 mag, the
bias in the bright end of the UVLFs is on par with other systematic effects
such as the lensing magnification bias. When tested on HLF images, all methods
yield UVLFs that are consistent within 2 sigma confidence, suggesting that UVLF
uncertainties in the literature are still dominated by small number statistics
and cosmic variance. The completeness simulation code used in this study
(GLACiaR2) is publicly released with this paper as a tool to analyze future
higher precision datasets such as those expected from the James Webb Space
Telescope.
