Skip to content

Commit

Permalink
fixes for spacetelescope#165, spacetelescope#238, spacetelescope#210, s…
Browse files Browse the repository at this point in the history 
…pacetelescope#205
  • Loading branch information
@sosey
sosey committed Oct 18, 2021
1 parent 9a4e7a4 commit cd21c7a
Show file tree
Hide file tree
Showing 4 changed files with 70 additions and 57 deletions.
4 changes: 3 additions & 1 deletion CHANGES.rst
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
version 0.9.2.dev (unreleased)
------------------------------

- updated the gauss1d to set fixed fitting parameters and improve FWHM estimate(#165, #238)
- fixed background sky subtraction (#210)
- formatting fix for print statement in aper phot plot (#205)

version 0.9.1 (2020-03-30)
--------------------------
Expand Down
1 change: 1 addition & 0 deletions imexam/imexam_defpars.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@
"fitplot": [True, "Overplot model fit?"],
"func": ["Gaussian1D", "Model form to fit"],
"center": [True, "Solve for center using 2d Gaussian? [bool]"],
"delta": [10, "bounding box for centering measurement"],
"background": [True, "Subtract background? [bool]"],
"skyrad": [10., "Background inner radius in pixels, from center of object"],
"width": [5., "Background annulus width in pixels"],
Expand Down
94 changes: 48 additions & 46 deletions imexam/imexamine.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@ def __init__(self):
"Moffat1D",
"MexicanHat1D",
"AiryDisk2D",
"Gaussian2D",
"Polynomial1D"]

# see if the package logger was already started
Expand Down Expand Up @@ -633,7 +634,7 @@ def aper_phot(self, x, y, data=None,

# Construct the output strings (header and parameter values)
pheader = f"x\ty\tradius\tflux\tmag(zpt={magzero:0.2})\t"
pstr = f"\n{x:.2f}\t{y:0.2f}\t{radius}\t{total_flux:0.2}\t{mag:0.2}\t"
pstr = f"\n{xx:.2f}\t{yy:0.2f}\t{radius}\t{total_flux:0.2}\t{mag:0.2}\t"

if subsky:
pheader += "sky/pix\t"
Expand Down Expand Up @@ -680,7 +681,6 @@ def aper_phot(self, x, y, data=None,
self.aper_phot_pars['color_max'][0]]

pad = outer * 1.2 # XXX TODO: Bad magic number
print(xx, yy, pad)
ax.imshow(data[int(yy - pad):int(yy + pad),
int(xx - pad):int(xx + pad)],
vmin=color_range[0], vmax=color_range[1],
Expand Down Expand Up @@ -726,7 +726,7 @@ def line_fit(self, x, y, data=None, form=None, genplot=True, fig=None, col=False
The background is currently ignored
If centering is True in the parameter set, then the center
is fit with a 2d gaussian, not performed for Polynomial1D
is fit with a 2d gaussian, center is not performed for Polynomial1D
"""

# Set which parameters to use
Expand Down Expand Up @@ -949,8 +949,6 @@ def com_center(self, x, y, data=None, delta=None, oversampling=1.):
if delta >= len(data) / 4:
delta = delta // 2

if delta is None:
delta = int(self.com_center_pars['delta'][0])
xx = int(x)
yy = int(y)

Expand Down Expand Up @@ -1046,18 +1044,32 @@ def radial_profile(self, x, y, data=None, form=None,
"""
pars = self.radial_profile_pars

if data is None:
data = self._data

fitplot = bool(pars["fitplot"][0])
if fitplot:
if form is None:
form = getattr(models, pars["func"][0])
if form.name not in self._fit_models:
raise ValueError(f"Functional {form.name} not in available: {self._fit_models}")

self.log.info(f"using model: {form}")

subtract_background = bool(pars["background"][0])
if not photutils_installed and subtract_background:
self.log.warning("Install photutils to enable "
"background subtraction")
subtract_background = False

if data is None:
data = self._data


getdata = bool(pars["getdata"][0])
center = pars["center"][0]
delta = int(pars["delta"][0])

# reset delta for small arrays
if delta >= len(data) / 4:
delta = delta // 2

# be careful with the clipping since most
# of the data will be near the low value
Expand All @@ -1067,41 +1079,30 @@ def radial_profile(self, x, y, data=None, form=None,
else:
sig_factor = 0

fitplot = bool(pars["fitplot"][0])

if fitplot:
if form is None:
fitform = getattr(models, pars["func"][0])
else:
if form not in self._fit_models:
msg = f"{form} not supported for fitting"
self.log.info(msg)
raise ValueError(msg)
else:
fitform = getattr(models, form)

# cut the data down to size and center cutout
datasize = int(pars["rplot"][0])
if datasize < 3:
self.log.info("Insufficient pixels, resetting chunk size to 3.")
datasize = 3

if center:
# reset delta for small arrays
# make it odd if it's even
if ((datasize % 2) == 0):
datasize = datasize + 1
xx = int(x)
yy = int(y)
# flipped from xpa
if pars["center"][0]:
xx=int(x)
yy=int(y)
data_chunk = data[yy - datasize:yy + datasize,
xx - datasize:xx + datasize]
amp, centerx, centery, sigmax, sigmay = self.gauss_center(xx, yy, data, delta=datasize)

amp, xx, yy, sigma, sigmay = self.gauss_center(xx,
yy,
data,
delta=delta)
if (xx < 0 or yy < 0 or xx > data.shape[1] or
yy > data.shape[0]):
self.log.warning("Problem with centering, pixel coords")
else:
centerx = xx
centery = yy
else:
centery = y
centerx = x

centery = y
icenterx = int(centerx)
icentery = int(centery)

Expand Down Expand Up @@ -1157,7 +1158,7 @@ def radial_profile(self, x, y, data=None, form=None,
sky_per_pix = 0.
self.log.info("Sky background negative, setting to zero")
self.log.info(f"Background per pixel: {sky_per_pix}")
flux -= sky_per_pix
flux = flux - sky_per_pix

if getdata:
self.log.info(f"Sky per pixel: {sky_per_pix} using "
Expand All @@ -1167,17 +1168,18 @@ def radial_profile(self, x, y, data=None, form=None,
self.log.info(info)
self.log.info(radius, flux)


# Fit the functional form to the radial profile flux
# TODO: Ignore sky subtracted pixels that push flux
# below zero?
if fitplot:
fline = np.linspace(0, datasize, 100) # finer sample
# fit model to data
if fitform.name == "Gaussian1D":
if form.name == "Gaussian1D":# make center radii at max ********

fitted = math_helper.fit_gauss_1d(radius, flux,
sigma_factor=sig_factor,
center_at=0,
center_at=radius[0],
weighted=True)
print(f"fitline: {fitted}")

fwhmx, fwhmy = math_helper.gfwhm(fitted.stddev_0.value)
legend = (f"Max. pix. flux = {np.max(flux):9.3f}\n"
Expand All @@ -1187,10 +1189,10 @@ def radial_profile(self, x, y, data=None, form=None,
legendx = datasize / 2
legendy = np.max(flux) / 2

elif fitform.name == "Moffat1D":
elif form.name == "Moffat1D":
fitted = math_helper.fit_moffat_1d(flux,
sigma_factor=sig_factor,
center_at=0,
center_at=radius[0],
weighted=True)
mfwhm = math_helper.mfwhm(fitted.alpha_0.value,
fitted.gamma_0.value)
Expand All @@ -1200,10 +1202,10 @@ def radial_profile(self, x, y, data=None, form=None,
legendx = datasize / 2
legendy = np.max(flux) / 2

elif fitform.name == "MexicanHat1D":
elif form.name == "MexicanHat1D":
fitted = math_helper.fit_mex_hat_1d(flux,
sigma_factor=sig_factor,
center_at=0,
center_at=radius[0],
weighted=True)
legend = (f"Max. pix. flux = {np.max(flux):9.3f}\n")
legendx = datasize / 2
Expand Down Expand Up @@ -1242,14 +1244,14 @@ def radial_profile(self, x, y, data=None, form=None,

if pars["title"][0] is None:
if fitplot:
title = f"Radial Profile at ({icenterx},{icentery}) with {fitform.name}"
title = f"Radial Profile at ({icenterx},{icentery}) with {form.name}"
else:
title = f"Radial Profile for {icenterx} {icentery}"
else:
title = pars["title"][0]

if fitplot:
ax.plot(fline, yfit, linestyle='-', c='r', label=fitform.name)
ax.plot(fline, yfit, linestyle='-', c='r', label=form.name)
ax.set_xlim(0, datasize, 0.5)
ax.text(legendx, legendy, legend)

Expand Down Expand Up @@ -1622,7 +1624,7 @@ def surface(self, x, y, data=None, fig=None):
if fig is None:
fig = plt.figure(self._figure_name)
fig.clf()
fig.add_subplot(111)
fig.add_subplot(111, projection='3d')
ax = fig.gca(projection='3d')

title = self.surface_pars["title"][0]
Expand Down
28 changes: 18 additions & 10 deletions imexam/math_helper.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,10 +93,13 @@ def fit_moffat_1d(data, gamma=2., alpha=1., sigma_factor=0.,
ldata = len(data)

# guesstimate mean
if center_at:
x0 = 0
else:
if center_at is None:
x0 = int(ldata / 2.)
fluxmax = max(data)
else:
x0 = center_at
fluxmax = max(data)


# assumes negligable background
if weighted:
Expand All @@ -115,7 +118,7 @@ def fit_moffat_1d(data, gamma=2., alpha=1., sigma_factor=0.,
fit = fitter

# Moffat1D + constant
model = (models.Moffat1D(amplitude=max(data),
model = (models.Moffat1D(amplitude=fluxmax,
x_0=x0,
gamma=gamma,
alpha=alpha) +
Expand Down Expand Up @@ -165,12 +168,15 @@ def fit_gauss_1d(radius, flux, sigma_factor=0, center_at=None, weighted=False):
if radius.shape != flux.shape:
raise ValueError("Expected same sizes for radius and flux")

# guesstimate the mean
# assumes ordered radius

if center_at is None:
delta = int(len(radius) / 2.)
center_mean = sum(radius)/len(radius)
fixed = {'amplitude': False}
bounds = {'amplitude': (flux.max()/2.,flux.max())}
else:
delta = center_at
center_mean = center_at
fixed = {'mean': True, 'amplitude': True}
bounds = {'amplitude': (flux.max()/2., flux.max())}

# assumes negligable background
if weighted:
Expand All @@ -187,8 +193,10 @@ def fit_gauss_1d(radius, flux, sigma_factor=0, center_at=None, weighted=False):
fit = fitter

# Gaussian1D + a constant
model = (models.Gaussian1D(amplitude=flux.max() - flux.min(),
mean=delta, stddev=1.) +
model = (models.Gaussian1D(amplitude=flux.max(),
mean=center_mean, stddev=1.,
fixed=fixed,
bounds=bounds) +
models.Polynomial1D(c0=flux.min(), degree=0))

with warnings.catch_warnings():
Expand Down

0 comments on commit cd21c7a

Please sign in to comment.