# Source code for astropy.visualization.lupton_rgb

"""
Combine 3 images to produce a properly-scaled RGB image following Lupton et al. (2004).

The three images must be aligned and have the same pixel scale and size.

"""

import numpy as np
from . import ZScaleInterval

__all__ = ['make_lupton_rgb']

def compute_intensity(image_r, image_g=None, image_b=None):
"""
Return a naive total intensity from the red, blue, and green intensities.

Parameters
----------
image_r : ndarray
Intensity of image to be mapped to red; or total intensity if ``image_g``
and ``image_b`` are None.
image_g : ndarray, optional
Intensity of image to be mapped to green.
image_b : ndarray, optional
Intensity of image to be mapped to blue.

Returns
-------
intensity : ndarray
Total intensity from the red, blue and green intensities, or ``image_r``
if green and blue images are not provided.
"""
if image_g is None or image_b is None:
if not (image_g is None and image_b is None):
raise ValueError("please specify either a single image "
"or red, green, and blue images.")
return image_r

intensity = (image_r + image_g + image_b)/3.0

# Repack into whatever type was passed to us
return np.asarray(intensity, dtype=image_r.dtype)

class Mapping:
"""
Baseclass to map red, blue, green intensities into uint8 values.

Parameters
----------
minimum : float or sequence(3)
Intensity that should be mapped to black (a scalar or array for R, G, B).
image : ndarray, optional
An image used to calculate some parameters of some mappings.
"""

def __init__(self, minimum=None, image=None):
self._uint8Max = float(np.iinfo(np.uint8).max)

try:
len(minimum)
except TypeError:
minimum = 3*[minimum]
if len(minimum) != 3:
raise ValueError("please provide 1 or 3 values for minimum.")

self.minimum = minimum
self._image = np.asarray(image)

def make_rgb_image(self, image_r, image_g, image_b):
"""
Convert 3 arrays, image_r, image_g, and image_b into an 8-bit RGB image.

Parameters
----------
image_r : ndarray
Image to map to red.
image_g : ndarray
Image to map to green.
image_b : ndarray
Image to map to blue.

Returns
-------
RGBimage : ndarray
RGB (integer, 8-bits per channel) color image as an NxNx3 numpy array.
"""
image_r = np.asarray(image_r)
image_g = np.asarray(image_g)
image_b = np.asarray(image_b)

if (image_r.shape != image_g.shape) or (image_g.shape != image_b.shape):
msg = "The image shapes must match. r: {}, g: {} b: {}"
raise ValueError(msg.format(image_r.shape, image_g.shape, image_b.shape))

return np.dstack(self._convert_images_to_uint8(image_r, image_g, image_b)).astype(np.uint8)

def intensity(self, image_r, image_g, image_b):
"""
Return the total intensity from the red, blue, and green intensities.
This is a naive computation, and may be overridden by subclasses.

Parameters
----------
image_r : ndarray
Intensity of image to be mapped to red; or total intensity if
``image_g`` and ``image_b`` are None.
image_g : ndarray, optional
Intensity of image to be mapped to green.
image_b : ndarray, optional
Intensity of image to be mapped to blue.

Returns
-------
intensity : ndarray
Total intensity from the red, blue and green intensities, or
``image_r`` if green and blue images are not provided.
"""
return compute_intensity(image_r, image_g, image_b)

def map_intensity_to_uint8(self, I):
"""
Return an array which, when multiplied by an image, returns that image
mapped to the range of a uint8, [0, 255] (but not converted to uint8).

The intensity is assumed to have had minimum subtracted (as that can be
done per-band).

Parameters
----------
I : ndarray
Intensity to be mapped.

Returns
-------
mapped_I : ndarray
``I`` mapped to uint8
"""
with np.errstate(invalid='ignore', divide='ignore'):
return np.clip(I, 0, self._uint8Max)

def _convert_images_to_uint8(self, image_r, image_g, image_b):
"""Use the mapping to convert images image_r, image_g, and image_b to a triplet of uint8 images"""
image_r = image_r - self.minimum[0]  # n.b. makes copy
image_g = image_g - self.minimum[1]
image_b = image_b - self.minimum[2]

fac = self.map_intensity_to_uint8(self.intensity(image_r, image_g, image_b))

image_rgb = [image_r, image_g, image_b]
for c in image_rgb:
c *= fac
with np.errstate(invalid='ignore'):
c[c < 0] = 0                # individual bands can still be < 0, even if fac isn't

pixmax = self._uint8Max
r0, g0, b0 = image_rgb           # copies -- could work row by row to minimise memory usage

with np.errstate(invalid='ignore', divide='ignore'):  # n.b. np.where can't and doesn't short-circuit
for i, c in enumerate(image_rgb):
c = np.where(r0 > g0,
np.where(r0 > b0,
np.where(r0 >= pixmax, c*pixmax/r0, c),
np.where(b0 >= pixmax, c*pixmax/b0, c)),
np.where(g0 > b0,
np.where(g0 >= pixmax, c*pixmax/g0, c),
np.where(b0 >= pixmax, c*pixmax/b0, c))).astype(np.uint8)
c[c > pixmax] = pixmax

image_rgb[i] = c

return image_rgb

class LinearMapping(Mapping):
"""
A linear map map of red, blue, green intensities into uint8 values.

A linear stretch from [minimum, maximum].
If one or both are omitted use image min and/or max to set them.

Parameters
----------
minimum : float
Intensity that should be mapped to black (a scalar or array for R, G, B).
maximum : float
Intensity that should be mapped to white (a scalar).
"""

def __init__(self, minimum=None, maximum=None, image=None):
if minimum is None or maximum is None:
if image is None:
raise ValueError("you must provide an image if you don't "
"set both minimum and maximum")
if minimum is None:
minimum = image.min()
if maximum is None:
maximum = image.max()

Mapping.__init__(self, minimum=minimum, image=image)
self.maximum = maximum

if maximum is None:
self._range = None
else:
if maximum == minimum:
raise ValueError("minimum and maximum values must not be equal")
self._range = float(maximum - minimum)

def map_intensity_to_uint8(self, I):
with np.errstate(invalid='ignore', divide='ignore'):  # n.b. np.where can't and doesn't short-circuit
return np.where(I <= 0, 0,
np.where(I >= self._range, self._uint8Max/I, self._uint8Max/self._range))

class AsinhMapping(Mapping):
"""
A mapping for an asinh stretch (preserving colours independent of brightness)

x = asinh(Q (I - minimum)/stretch)/Q

This reduces to a linear stretch if Q == 0

Parameters
----------

minimum : float
Intensity that should be mapped to black (a scalar or array for R, G, B).
stretch : float
The linear stretch of the image.
Q : float
The asinh softening parameter.
"""

def __init__(self, minimum, stretch, Q=8):
Mapping.__init__(self, minimum)

epsilon = 1.0/2**23            # 32bit floating point machine epsilon; sys.float_info.epsilon is 64bit
if abs(Q) < epsilon:
Q = 0.1
else:
Qmax = 1e10
if Q > Qmax:
Q = Qmax

frac = 0.1                  # gradient estimated using frac*stretch is _slope
self._slope = frac*self._uint8Max/np.arcsinh(frac*Q)

self._soften = Q/float(stretch)

def map_intensity_to_uint8(self, I):
with np.errstate(invalid='ignore', divide='ignore'):  # n.b. np.where can't and doesn't short-circuit
return np.where(I <= 0, 0, np.arcsinh(I*self._soften)*self._slope/I)

class AsinhZScaleMapping(AsinhMapping):
"""
A mapping for an asinh stretch, estimating the linear stretch by zscale.

x = asinh(Q (I - z1)/(z2 - z1))/Q

Parameters
----------
image1 : ndarray or a list of arrays
The image to analyse, or a list of 3 images to be converted to
an intensity image.
image2 : ndarray, optional
the second image to analyse (must be specified with image3).
image3 : ndarray, optional
the third image to analyse (must be specified with image2).
Q : float, optional
The asinh softening parameter. Default is 8.
pedestal : float or sequence(3), optional
The value, or array of 3 values, to subtract from the images; or None.

Notes
-----
pedestal, if not None, is removed from the images when calculating the
zscale stretch, and added back into Mapping.minimum[]
"""

def __init__(self, image1, image2=None, image3=None, Q=8, pedestal=None):
"""
"""

if image2 is None or image3 is None:
if not (image2 is None and image3 is None):
raise ValueError("please specify either a single image "
"or three images.")
image = [image1]
else:
image = [image1, image2, image3]

if pedestal is not None:
try:
len(pedestal)
except TypeError:
pedestal = 3*[pedestal]

if len(pedestal) != 3:
raise ValueError("please provide 1 or 3 pedestals.")

image = list(image)        # needs to be mutable
for i, im in enumerate(image):
if pedestal[i] != 0.0:
image[i] = im - pedestal[i]  # n.b. a copy
else:
pedestal = len(image)*[0.0]

image = compute_intensity(*image)

zscale_limits = ZScaleInterval().get_limits(image)
zscale = LinearMapping(*zscale_limits, image=image)
stretch = zscale.maximum - zscale.minimum[0]  # zscale.minimum is always a triple
minimum = zscale.minimum

for i, level in enumerate(pedestal):
minimum[i] += level

AsinhMapping.__init__(self, minimum, stretch, Q)
self._image = image

[docs]def make_lupton_rgb(image_r, image_g, image_b, minimum=0, stretch=5, Q=8, filename=None): """ Return a Red/Green/Blue color image from up to 3 images using an asinh stretch. The input images can be int or float, and in any range or bit-depth. For a more detailed look at the use of this method, see the document :ref:`astropy:astropy-visualization-rgb`. Parameters ---------- image_r : ndarray Image to map to red. image_g : ndarray Image to map to green. image_b : ndarray Image to map to blue. minimum : float Intensity that should be mapped to black (a scalar or array for R, G, B). stretch : float The linear stretch of the image. Q : float The asinh softening parameter. filename : str Write the resulting RGB image to a file (file type determined from extension). Returns ------- rgb : ndarray RGB (integer, 8-bits per channel) color image as an NxNx3 numpy array. """ asinhMap = AsinhMapping(minimum, stretch, Q) rgb = asinhMap.make_rgb_image(image_r, image_g, image_b) if filename: import matplotlib.image matplotlib.image.imsave(filename, rgb, origin='lower') return rgb