LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
_rgbContinued.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2015-2016 LSST/AURA
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <https://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 import numpy as np
24 
25 import lsst.afw.image as afwImage
26 import lsst.afw.math as afwMath
27 from ._rgb import replaceSaturatedPixels, getZScale
28 
29 
30 def computeIntensity(imageR, imageG=None, imageB=None):
31  """Return a naive total intensity from the red, blue, and green intensities
32 
33  Parameters
34  ----------
35  imageR : `lsst.afw.image.MaskedImage`, `lsst.afw.image.Image`, or `numpy.ndarray`, (Nx, Ny)
36  intensity of image that'll be mapped to red; or intensity if imageG and imageB are None
37  imageG : `lsst.afw.image.MaskedImage`, `lsst.afw.image.Image`, or `numpy.ndarray`, (Nx, Ny)
38  intensity of image that'll be mapped to green; or None
39  imageB : `lsst.afw.image.MaskedImage`, `lsst.afw.image.Image`, or `numpy.ndarray`, (Nx, Ny)
40  intensity of image that'll be mapped to blue; or None
41 
42  Returns
43  -------
44  image : type of ``imageR``, ``imageG``, and `imageB``
45  """
46  if imageG is None or imageB is None:
47  assert imageG is None and imageB is None, \
48  "Please specify either a single image or red, green, and blue images"
49  return imageR
50 
51  imageRGB = [imageR, imageG, imageB]
52 
53  for i, c in enumerate(imageRGB):
54  if hasattr(c, "getImage"):
55  c = imageRGB[i] = c.getImage()
56  if hasattr(c, "getArray"):
57  imageRGB[i] = c.getArray()
58 
59  intensity = (imageRGB[0] + imageRGB[1] + imageRGB[2])/float(3)
60  #
61  # Repack into whatever type was passed to us
62  #
63  Image = afwImage.ImageU if intensity.dtype == 'uint16' else afwImage.ImageF
64 
65  if hasattr(imageR, "getImage"): # a maskedImage
66  intensity = afwImage.makeMaskedImage(Image(intensity))
67  elif hasattr(imageR, "getArray"):
68  intensity = Image(intensity)
69 
70  return intensity
71 
72 
73 class Mapping:
74  """Base class to map red, blue, green intensities into uint8 values
75 
76  Parameters
77  ----------
78  minimum : `float` or sequence of `float`
79  Intensity that should be mapped to black. If an array, has three
80  elements for R, G, B.
81  image
82  The image to be used to calculate the mapping.
83  If provided, also the default for makeRgbImage()
84  """
85 
86  def __init__(self, minimum=None, image=None):
87  self._uint8Max_uint8Max = float(np.iinfo(np.uint8).max)
88 
89  try:
90  len(minimum)
91  except TypeError:
92  minimum = 3*[minimum]
93  assert len(minimum) == 3, "Please provide 1 or 3 values for minimum"
94 
95  self.minimumminimum = minimum
96  self._image_image = image
97 
98  def makeRgbImage(self, imageR=None, imageG=None, imageB=None,
99  xSize=None, ySize=None, rescaleFactor=None):
100  """Convert 3 arrays, imageR, imageG, and imageB into a numpy RGB image
101 
102  imageR : `lsst.afw.image.Image` or `numpy.ndarray`, (Nx, Ny)
103  Image to map to red (if `None`, use the image passed to the ctor)
104  imageG : `lsst.afw.image.Image` or `numpy.ndarray`, (Nx, Ny), optional
105  Image to map to green (if `None`, use imageR)
106  imageB : `lsst.afw.image.Image` or `numpy.ndarray`, (Nx, Ny), optional
107  Image to map to blue (if `None`, use imageR)
108  xSize : `int`, optional
109  Desired width of RGB image. If ``ySize`` is `None`, preserve aspect ratio
110  ySize : `int`, optional
111  Desired height of RGB image
112  rescaleFactor : `float`, optional
113  Make size of output image ``rescaleFactor*size`` of the input image
114  """
115  if imageR is None:
116  if self._image_image is None:
117  raise RuntimeError(
118  "You must provide an image (or pass one to the constructor)")
119  imageR = self._image_image
120 
121  if imageG is None:
122  imageG = imageR
123  if imageB is None:
124  imageB = imageR
125 
126  imageRGB = [imageR, imageG, imageB]
127  for i, c in enumerate(imageRGB):
128  if hasattr(c, "getImage"):
129  c = imageRGB[i] = c.getImage()
130  if hasattr(c, "getArray"):
131  imageRGB[i] = c.getArray()
132 
133  if xSize is not None or ySize is not None:
134  assert rescaleFactor is None, "You may not specify a size and rescaleFactor"
135  h, w = imageRGB[0].shape
136  if ySize is None:
137  ySize = int(xSize*h/float(w) + 0.5)
138  elif xSize is None:
139  xSize = int(ySize*w/float(h) + 0.5)
140 
141  size = (ySize, xSize) # n.b. y, x order for scipy
142  elif rescaleFactor is not None:
143  size = float(rescaleFactor) # an int is intepreted as a percentage
144  else:
145  size = None
146 
147  if size is not None:
148  try:
149  import scipy.misc
150  except ImportError as e:
151  raise RuntimeError(
152  f"Unable to rescale as scipy.misc is unavailable: {e}")
153 
154  for i, im in enumerate(imageRGB):
155  imageRGB[i] = scipy.misc.imresize(
156  im, size, interp='bilinear', mode='F')
157 
158  return np.dstack(self._convertImagesToUint8_convertImagesToUint8(*imageRGB)).astype(np.uint8)
159 
160  def intensity(self, imageR, imageG, imageB):
161  """Return the total intensity from the red, blue, and green intensities
162 
163  Notes
164  -----
165  This is a naive computation, and may be overridden by subclasses
166  """
167  return computeIntensity(imageR, imageG, imageB)
168 
169  def mapIntensityToUint8(self, intensity):
170  """Map an intensity into the range of a uint8, [0, 255] (but not converted to uint8)
171  """
172  with np.errstate(invalid='ignore', divide='ignore'): # n.b. np.where can't and doesn't short-circuit
173  return np.where(intensity <= 0, 0,
174  np.where(intensity < self._uint8Max_uint8Max, intensity, self._uint8Max_uint8Max))
175 
176  def _convertImagesToUint8(self, imageR, imageG, imageB):
177  """Use the mapping to convert images imageR, imageG, and imageB to a triplet of uint8 images
178  """
179  imageR = imageR - self.minimumminimum[0] # n.b. makes copy
180  imageG = imageG - self.minimumminimum[1]
181  imageB = imageB - self.minimumminimum[2]
182 
183  fac = self.mapIntensityToUint8mapIntensityToUint8(self.intensityintensity(imageR, imageG, imageB))
184 
185  imageRGB = [imageR, imageG, imageB]
186  with np.errstate(invalid="ignore"): # suppress NAN warnings
187  for c in imageRGB:
188  c *= fac
189  # individual bands can still be < 0, even if fac isn't
190  c[c < 0] = 0
191 
192  pixmax = self._uint8Max_uint8Max
193  # copies -- could work row by row to minimise memory usage
194  r0, g0, b0 = imageRGB
195 
196  # n.b. np.where can't and doesn't short-circuit
197  with np.errstate(invalid='ignore', divide='ignore'):
198  for i, c in enumerate(imageRGB):
199  c = np.where(r0 > g0,
200  np.where(r0 > b0,
201  np.where(r0 >= pixmax, c*pixmax/r0, c),
202  np.where(b0 >= pixmax, c*pixmax/b0, c)),
203  np.where(g0 > b0,
204  np.where(g0 >= pixmax, c*pixmax/g0, c),
205  np.where(b0 >= pixmax, c*pixmax/b0, c))).astype(np.uint8)
206  c[c > pixmax] = pixmax
207 
208  imageRGB[i] = c
209 
210  return imageRGB
211 
212 
214  """A linear map of red, blue, green intensities into uint8 values
215 
216  Parameters
217  ----------
218  minimum : `float` or sequence of `float`
219  Intensity that should be mapped to black. If an array, has three
220  elements for R, G, B.
221  maximum : `float`
222  Intensity that should be mapped to white
223  image
224  Image to estimate minimum/maximum if not explicitly set
225  """
226 
227  def __init__(self, minimum=None, maximum=None, image=None):
228  if minimum is None or maximum is None:
229  assert image is not None, "You must provide an image if you don't set both minimum and maximum"
230 
231  stats = afwMath.makeStatistics(image, afwMath.MIN | afwMath.MAX)
232  if minimum is None:
233  minimum = stats.getValue(afwMath.MIN)
234  if maximum is None:
235  maximum = stats.getValue(afwMath.MAX)
236 
237  Mapping.__init__(self, minimum, image)
238  self.maximummaximum = maximum
239 
240  if maximum is None:
241  self._range_range = None
242  else:
243  assert maximum - minimum != 0, "minimum and maximum values must not be equal"
244  self._range_range = float(maximum - minimum)
245 
246  def mapIntensityToUint8(self, intensity):
247  """Return an array which, when multiplied by an image, returns that
248  image mapped to the range of a uint8, [0, 255] (but not converted to uint8)
249 
250  The intensity is assumed to have had ``minimum`` subtracted (as that
251  can be done per-band)
252  """
253  with np.errstate(invalid='ignore', divide='ignore'): # n.b. np.where can't and doesn't short-circuit
254  return np.where(intensity <= 0, 0,
255  np.where(intensity >= self._range_range,
256  self._uint8Max_uint8Max/intensity, self._uint8Max_uint8Max/self._range_range))
257 
258 
260  """A mapping for a linear stretch chosen by the zscale algorithm
261  (preserving colours independent of brightness)
262 
263  x = (I - minimum)/range
264 
265  Parameters
266  ----------
267  image
268  Image whose parameters are desired
269  nSamples : `int`
270  The number of samples to use to estimate the zscale parameters
271  contrast : `float`
272  """
273 
274  def __init__(self, image, nSamples=1000, contrast=0.25):
275  if not hasattr(image, "getArray"):
276  image = afwImage.ImageF(image)
277  z1, z2 = getZScale(image, nSamples, contrast)
278 
279  LinearMapping.__init__(self, z1, z2, image)
280 
281 
283  """A mapping for an asinh stretch (preserving colours independent of brightness)
284 
285  x = asinh(Q (I - minimum)/range)/Q
286 
287  Notes
288  -----
289  This reduces to a linear stretch if Q == 0
290 
291  See http://adsabs.harvard.edu/abs/2004PASP..116..133L
292  """
293 
294  def __init__(self, minimum, dataRange, Q=8):
295  Mapping.__init__(self, minimum)
296 
297  # 32bit floating point machine epsilon; sys.float_info.epsilon is 64bit
298  epsilon = 1.0/2**23
299  if abs(Q) < epsilon:
300  Q = 0.1
301  else:
302  Qmax = 1e10
303  if Q > Qmax:
304  Q = Qmax
305 
306  if False:
307  self._slope_slope = self._uint8Max_uint8Max/Q # gradient at origin is self._slope
308  else:
309  frac = 0.1 # gradient estimated using frac*range is _slope
310  self._slope_slope = frac*self._uint8Max_uint8Max/np.arcsinh(frac*Q)
311 
312  self._soften_soften = Q/float(dataRange)
313 
314  def mapIntensityToUint8(self, intensity):
315  """Return an array which, when multiplied by an image, returns that image mapped to the range of a
316  uint8, [0, 255] (but not converted to uint8)
317 
318  The intensity is assumed to have had minimum subtracted (as that can be done per-band)
319  """
320  with np.errstate(invalid='ignore', divide='ignore'): # n.b. np.where can't and doesn't short-circuit
321  return np.where(intensity <= 0, 0, np.arcsinh(intensity*self._soften_soften)*self._slope_slope/intensity)
322 
323 
325  """A mapping for an asinh stretch, estimating the linear stretch by zscale
326 
327  x = asinh(Q (I - z1)/(z2 - z1))/Q
328 
329  Parameters
330  ----------
331  image
332  The image to analyse, or a list of 3 images to be converted to an intensity image
333  Q : `int`
334  The asinh softening parameter
335  pedestal : `float` or sequence of `float`, optional
336  The value, or array of 3 values, to subtract from the images
337 
338  N.b. pedestal, if not None, is removed from the images when calculating the zscale
339  stretch, and added back into Mapping.minimum[]
340 
341  See also
342  --------
343  AsinhMapping
344  """
345 
346  def __init__(self, image, Q=8, pedestal=None):
347  try:
348  assert len(image) in (1, 3,), "Please provide 1 or 3 images"
349  except TypeError:
350  image = [image]
351 
352  if pedestal is not None:
353  try:
354  assert len(pedestal) in (
355  1, 3,), "Please provide 1 or 3 pedestals"
356  except TypeError:
357  pedestal = 3*[pedestal]
358 
359  image = list(image) # needs to be mutable
360  for i, im in enumerate(image):
361  if pedestal[i] != 0.0:
362  if hasattr(im, "getImage"):
363  im = im.getImage()
364  if hasattr(im, "getArray"):
365  im = im.getArray()
366 
367  image[i] = im - pedestal[i] # n.b. a copy
368  else:
369  pedestal = len(image)*[0.0]
370 
371  image = computeIntensity(*image)
372 
373  zscale = ZScaleMapping(image)
374  # zscale.minimum is always a triple
375  dataRange = zscale.maximum - zscale.minimum[0]
376  minimum = zscale.minimum
377 
378  for i, level in enumerate(pedestal):
379  minimum[i] += level
380 
381  AsinhMapping.__init__(self, minimum, dataRange, Q)
382  self._image_image_image = image # support self.makeRgbImage()
383 
384 
385 def makeRGB(imageR, imageG=None, imageB=None, minimum=0, dataRange=5, Q=8, fileName=None,
386  saturatedBorderWidth=0, saturatedPixelValue=None,
387  xSize=None, ySize=None, rescaleFactor=None):
388  """Make a set of three images into an RGB image using an asinh stretch and
389  optionally write it to disk
390 
391  Parameters
392  ----------
393  imageR
394  imageG
395  imageB
396  minimum : `float` or sequence of `float`
397  dataRange
398  Q : `int`
399  fileName : `str`
400  The output file. The suffix defines the format, and must be supported by matplotlib
401  saturatedBorderWidth
402  If saturatedBorderWidth is non-zero, replace saturated pixels with
403  ``saturatedPixelValue``. Note that replacing saturated pixels requires
404  that the input images be `lsst.afw.image.MaskedImage`.
405  saturatedPixelValue
406  xSize
407  ySize
408  rescaleFactor
409  """
410  if imageG is None:
411  imageG = imageR
412  if imageB is None:
413  imageB = imageR
414 
415  if saturatedBorderWidth:
416  if saturatedPixelValue is None:
417  raise ValueError(
418  "saturatedPixelValue must be set if saturatedBorderWidth is set")
419  replaceSaturatedPixels(imageR, imageG, imageB,
420  saturatedBorderWidth, saturatedPixelValue)
421 
422  asinhMap = AsinhMapping(minimum, dataRange, Q)
423  rgb = asinhMap.makeRgbImage(imageR, imageG, imageB,
424  xSize=xSize, ySize=ySize, rescaleFactor=rescaleFactor)
425 
426  if fileName:
427  writeRGB(fileName, rgb)
428 
429  return rgb
430 
431 
432 def displayRGB(rgb, show=True):
433  """Display an rgb image using matplotlib
434 
435  Parameters
436  ----------
437  rgb
438  The RGB image in question
439  show : `bool`
440  If `True`, call `matplotlib.pyplot.show()`
441  """
442  import matplotlib.pyplot as plt
443  plt.imshow(rgb, interpolation='nearest', origin="lower")
444  if show:
445  plt.show()
446  return plt
447 
448 
449 def writeRGB(fileName, rgbImage):
450  """Write an RGB image to disk
451 
452  Parameters
453  ----------
454  fileName : `str`
455  The output file. The suffix defines the format, and must be supported by matplotlib
456 
457  Most versions of matplotlib support png and pdf (although the eps/pdf/svg writers may be buggy,
458  possibly due an interaction with useTeX=True in the matplotlib settings).
459 
460  If your matplotlib bundles pil/pillow you should also be able to write jpeg and tiff files.
461  rgbImage
462  The image, as made by e.g. makeRGB
463  """
464  import matplotlib.image
465  matplotlib.image.imsave(fileName, rgbImage)
def __init__(self, minimum, dataRange, Q=8)
def __init__(self, image, Q=8, pedestal=None)
def __init__(self, minimum=None, maximum=None, image=None)
def makeRgbImage(self, imageR=None, imageG=None, imageB=None, xSize=None, ySize=None, rescaleFactor=None)
def _convertImagesToUint8(self, imageR, imageG, imageB)
def intensity(self, imageR, imageG, imageB)
def __init__(self, minimum=None, image=None)
def __init__(self, image, nSamples=1000, contrast=0.25)
daf::base::PropertyList * list
Definition: fits.cc:913
def computeIntensity(imageR, imageG=None, imageB=None)
def writeRGB(fileName, rgbImage)
def makeRGB(imageR, imageG=None, imageB=None, minimum=0, dataRange=5, Q=8, fileName=None, saturatedBorderWidth=0, saturatedPixelValue=None, xSize=None, ySize=None, rescaleFactor=None)
void replaceSaturatedPixels(ImageT &rim, ImageT &gim, ImageT &bim, int borderWidth, float saturatedPixelValue)
Definition: _saturated.cc:32
std::pair< double, double > getZScale(image::Image< T > const &image, int const nSamples, double const contrast)
Calculate an IRAF/ds9-style zscaling.
Definition: _scaling.cc:167
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT >> image, typename std::shared_ptr< Mask< MaskPixelT >> mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT >> variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
Definition: MaskedImage.h:1268
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:354
Angle abs(Angle const &a)
Definition: Angle.h:106