import matplotlib.pyplot as plt
import numpy as np
from photutils.datasets import make_noise_image
from photutils.psf import (CircularGaussianPRF, SourceGrouper,
                           make_psf_model_image)

shape = (256, 256)
psf_shape = (11, 11)
border_size = (7, 7)
flux = (500, 1000)
fwhm = 4.7
psf_model = CircularGaussianPRF(fwhm=fwhm)
n_sources = 100
data, stars = make_psf_model_image(shape, psf_model, n_sources,
                                   flux=flux,
                                   model_shape=psf_shape,
                                   border_size=border_size, seed=123)
noise = make_noise_image(shape, mean=0, stddev=2, seed=123)
data += noise

min_separation = 2.5 * fwhm
grouper = SourceGrouper(min_separation)

x = np.array(stars['x_0'])
y = np.array(stars['y_0'])
groups = grouper(x, y, return_groups_object=True)

fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(data, origin='lower')
groups.plot(radius=fwhm, ax=ax, lw=2, seed=123)