We have developed a novel method of determining 2D radial density profiles for astronomical systems of discrete objects using Voronoi tessellations. This Voronoi-based method was tested against the standard annulus-based method on five simulated systems of objects, following known Hubble density profiles of varying parameters and sizes. It was found that the Voronoi-based method returned radial density fits with lower uncertainties on the fitting parameters across all five systems compared to the annulus-based method. The Voronoi-based method also consistently returned more accurate estimates of the total number of objects in each system than the annulus-based method, and this accuracy increased with increasing system size. Finally, the Voronoi-based method was applied to two observed globular cluster systems around brightest cluster galaxies ESO 444-G046 and 2MASX J13272961-3123237 and the results were compared to previous results for these galaxies obtained with the annulus-based method. Again, it was found that the Voronoi-based method returned fits with lower uncertainties on the fitting parameters, and the total number of globular clusters returned are within errors of the annulus-based method estimates, however also with lower uncertainties.