From cd21c7acfe29a40d1ff170f93fdd87b0540151af Mon Sep 17 00:00:00 2001 From: Megan Sosey Date: Mon, 18 Oct 2021 15:59:51 -0400 Subject: [PATCH 1/3] fixes for #165, #238, #210, #205 --- CHANGES.rst | 4 +- imexam/imexam_defpars.py | 1 + imexam/imexamine.py | 94 ++++++++++++++++++++-------------------- imexam/math_helper.py | 28 +++++++----- 4 files changed, 70 insertions(+), 57 deletions(-) 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/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..472671f 100644 --- a/imexam/imexamine.py +++ b/imexam/imexamine.py @@ -101,6 +101,7 @@ def __init__(self): "Moffat1D", "MexicanHat1D", "AiryDisk2D", + "Gaussian2D", "Polynomial1D"] # see if the package logger was already started @@ -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" @@ -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], @@ -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 @@ -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) @@ -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 @@ -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) @@ -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 " @@ -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" @@ -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) @@ -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 @@ -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) @@ -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] diff --git a/imexam/math_helper.py b/imexam/math_helper.py index 826af02..05f0c8c 100644 --- a/imexam/math_helper.py +++ b/imexam/math_helper.py @@ -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: @@ -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) + @@ -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: @@ -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(): From f9cc6ca8ee01584168ff5bd2a5574cca1ec8f700 Mon Sep 17 00:00:00 2001 From: Megan Sosey Date: Mon, 18 Oct 2021 18:24:07 -0400 Subject: [PATCH 2/3] fix plot saving bug in imexam loop --- imexam/connect.py | 2 +- imexam/imexamine.py | 38 +++++++++++++++++++++----------------- 2 files changed, 22 insertions(+), 18 deletions(-) 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/imexamine.py b/imexam/imexamine.py index 472671f..1e2d95a 100644 --- a/imexam/imexamine.py +++ b/imexam/imexamine.py @@ -91,11 +91,11 @@ 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", @@ -130,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.""" @@ -296,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): @@ -353,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. @@ -409,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. @@ -479,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. @@ -510,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.savefig(self.plot_name) - pstr = f"plot saved to {self.plot_name}" + fig = self._plot_windows[-1] + + fig.savefig(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, @@ -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) @@ -880,6 +877,8 @@ 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 @@ -1044,7 +1043,7 @@ def radial_profile(self, x, y, data=None, form=None, """ pars = self.radial_profile_pars - + if data is None: data = self._data @@ -1262,6 +1261,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 @@ -1358,6 +1358,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 @@ -1527,6 +1528,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, @@ -1600,6 +1602,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. @@ -1706,6 +1709,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. From 9fe3d3397f23c10e9e0204791914fd61b197e12d Mon Sep 17 00:00:00 2001 From: Megan Sosey Date: Mon, 18 Oct 2021 18:36:16 -0400 Subject: [PATCH 3/3] pep8 fixes --- imexam/imexamine.py | 17 +++++++---------- imexam/math_helper.py | 10 ++++------ 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/imexam/imexamine.py b/imexam/imexamine.py index 1e2d95a..db8acdf 100644 --- a/imexam/imexamine.py +++ b/imexam/imexamine.py @@ -507,8 +507,8 @@ def save(self, filename=None, fig=None): if fig is None: fig = self._plot_windows[-1] - - fig.savefig(self.plot_name) + + fig.savefig(self.plot_name) pstr = f"plot {fig.number} saved to {self.plot_name}" self.log.info(pstr) @@ -879,7 +879,6 @@ def line_fit(self, x, y, data=None, form=None, genplot=True, fig=None, col=False fig.canvas.draw_idle() self._plot_windows[-1] = plt.gcf() - else: return fitted @@ -1053,7 +1052,7 @@ def radial_profile(self, x, y, data=None, form=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]) @@ -1061,7 +1060,6 @@ def radial_profile(self, x, y, data=None, form=None, self.log.warning("Install photutils to enable " "background subtraction") subtract_background = False - getdata = bool(pars["getdata"][0]) delta = int(pars["delta"][0]) @@ -1085,8 +1083,8 @@ def radial_profile(self, x, y, data=None, form=None, datasize = 3 if pars["center"][0]: - xx=int(x) - yy=int(y) + xx = int(x) + yy = int(y) data_chunk = data[yy - datasize:yy + datasize, xx - datasize:xx + datasize] amp, xx, yy, sigma, sigmay = self.gauss_center(xx, @@ -1167,12 +1165,11 @@ 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 if fitplot: fline = np.linspace(0, datasize, 100) # finer sample # fit model to data - if form.name == "Gaussian1D":# make center radii at max ******** + if form.name == "Gaussian1D": # make center radii at max ******** fitted = math_helper.fit_gauss_1d(radius, flux, sigma_factor=sig_factor, @@ -1211,7 +1208,7 @@ def radial_profile(self, x, y, data=None, form=None, 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) diff --git a/imexam/math_helper.py b/imexam/math_helper.py index 05f0c8c..cc8ac23 100644 --- a/imexam/math_helper.py +++ b/imexam/math_helper.py @@ -100,7 +100,6 @@ def fit_moffat_1d(data, gamma=2., alpha=1., sigma_factor=0., x0 = center_at fluxmax = max(data) - # assumes negligable background if weighted: z = np.nan_to_num(1. / np.sqrt(data)) # use as weight @@ -168,15 +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") - if center_at is None: - center_mean = sum(radius)/len(radius) + center_mean = sum(radius) / len(radius) fixed = {'amplitude': False} - bounds = {'amplitude': (flux.max()/2.,flux.max())} + bounds = {'amplitude': (flux.max() / 2., flux.max())} else: center_mean = center_at fixed = {'mean': True, 'amplitude': True} - bounds = {'amplitude': (flux.max()/2., flux.max())} + bounds = {'amplitude': (flux.max() / 2., flux.max())} # assumes negligable background if weighted: @@ -194,7 +192,7 @@ def fit_gauss_1d(radius, flux, sigma_factor=0, center_at=None, weighted=False): # Gaussian1D + a constant model = (models.Gaussian1D(amplitude=flux.max(), - mean=center_mean, stddev=1., + mean=center_mean, stddev=1., fixed=fixed, bounds=bounds) + models.Polynomial1D(c0=flux.min(), degree=0))