/usr/lib/python3/dist-packages/ginga/util/zscale.py is in python3-ginga 2.6.1-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 | """
This file is part of the STScI numdisplay package:
https://www.stsci.edu/trac/ssb/stsci_python/browser/stsci_python/trunk/numdisplay/lib/stsci/numdisplay/zscale.py?rev=13225
under the following license:
Copyright (C) 2005 Association of Universities for Research in Astronomy (AURA)
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following
disclaimer in the documentation and/or other materials provided
with the distribution.
3. The name of AURA and its representatives may not be used to
endorse or promote products derived from this software without
specific prior written permission.
THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR IMPLIED
WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
"""
import math
import numpy
MAX_REJECT = 0.5
MIN_NPIXELS = 5
GOOD_PIXEL = 0
BAD_PIXEL = 1
KREJ = 2.5
MAX_ITERATIONS = 5
def zscale(image, nsamples=1000, contrast=0.25):
"""Implement IRAF zscale algorithm
nsamples=1000 and contrast=0.25 are the IRAF display task defaults
image is a 2-d numpy array
returns (z1, z2)
"""
# Sample the image
samples = zsc_sample(image, nsamples)
return zscale_samples(samples, contrast=contrast)
def zsc_sample(image, maxpix, bpmask=None, zmask=None):
# Figure out which pixels to use for the zscale algorithm
# Returns the 1-d array samples
# Don't worry about the bad pixel mask or zmask for the moment
# Sample in a square grid, and return the first maxpix in the sample
nc = image.shape[0]
nl = image.shape[1]
stride = max(1.0, math.sqrt((nc - 1) * (nl - 1) / float(maxpix)))
stride = int(stride)
samples = image[::stride,::stride].flatten()
# remove NaN and Inf
samples = samples[numpy.isfinite(samples)]
return samples[:maxpix]
def zscale_samples(samples, contrast=0.25):
npix = len(samples)
samples.sort()
zmin = samples[0]
zmax = samples[-1]
# For a zero-indexed array
center_pixel = int((npix - 1) // 2)
if npix%2 == 1:
median = samples[center_pixel]
else:
median = 0.5 * (samples[center_pixel] + samples[center_pixel + 1])
#
# Fit a line to the sorted array of samples
minpix = max(MIN_NPIXELS, int(npix * MAX_REJECT))
ngrow = max(1, int (npix * 0.01))
ngoodpix, zstart, zslope = zsc_fit_line(samples, npix, KREJ, ngrow,
MAX_ITERATIONS)
#print "slope=%f intercept=%f" % (zslope, zstart)
if ngoodpix < minpix:
z1 = zmin
z2 = zmax
else:
if contrast > 0: zslope = zslope / contrast
z1 = max(zmin, median - (center_pixel - 1) * zslope)
z2 = min(zmax, median + (npix - center_pixel) * zslope)
return z1, z2
def zsc_fit_line(samples, npix, krej, ngrow, maxiter):
if npix <= 1:
return npix, 0, 1
#
# First re-map indices from -1.0 to 1.0
xscale = 2.0 / (npix - 1)
xnorm = numpy.arange(npix)
xnorm = xnorm * xscale - 1.0
ngoodpix = npix
minpix = max(MIN_NPIXELS, int(npix*MAX_REJECT))
last_ngoodpix = npix + 1
# This is the mask used in k-sigma clipping. 0 is good, 1 is bad
badpix = numpy.zeros(npix, dtype="int32")
#
# Iterate
for niter in range(maxiter):
if (ngoodpix >= last_ngoodpix) or (ngoodpix < minpix):
break
# Accumulate sums to calculate straight line fit
goodpixels = numpy.where(badpix == GOOD_PIXEL)
sumx = xnorm[goodpixels].sum()
sumxx = (xnorm[goodpixels]*xnorm[goodpixels]).sum()
sumxy = (xnorm[goodpixels]*samples[goodpixels]).sum()
sumy = samples[goodpixels].sum()
sum = len(goodpixels[0])
delta = sum * sumxx - sumx * sumx
# Slope and intercept
intercept = (sumxx * sumy - sumx * sumxy) / delta
slope = (sum * sumxy - sumx * sumy) / delta
# Subtract fitted line from the data array
fitted = xnorm*slope + intercept
flat = samples - fitted
# Compute the k-sigma rejection threshold
ngoodpix, mean, sigma = zsc_compute_sigma (flat, badpix, npix)
threshold = sigma * krej
# Detect and reject pixels further than k*sigma from the fitted line
lcut = -threshold
hcut = threshold
below = numpy.where(flat < lcut)
above = numpy.where(flat > hcut)
badpix[below] = BAD_PIXEL
badpix[above] = BAD_PIXEL
# Convolve with a kernel of length ngrow
kernel = numpy.ones(ngrow,dtype="int32")
badpix = numpy.convolve(badpix, kernel, mode='same')
ngoodpix = len(numpy.where(badpix == GOOD_PIXEL)[0])
niter += 1
# Transform the line coefficients back to the X range [0:npix-1]
zstart = intercept - slope
zslope = slope * xscale
return ngoodpix, zstart, zslope
def zsc_compute_sigma (flat, badpix, npix):
# Compute the rms deviation from the mean of a flattened array.
# Ignore rejected pixels
# Accumulate sum and sum of squares
goodpixels = numpy.where(badpix == GOOD_PIXEL)
sumz = flat[goodpixels].sum()
sumsq = (flat[goodpixels]*flat[goodpixels]).sum()
ngoodpix = len(goodpixels[0])
if ngoodpix == 0:
mean = None
sigma = None
elif ngoodpix == 1:
mean = sumz
sigma = None
else:
mean = sumz / ngoodpix
temp = sumsq / (ngoodpix - 1) - sumz*sumz / (ngoodpix * (ngoodpix - 1))
if temp < 0:
sigma = 0.0
else:
sigma = math.sqrt (temp)
return ngoodpix, mean, sigma
|