Coverage for python / lsst / pipe / tasks / _fallback_asinhmapping.py: 0%
79 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 09:04 +0000
1__all__ = ["AsinhMapping",]
3import numpy as np
5from astropy.visualization.lupton_rgb import compute_intensity
8class Mapping:
9 """
10 Baseclass to map red, blue, green intensities into uint8 values.
11 Parameters
12 ----------
13 minimum : float or sequence(3)
14 Intensity that should be mapped to black (a scalar or array for R, G, B).
15 image : ndarray, optional
16 An image used to calculate some parameters of some mappings.
17 """
19 def __init__(self, minimum=None, image=None):
20 self._uint8Max = float(np.iinfo(np.uint8).max)
22 try:
23 len(minimum)
24 except TypeError:
25 minimum = 3 * [minimum]
26 if len(minimum) != 3:
27 raise ValueError("please provide 1 or 3 values for minimum.")
29 self.minimum = minimum
30 self._image = np.asarray(image)
32 def make_rgb_image(self, image_r, image_g, image_b):
33 """
34 Convert 3 arrays, image_r, image_g, and image_b into an 8-bit RGB image.
35 Parameters
36 ----------
37 image_r : ndarray
38 Image to map to red.
39 image_g : ndarray
40 Image to map to green.
41 image_b : ndarray
42 Image to map to blue.
43 Returns
44 -------
45 RGBimage : ndarray
46 RGB (integer, 8-bits per channel) color image as an NxNx3 numpy array.
47 """
48 image_r = np.asarray(image_r)
49 image_g = np.asarray(image_g)
50 image_b = np.asarray(image_b)
52 if (image_r.shape != image_g.shape) or (image_g.shape != image_b.shape):
53 msg = "The image shapes must match. r: {}, g: {} b: {}"
54 raise ValueError(msg.format(image_r.shape, image_g.shape, image_b.shape))
56 return np.dstack(
57 self._convert_images_to_uint8(image_r, image_g, image_b)
58 ).astype(np.uint8)
60 def intensity(self, image_r, image_g, image_b):
61 """
62 Return the total intensity from the red, blue, and green intensities.
63 This is a naive computation, and may be overridden by subclasses.
64 Parameters
65 ----------
66 image_r : ndarray
67 Intensity of image to be mapped to red; or total intensity if
68 ``image_g`` and ``image_b`` are None.
69 image_g : ndarray, optional
70 Intensity of image to be mapped to green.
71 image_b : ndarray, optional
72 Intensity of image to be mapped to blue.
73 Returns
74 -------
75 intensity : ndarray
76 Total intensity from the red, blue and green intensities, or
77 ``image_r`` if green and blue images are not provided.
78 """
79 return compute_intensity(image_r, image_g, image_b)
81 def map_intensity_to_uint8(self, I): # noqa: E741
82 """
83 Return an array which, when multiplied by an image, returns that image
84 mapped to the range of a uint8, [0, 255] (but not converted to uint8).
85 The intensity is assumed to have had minimum subtracted (as that can be
86 done per-band).
87 Parameters
88 ----------
89 I : ndarray
90 Intensity to be mapped.
91 Returns
92 -------
93 mapped_I : ndarray
94 ``I`` mapped to uint8
95 """
96 with np.errstate(invalid="ignore", divide="ignore"):
97 return np.clip(I, 0, self._uint8Max)
99 def _convert_images_to_uint8(self, image_r, image_g, image_b):
100 """
101 Use the mapping to convert images image_r, image_g, and image_b to a triplet of uint8 images.
102 """
103 image_r = image_r - self.minimum[0] # n.b. makes copy
104 image_g = image_g - self.minimum[1]
105 image_b = image_b - self.minimum[2]
107 fac = self.map_intensity_to_uint8(self.intensity(image_r, image_g, image_b))
109 image_rgb = [image_r, image_g, image_b]
110 for c in image_rgb:
111 c *= fac
112 with np.errstate(invalid="ignore"):
113 c[c < 0] = 0 # individual bands can still be < 0, even if fac isn't
115 pixmax = self._uint8Max
116 # copies -- could work row by row to minimise memory usage
117 r0, g0, b0 = image_rgb
119 # n.b. np.where can't and doesn't short-circuit
120 with np.errstate(invalid="ignore", divide="ignore"):
121 for i, c in enumerate(image_rgb):
122 c = np.where(
123 r0 > g0,
124 np.where(
125 r0 > b0,
126 np.where(r0 >= pixmax, c * pixmax / r0, c),
127 np.where(b0 >= pixmax, c * pixmax / b0, c),
128 ),
129 np.where(
130 g0 > b0,
131 np.where(g0 >= pixmax, c * pixmax / g0, c),
132 np.where(b0 >= pixmax, c * pixmax / b0, c),
133 ),
134 ).astype(np.uint8)
135 c[c > pixmax] = pixmax
137 image_rgb[i] = c
139 return image_rgb
142class LinearMapping(Mapping):
143 """
144 A linear map map of red, blue, green intensities into uint8 values.
145 A linear stretch from [minimum, maximum].
146 If one or both are omitted use image min and/or max to set them.
147 Parameters
148 ----------
149 minimum : float
150 Intensity that should be mapped to black (a scalar or array for R, G, B).
151 maximum : float
152 Intensity that should be mapped to white (a scalar).
153 """
155 def __init__(self, minimum=None, maximum=None, image=None):
156 if minimum is None or maximum is None:
157 if image is None:
158 raise ValueError(
159 "you must provide an image if you don't "
160 "set both minimum and maximum"
161 )
162 if minimum is None:
163 minimum = image.min()
164 if maximum is None:
165 maximum = image.max()
167 Mapping.__init__(self, minimum=minimum, image=image)
168 self.maximum = maximum
170 if maximum is None:
171 self._range = None
172 else:
173 if maximum == minimum:
174 raise ValueError("minimum and maximum values must not be equal")
175 self._range = float(maximum - minimum)
177 def map_intensity_to_uint8(self, I): # noqa: E741
178 # n.b. np.where can't and doesn't short-circuit
179 with np.errstate(invalid="ignore", divide="ignore"):
180 return np.where(
181 I <= 0,
182 0,
183 np.where(
184 I >= self._range, self._uint8Max / I, self._uint8Max / self._range
185 ),
186 )
189class AsinhMapping(Mapping):
190 """
191 A mapping for an asinh stretch (preserving colours independent of brightness).
192 x = asinh(Q (I - minimum)/stretch)/Q
193 This reduces to a linear stretch if Q == 0
194 See https://ui.adsabs.harvard.edu/abs/2004PASP..116..133L
195 Parameters
196 ----------
197 minimum : float
198 Intensity that should be mapped to black (a scalar or array for R, G, B).
199 stretch : float
200 The linear stretch of the image.
201 Q : float
202 The asinh softening parameter.
203 """
205 def __init__(self, minimum, stretch, Q=8):
206 Mapping.__init__(self, minimum)
208 # 32bit floating point machine epsilon; sys.float_info.epsilon is 64bit
209 epsilon = 1.0 / 2**23
210 if abs(Q) < epsilon:
211 Q = 0.1
212 else:
213 Qmax = 1e10
214 if Q > Qmax:
215 Q = Qmax
217 frac = 0.1 # gradient estimated using frac*stretch is _slope
218 self._slope = frac * self._uint8Max / np.arcsinh(frac * Q)
220 self._soften = Q / float(stretch)
222 def map_intensity_to_uint8(self, I): # noqa: E741
223 # n.b. np.where can't and doesn't short-circuit
224 with np.errstate(invalid="ignore", divide="ignore"):
225 return np.where(I <= 0, 0, np.arcsinh(I * self._soften) * self._slope / I)