import matplotlib.pyplot as plt
from astropy.nddata import NDData
from astropy.stats import sigma_clipped_stats
from astropy.table import Table
from astropy.visualization import simple_norm
from photutils.datasets import (load_simulated_hst_star_image,
                                make_noise_image)
from photutils.detection import DAOStarFinder
from photutils.psf import extract_stars

hdu = load_simulated_hst_star_image()
data = hdu.data
data += make_noise_image(data.shape, distribution='gaussian', mean=10.0,
                         stddev=5.0, seed=0)
finder = DAOStarFinder(threshold=100.0, fwhm=1.5)
sources = finder(data)

size = 25
hsize = (size - 1) / 2
x = sources['x_centroid']
y = sources['y_centroid']
mask = ((x > hsize) & (x < (data.shape[1] - 1 - hsize))
        & (y > hsize) & (y < (data.shape[0] - 1 - hsize)))
stars_tbl = Table()
stars_tbl['x'] = x[mask]
stars_tbl['y'] = y[mask]

mean_val, median_val, std_val = sigma_clipped_stats(data, sigma=2.0)
data -= median_val

nddata = NDData(data=data)

stars = extract_stars(nddata, stars_tbl, size=25)

nrows = 5
ncols = 5
fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(20, 20),
                       squeeze=True)
ax = ax.ravel()
for i in range(nrows * ncols):
    norm = simple_norm(stars[i], 'log', percent=99.0)
    ax[i].imshow(stars[i], norm=norm, origin='lower')