Searches for continuous gravitational waves target nearly monochromatic gravitational wave emission from e.g. non-axysmmetric fast-spinning neutron stars. Broad surveys often require to explicitly search for a very large number of different waveforms, easily exceeding ∼ 10 17 templates. In such cases, for practical reasons, only the top, say ∼ 10 10 , results are saved and followed-up through a hierarchy of stages. Most of these candidates are not completely independent of neighbouring ones, but arise due to some common cause: a fluctuation, a signal or a disturbance. By judiciously clustering together candidates stemming from the same root cause, the subsequent follow-ups become more effective. A number of clustering algorithms have been employed in past searches based on iteratively finding symmetric and compact over-densities around candidates with high detection statistic values. The new clustering method presented in this paper is a significant improvement over previous methods: it is agnostic about the shape of the over-densities, is very efficient and it is effective: at a very high detection efficiency, it has a noise rejection of 99.99% , is capable of clustering two orders of magnitude more candidates than attainable before and, at fixed sensitivity it enables more than a factor of 30 faster follow-ups. We also demonstrate how to optimally choose the clustering parameters.