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 EPSFBuilder, 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)

epsf_builder = EPSFBuilder(oversampling=4, maxiters=3,
                           progress_bar=False)
epsf, fitted_stars = epsf_builder(stars)

fig, ax = plt.subplots(figsize=(8, 8))
norm = simple_norm(epsf.data, 'log', percent=99.0)
axim = ax.imshow(epsf.data, norm=norm, origin='lower')
fig.colorbar(axim)