LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
rgb.py
Go to the documentation of this file.
1 import numpy as np
2 
3 from lsst.afw.display.displayLib import replaceSaturatedPixels
4 
5 class Mapping(object):
6  """!Baseclass to map red, blue, green intensities into uint8 values"""
7 
8  def __init__(self, min):
9  """!Create a mapping
10  \param min Intensity that should be mapped to black (a scalar or array for R, G, B)
11  """
12 
13  self._uint8Max = float(np.iinfo(np.uint8).max)
14 
15  try:
16  len(min)
17  except:
18  min = 3*[min]
19  assert len(min) == 3, "Please provide 1 or 3 values for min"
20 
21  self._min = min
22 
23  def makeRgbImage(self, imageR, imageG, imageB):
24  """!Convert 3 arrays, imageR, imageG, and imageB into a numpy RGB image
25 
26  N.b. images may be afwImages or numpy arrays
27  """
28  imageRGB = [imageR, imageG, imageB]
29  for i, c in enumerate(imageRGB):
30  if hasattr(c, "getImage"):
31  c = imageRGB[i] = c.getImage()
32  if hasattr(c, "getArray"):
33  imageRGB[i] = c.getArray()
34 
35  return np.flipud(np.dstack(self._convertImagesToUint8(*imageRGB)).astype(np.uint8))
36 
37  def intensity(self, imageR, imageG, imageB):
38  """!Return the total intensity from the red, blue, and green intensities"""
39  return (imageR + imageG + imageB)/float(3);
40 
41  def mapIntensityToUint8(self, I):
42  """Map an intensity into the range of a uint8, [0, 255] (but not converted to uint8)"""
43  return np.where(I <= 0, 0, np.where(I < self._uint8Max, I, self._uint8Max))
44 
45  def _convertImagesToUint8(self, imageR, imageG, imageB):
46  """Use the mapping to convert images imageR, imageG, and imageB to a triplet of uint8 images"""
47  imageR = imageR - self._min[0] # n.b. makes copy
48  imageG = imageG - self._min[1]
49  imageB = imageB - self._min[2]
50 
51  fac = self.mapIntensityToUint8(self.intensity(imageR, imageG, imageB))
52 
53  imageRGB = [imageR, imageG, imageB]
54 
55  for c in imageRGB:
56  c *= fac
57  c[c <= 0] = 0
58 
59  pixmax = self._uint8Max
60  r0, g0, b0 = imageRGB # copies -- could work row by row to minimise memory usage
61 
62  with np.errstate(invalid='ignore', divide='ignore'): # n.b. np.where doesn't (and can't) short-circuit
63  for i, c in enumerate(imageRGB):
64  c = np.where(r0 > g0,
65  np.where(r0 > b0,
66  np.where(r0 >= pixmax, c*pixmax/r0, c),
67  np.where(b0 >= pixmax, c*pixmax/b0, c)),
68  np.where(g0 > b0,
69  np.where(g0 >= pixmax, c*pixmax/g0, c),
70  np.where(b0 >= pixmax, c*pixmax/b0, c))).astype(np.uint8)
71  c[c > pixmax] = pixmax
72 
73  imageRGB[i] = c
74 
75  return imageRGB
76 
78  """!A mapping for an asinh stretch (preserving colours independent of brightness)
79 
80  x = asinh(Q (I - min)/range)/Q
81 
82  This reduces to a linear stretch if Q == 0
83 
84  See http://adsabs.harvard.edu/abs/2004PASP..116..133L
85  """
86 
87  def __init__(self, min, range, Q):
88  Mapping.__init__(self, min)
89 
90  epsilon = 1.0/2**23 # 32bit floating point machine epsilon; sys.float_info.epsilon is 64bit
91  if abs(Q) < epsilon:
92  Q = 0.1
93  else:
94  Qmax = 1e10
95  if Q > Qmax:
96  Q = Qmax
97 
98  if False:
99  self._slope = self._uint8Max/Q # gradient at origin is self._slope
100  else:
101  frac = 0.1 # gradient estimated using frac*range is _slope
102  self._slope = frac*self._uint8Max/np.arcsinh(frac*Q)
103 
104  self._soften = Q/float(range);
105 
106  def mapIntensityToUint8(self, I):
107  return np.where(I <= 0, 0, np.arcsinh(I*self._soften)*self._slope/I)
108 
109 def makeRGB(imageR, imageG, imageB, min=0, range=5, Q=20, fileName=None,
110  saturatedBorderWidth=0, saturatedPixelValue=None):
111  """Make a set of three images into an RGB image using an asinh stretch and optionally write it to disk
112 
113  If saturatedBorderWidth is non-zero, replace saturated pixels with saturatedPixelValue. Note
114  that replacing saturated pixels requires that the input images be MaskedImages.
115  """
116  if saturatedBorderWidth:
117  if saturatedPixelValue is None:
118  raise ValueError("saturatedPixelValue must be set if saturatedBorderWidth is set")
119  replaceSaturatedPixels(imageR, imageG, imageB, saturatedBorderWidth, saturatedPixelValue)
120 
121  asinhMap = AsinhMapping(min, range, Q)
122  rgb = asinhMap.makeRgbImage(imageR, imageG, imageB)
123  if fileName:
124  writeRGB(fileName, rgb)
125 
126  return rgb
127 
128 def displayRGB(rgb):
129  """Display an rgb image using matplotlib"""
130  import matplotlib.pyplot as plt
131  plt.imshow(rgb, interpolation='nearest')
132  plt.show()
133 
134 def writeRGB(fileName, rgbImage):
135  """Write an RGB image (made by e.g. makeRGB) to fileName"""
136  import matplotlib.image
137  matplotlib.image.imsave(fileName, rgbImage)
138 
139 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
140 #
141 # Support the legacy API
142 #
143 class asinhMappingF(object):
144  def __init__(self, min, range, Q):
145  self.min = min
146  self.range = range
147  self.Q = Q
148 
149 class _RgbImageF(object):
150  def __init__(self, imageR, imageG, imageB, mapping):
151  asinh = AsinhMapping(mapping.min, mapping.range, mapping.Q)
152  self.rgb = asinh.makeRgbImage(imageR, imageG, imageB)
153 
154  def write(self, fileName):
155  writeRGB(fileName, self.rgb)
156 
157 def RgbImageF(imageR, imageG, imageB, mapping):
158  return _RgbImageF(imageR, imageG, imageB, mapping)
def intensity
Return the total intensity from the red, blue, and green intensities.
Definition: rgb.py:37
Baseclass to map red, blue, green intensities into uint8 values.
Definition: rgb.py:5
def makeRgbImage
Convert 3 arrays, imageR, imageG, and imageB into a numpy RGB image.
Definition: rgb.py:23
def __init__
Create a mapping.
Definition: rgb.py:8
A mapping for an asinh stretch (preserving colours independent of brightness)
Definition: rgb.py:77
template void replaceSaturatedPixels(image::MaskedImage< float > &rim, image::MaskedImage< float > &gim, image::MaskedImage< float > &bim, int borderWidth, float saturatedPixelValue)