From 91e53ec7c6da162ec61289749838e62e2f2ea31b Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Mon, 28 Nov 2016 14:48:44 -0500 Subject: [PATCH 1/3] Move zscale code into ZScaleInterval --- astropy/visualization/interval.py | 91 +++++++++++++++++++++++++------ 1 file changed, 74 insertions(+), 17 deletions(-) diff --git a/astropy/visualization/interval.py b/astropy/visualization/interval.py index 34b478c8ae9..955c7062c76 100644 --- a/astropy/visualization/interval.py +++ b/astropy/visualization/interval.py @@ -13,7 +13,6 @@ from ..extern import six from ..utils.misc import InheritDocstrings from .transform import BaseTransform -from .zscale import zscale __all__ = ['BaseInterval', 'ManualInterval', 'MinMaxInterval', @@ -189,28 +188,32 @@ class ZScaleInterval(BaseInterval): http://iraf.net/forum/viewtopic.php?showtopic=134139 + Original implementation: + https://trac.stsci.edu/ssb/stsci_python/browser/stsci_python/trunk/numdisplay/lib/stsci/numdisplay/zscale.py?rev=19347 + Parameters ---------- - image : array_like - Input array. nsamples : int, optional - Number of points in array to sample for determining scaling - factors. Defaults to 1000. + The number of points in the array to sample for determining + scaling factors. Defaults to 1000. contrast : float, optional - Scaling factor (between 0 and 1) for determining min and max. - Larger values increase the difference between min and max values - used for display. Defaults to 0.25. + The scaling factor (between 0 and 1) for determining the minimum + and maximum value. Larger values increase the difference + between the minimum and maximum values used for display. + Defaults to 0.25. max_reject : float, optional If more than ``max_reject * npixels`` pixels are rejected, then - the returned values are the min and max of the data. Defaults to - 0.5. + the returned values are the minimum and maximum of the data. + Defaults to 0.5. min_npixels : int, optional If less than ``min_npixels`` pixels are rejected, then the - returned values are the min and max of the data. Defaults to 5. + returned values are the minimum and maximum of the data. + Defaults to 5. krej : float, optional - Number of sigma used for the rejection. Defaults to 2.5. + The number of sigma used for the rejection. Defaults to 2.5. max_iterations : int, optional - Maximum number of iterations for the rejection. Defaults to 5. + The maximum number of iterations for the rejection. Defaults to + 5. """ def __init__(self, nsamples=1000, contrast=0.25, max_reject=0.5, @@ -223,7 +226,61 @@ def __init__(self, nsamples=1000, contrast=0.25, max_reject=0.5, self.max_iterations = max_iterations def get_limits(self, values): - return zscale(values, nsamples=self.nsamples, contrast=self.contrast, - max_reject=self.max_reject, - min_npixels=self.min_npixels, krej=self.krej, - max_iterations=self.max_iterations) + # Sample the image + values = np.asarray(values) + values = values[np.isfinite(values)] + stride = int(max(1.0, values.size / self.nsamples)) + samples = values[::stride][:self.nsamples] + samples.sort() + + npix = len(samples) + vmin = samples[0] + vmax = samples[-1] + + # Fit a line to the sorted array of samples + minpix = max(self.min_npixels, int(npix * self.max_reject)) + x = np.arange(npix) + ngoodpix = npix + last_ngoodpix = npix + 1 + + # Bad pixels mask used in k-sigma clipping + badpix = np.zeros(npix, dtype=bool) + + # Kernel used to dilate the bad pixels mask + ngrow = max(1, int(npix * 0.01)) + kernel = np.ones(ngrow, dtype=bool) + + for niter in six.moves.range(self.max_iterations): + if ngoodpix >= last_ngoodpix or ngoodpix < minpix: + break + + fit = np.polyfit(x, samples, deg=1, w=(~badpix).astype(int)) + fitted = np.poly1d(fit)(x) + + # Subtract fitted line from the data array + flat = samples - fitted + + # Compute the k-sigma rejection threshold + threshold = self.krej * flat[~badpix].std() + + # Detect and reject pixels further than k*sigma from the + # fitted line + badpix[(flat < - threshold) | (flat > threshold)] = True + + # Convolve with a kernel of length ngrow + badpix = np.convolve(badpix, kernel, mode='same') + + last_ngoodpix = ngoodpix + ngoodpix = np.sum(~badpix) + + slope, intercept = fit + + if ngoodpix >= minpix: + if self.contrast > 0: + slope = slope / self.contrast + center_pixel = (npix - 1) // 2 + median = np.median(samples) + vmin = max(vmin, median - (center_pixel - 1) * slope) + vmax = min(vmax, median + (npix - center_pixel) * slope) + + return vmin, vmax From 50f1c8270ef7fca47a04f2a9d494298493dc5b28 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Mon, 28 Nov 2016 14:52:01 -0500 Subject: [PATCH 2/3] Remove zscale.py module --- astropy/visualization/zscale.py | 111 -------------------------------- 1 file changed, 111 deletions(-) delete mode 100644 astropy/visualization/zscale.py diff --git a/astropy/visualization/zscale.py b/astropy/visualization/zscale.py deleted file mode 100644 index 4979583e9eb..00000000000 --- a/astropy/visualization/zscale.py +++ /dev/null @@ -1,111 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst - -# Adapted from code from the Numdisplay package (see NUMDISPLAY_LICENSE.rst) - -""" -Zscale implementation based on the one from the STScI Numdisplay package. - -Original implementation from Numdisplay is BSD licensed: -https://trac.stsci.edu/ssb/stsci_python/browser/stsci_python/trunk/numdisplay/lib/stsci/numdisplay/zscale.py?rev=19347 -https://trac.stsci.edu/ssb/stsci_python/browser/stsci_python/trunk/numdisplay/LICENSE.txt?rev=19347 -""" - -from __future__ import absolute_import, division - -import numpy as np - -from ..extern.six.moves import range - - -__all__ = ['zscale'] - - -def zscale(image, nsamples=1000, contrast=0.25, max_reject=0.5, min_npixels=5, - krej=2.5, max_iterations=5): - """Implement IRAF zscale algorithm. - - Parameters - ---------- - image : array_like - Input array. - nsamples : int, optional - Number of points in array to sample for determining scaling factors. - Defaults to 1000. - contrast : float, optional - Scaling factor (between 0 and 1) for determining min and max. Larger - values increase the difference between min and max values used for - display. Defaults to 0.25. - max_reject : float, optional - If more than ``max_reject * npixels`` pixels are rejected, then the - returned values are the min and max of the data. Defaults to 0.5. - min_npixels : int, optional - If less than ``min_npixels`` pixels are rejected, then the - returned values are the min and max of the data. Defaults to 5. - krej : float, optional - Number of sigma used for the rejection. Defaults to 2.5. - max_iterations : int, optional - Maximum number of iterations for the rejection. Defaults to 5. - - Returns - ------- - zmin, zmax: float - Computed min and max values. - - """ - - # Sample the image - image = np.asarray(image) - image = image[np.isfinite(image)] - stride = int(max(1.0, image.size / nsamples)) - samples = image[::stride][:nsamples] - samples.sort() - - npix = len(samples) - zmin = samples[0] - zmax = samples[-1] - - # Fit a line to the sorted array of samples - minpix = max(min_npixels, int(npix * max_reject)) - x = np.arange(npix) - ngoodpix = npix - last_ngoodpix = npix + 1 - - # Bad pixels mask used in k-sigma clipping - badpix = np.zeros(npix, dtype=bool) - - # Kernel used to dilate the bad pixels mask - ngrow = max(1, int(npix * 0.01)) - kernel = np.ones(ngrow, dtype=bool) - - for niter in range(max_iterations): - if ngoodpix >= last_ngoodpix or ngoodpix < minpix: - break - - fit = np.polyfit(x, samples, deg=1, w=(~badpix).astype(int)) - fitted = np.poly1d(fit)(x) - - # Subtract fitted line from the data array - flat = samples - fitted - - # Compute the k-sigma rejection threshold - threshold = krej * flat[~badpix].std() - - # Detect and reject pixels further than k*sigma from the fitted line - badpix[(flat < - threshold) | (flat > threshold)] = True - - # Convolve with a kernel of length ngrow - badpix = np.convolve(badpix, kernel, mode='same') - - last_ngoodpix = ngoodpix - ngoodpix = np.sum(~badpix) - - slope, intercept = fit - - if ngoodpix >= minpix: - if contrast > 0: - slope = slope / contrast - center_pixel = (npix - 1) // 2 - median = np.median(samples) - zmin = max(zmin, median - (center_pixel - 1) * slope) - zmax = min(zmax, median + (npix - center_pixel) * slope) - return zmin, zmax From d1c2937c8234c348d94e6f3b0ee99d521b249224 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Mon, 28 Nov 2016 15:50:33 -0500 Subject: [PATCH 3/3] Renamed license to AURA_LICENSE; added note about license in ZScaleInterval --- astropy/visualization/interval.py | 2 ++ licenses/{NUMDISPLAY_LICENSE.rst => AURA_LICENSE.rst} | 0 2 files changed, 2 insertions(+) rename licenses/{NUMDISPLAY_LICENSE.rst => AURA_LICENSE.rst} (100%) diff --git a/astropy/visualization/interval.py b/astropy/visualization/interval.py index 955c7062c76..f1249fbda8b 100644 --- a/astropy/visualization/interval.py +++ b/astropy/visualization/interval.py @@ -191,6 +191,8 @@ class ZScaleInterval(BaseInterval): Original implementation: https://trac.stsci.edu/ssb/stsci_python/browser/stsci_python/trunk/numdisplay/lib/stsci/numdisplay/zscale.py?rev=19347 + Licensed under a 3-clause BSD style license (see AURA_LICENSE.rst). + Parameters ---------- nsamples : int, optional diff --git a/licenses/NUMDISPLAY_LICENSE.rst b/licenses/AURA_LICENSE.rst similarity index 100% rename from licenses/NUMDISPLAY_LICENSE.rst rename to licenses/AURA_LICENSE.rst