diff --git a/CHANGES.rst b/CHANGES.rst index a624a80..ff96bb4 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -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) -------------------------- diff --git a/imexam/connect.py b/imexam/connect.py index 2cbad46..92e219c 100644 --- a/imexam/connect.py +++ b/imexam/connect.py @@ -233,7 +233,7 @@ def _run_imexam(self): # This loop now recognizes the arrow keys # for moving the cursor in the window, it calls # the windows cursor function. However, depending - # on how the use has their windowing focus setup + # on how the user has their windowing focus setup # they might loose focus on the DS9 window and have to # move the cursor to gain focus again, is there a way # around this? Cursor is not implemented in the ginga diff --git a/imexam/imexam_defpars.py b/imexam/imexam_defpars.py index df9f776..1132f0b 100644 --- a/imexam/imexam_defpars.py +++ b/imexam/imexam_defpars.py @@ -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"], diff --git a/imexam/imexamine.py b/imexam/imexamine.py index 0f1ae96..db8acdf 100644 --- a/imexam/imexamine.py +++ b/imexam/imexamine.py @@ -91,16 +91,17 @@ def __init__(self): self._define_default_pars() # default plot name saved with "s" key self.plot_name = "imexam_plot.pdf" - # let users have multiple plot windows, the list stores their names + # let users have multiple plot windows, the list stores their instance reference self._plot_windows = list() # this contains the name of the current plotting window self._figure_name = "imexam" - self._plot_windows.append(self._figure_name) + self._plot_windows.append(plt.gcf()) self._reserved_keys = ['q', '2'] # not to be changed with user funcs self._fit_models = ["Gaussian1D", "Moffat1D", "MexicanHat1D", "AiryDisk2D", + "Gaussian2D", "Polynomial1D"] # see if the package logger was already started @@ -129,7 +130,7 @@ def setlog(self, filename=None, on=True, level=logging.INFO): def _close_plots(self): """Make sure to release plot memory at end of exam loop.""" for plot in self._plot_windows: - plt.close() + plt.close(plot) def close(self): """For use with the Imexamine object standalone.""" @@ -295,8 +296,9 @@ def new_plot_window(self, x, y, data=None): """ if data is None: data = self._data + # save the old figure + self._plot_windows.append(plt.gcf()) self._figure_name = "imexam" + str(len(self._plot_windows) + 1) - self._plot_windows.append(self._figure_name) self.log.info(f"Plots now directed towards {self._figure_name}") def plot_line(self, x, y, data=None, fig=None): @@ -352,6 +354,7 @@ def plot_line(self, x, y, data=None, fig=None): plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() def plot_column(self, x, y, data=None, fig=None): """column plot of data at point y. @@ -408,6 +411,7 @@ def plot_column(self, x, y, data=None, fig=None): plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() def show_xy_coords(self, x, y, data=None): """print the x,y,value to the screen. @@ -478,17 +482,10 @@ def save_figure(self, fig=None): Parameters ---------- - data: numpy array - The data array to work on fig: figure for redirect Used for interaction with the ginga GUI """ - if fig is None: - fig = plt.figure(self._figure_name) - ax = fig.gca() - fig.savefig(self.plot_name) - pstr = f"plot saved to {self.plot_name}" - self.log.info(pstr) + self.save(fig=fig) def save(self, filename=None, fig=None): """Save to file the figure that's currently displayed. @@ -509,10 +506,10 @@ def save(self, filename=None, fig=None): self.set_plot_name(self._figure_name + ".pdf") if fig is None: - fig = plt.figure(self._figure_name) - ax = fig.gca() + fig = self._plot_windows[-1] + fig.savefig(self.plot_name) - pstr = f"plot saved to {self.plot_name}" + pstr = f"plot {fig.number} saved to {self.plot_name}" self.log.info(pstr) def aper_phot(self, x, y, data=None, @@ -633,7 +630,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" @@ -680,7 +677,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], @@ -697,6 +693,7 @@ def aper_phot(self, x, y, data=None, plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() else: return (apertures, annulus_apertures, rawflux_table, sky_per_pix) @@ -726,7 +723,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 @@ -880,6 +877,7 @@ def line_fit(self, x, y, data=None, form=None, genplot=True, fig=None, col=False plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() else: return fitted @@ -949,8 +947,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) @@ -1047,17 +1043,30 @@ 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 @@ -1067,41 +1076,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 + if pars["center"][0]: xx = int(x) yy = int(y) - # flipped from xpa 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) @@ -1157,7 +1155,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 " @@ -1168,16 +1166,16 @@ def radial_profile(self, x, y, data=None, form=None, 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" @@ -1187,10 +1185,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) @@ -1200,17 +1198,17 @@ 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 legendy = np.max(flux) / 2 if fitted is None: - msg = f"Problem with the {fitform.name} fit" + msg = f"Problem with the {form.name} fit" self.log.info(msg) raise ValueError(msg) @@ -1242,14 +1240,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) @@ -1260,6 +1258,7 @@ def radial_profile(self, x, y, data=None, form=None, plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() else: return radius, flux @@ -1356,6 +1355,7 @@ def curve_of_growth(self, x, y, data=None, genplot=True, fig=None): plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() else: return rapert, flux @@ -1525,6 +1525,7 @@ def histogram(self, x, y, data=None, genplot=True, fig=None): plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() else: hist, bin_edges = np.histogram(flat_data, num_bins, @@ -1598,6 +1599,7 @@ def contour(self, x, y, data=None, fig=None): plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() def surface(self, x, y, data=None, fig=None): """plot a surface around the specified location. @@ -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] @@ -1704,6 +1706,7 @@ def surface(self, x, y, data=None, fig=None): plt.pause(0.001) else: fig.canvas.draw_idle() + self._plot_windows[-1] = plt.gcf() def cutout(self, x, y, data=None, size=None, fig=None): """Make a fits cutout around the pointer location without wcs. diff --git a/imexam/math_helper.py b/imexam/math_helper.py index 826af02..cc8ac23 100644 --- a/imexam/math_helper.py +++ b/imexam/math_helper.py @@ -93,10 +93,12 @@ 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: @@ -115,7 +117,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) + @@ -165,12 +167,14 @@ 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: @@ -187,8 +191,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():