In [72]:

    

%matplotlib inline

import numpy as np
from astropy.io import fits
from time import time as tt
import matplotlib.pyplot as plt
from skimage import img_as_float, img_as_uint, img_as_ubyte
from skimage.filters import median, gaussian
from skimage.exposure import rescale_intensity
from skimage.morphology import remove_small_objects
from skimage.morphology.selem import disk
from skimage.transform import resize
from skimage.morphology import white_tophat, opening
from skimage.filters import threshold_otsu, threshold_li, threshold_yen


    

In [86]:

    

# Load and downsample the test image
fname = 'M16_Halpha.fit'
with fits.open(fname) as hdu_list:
    img = hdu_list[0].data

ny, nx = img.shape

# FWHM in pixels (from separate estimate)
FWHM = 7.4

# Matched Gaussian filter (sigma = FWHM/2)
sigma_g = FWHM * 0.5
img = gaussian(img, sigma_g)

# Resampling scale factor to give FWHM = 4 pixels
sf = 4.0 / FWHM

# Downsample image (biquad, antialiasing)
nxd = int(nx * sf)
nyd = int(ny * sf)
imgd = resize(img, [nyd, nxd], order=3, anti_aliasing=False, mode='reflect')

print('Resized dtype  : %s' % imgd.dtype)
print('Resized size   : %d x %d' % (nxd, nyd))


    

    




Resized dtype  : float64
Resized size   : 1673 x 1124

In [93]:

    

# Structuring element on scale of typical star
# Radius = 5 works well after matched downsampling (empirical)
star_selem = disk(radius=5)

# White tophat filter
# - suppress smooth background
# - highlight bright objects smaller than selem
imgd_wth = white_tophat(imgd, star_selem)

# Global threshold and connected regions
star_mask = imgd_wth > threshold_yen(imgd_wth)

# Remove small objects (area < 5)
star_mask = remove_small_objects(star_mask, min_size=5)


    

In [96]:

    

# ROI
x0f, x1f = 0.40, 0.58
y0f, y1f = 0.5, 0.75
x0, x1 = x0f * nxd, x1f * nxd
y0, y1 = y0f * nyd, y1f * nyd

# Plot white tophat filtered star field
fig, ax = plt.subplots(1, 3,figsize=(16,16))

ax[0].imshow(imgd, origin='lower')
ax[0].set_xlim(x0, x1)
ax[0].set_ylim(y0, y1)

ax[1].imshow(imgd_wth, origin='lower')
ax[1].set_xlim(x0, x1)
ax[1].set_ylim(y0, y1)

ax[2].imshow(star_mask, origin='lower')
ax[2].set_xlim(x0, x1)
ax[2].set_ylim(y0, y1)

plt.show()