LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
_rgbContinued.py
Go to the documentation of this file.
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
23import numpy as np
24
25import lsst.afw.image as afwImage
26import lsst.afw.math as afwMath
27from ._rgb import replaceSaturatedPixels, getZScale
28
29
30def 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
73class 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 = 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.minimum = minimum
96 self._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 is None:
117 raise RuntimeError(
118 "You must provide an image (or pass one to the constructor)")
119 imageR = self._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(*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, intensity, self._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.minimum[0] # n.b. makes copy
180 imageG = imageG - self.minimum[1]
181 imageB = imageB - self.minimum[2]
182
183 fac = self.mapIntensityToUint8(self.intensity(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
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.maximum = maximum
239
240 if maximum is None:
241 self._range = None
242 else:
243 assert maximum - minimum != 0, "minimum and maximum values must not be equal"
244 self._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,
256 self._uint8Max/intensity, self._uint8Max/self._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 = self._uint8Max/Q # gradient at origin is self._slope
308 else:
309 frac = 0.1 # gradient estimated using frac*range is _slope
310 self._slope = frac*self._uint8Max/np.arcsinh(frac*Q)
311
312 self._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)*self._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 # support self.makeRgbImage()
383
384
385def 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
432def 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
449def 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)
__init__(self, minimum=None, maximum=None, image=None)
intensity(self, imageR, imageG, imageB)
__init__(self, minimum=None, image=None)
_convertImagesToUint8(self, imageR, imageG, imageB)
makeRgbImage(self, imageR=None, imageG=None, imageB=None, xSize=None, ySize=None, rescaleFactor=None)
__init__(self, image, nSamples=1000, contrast=0.25)
makeRGB(imageR, imageG=None, imageB=None, minimum=0, dataRange=5, Q=8, fileName=None, saturatedBorderWidth=0, saturatedPixelValue=None, xSize=None, ySize=None, rescaleFactor=None)
computeIntensity(imageR, imageG=None, imageB=None)
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.
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:361