Create a Delauney triangulation and assign a quantity (for example the particle density) to each vertex. From this locate the local minimas, i.e. vertices where all of the neighbors have a larger value of the quantity. Then we fill "water" into the minimas until they "overflow" into a nearby basin. This defines watershed basins: a collection of particles that belong to a given minima. We can also go further and continue filling water into the basins until they overflow into a basin with lower density and thereby get larger and larger basins (this is a simple merging that I've not bothered doing yet). This can be used as a void and halo finder (which also gives us sub-voids and sub-halos).
Figure 1: Watershed regions
// A general watershed algorithm template<class T, class U> void WatershedGeneral( MPIPeriodicDelaunay<T, VertexDataWatershed> & PdT, T *p, size_t NumPart, std::vector
& quantity, std::vector<U> & watershed_groups); // Same as above but where quantity is the density (mass of particle over Voronoi volume) // and where we create the triangulation within the method template<class T, class U> void WatershedDensity( T *p, size_t NumPart, std::vector<U> & watershed_groups, double buffer_fraction = 0.30, double random_fraction = 0.5); bool do_density_maximum = false);