Add the reinhard enhancements

doc/source/enhancements.rst

@@ -66,6 +66,19 @@ crefl_scaling
 cira_stretch
 ------------
 
+Logarithmic stretch based on a cira recipe.
+
+reinhard_to_srgb
+----------------
+
+Stretch method based on the Reinhard algorithm, using luminance.
+
+The function includes conversion to sRGB colorspace.
+
+    Reinhard, Erik & Stark, Michael & Shirley, Peter & Ferwerda, James. (2002).
+    Photographic Tone Reproduction For Digital Images. ACM Transactions on Graphics.
+    :doi: `21. 10.1145/566654.566575`
+
 lookup
 ------
 

satpy/enhancements/__init__.py

@@ -140,6 +140,54 @@ def func(band_data):
     return apply_enhancement(img.data, func)
 
 
+def reinhard_to_srgb(img, saturation=1.25, white=100, **kwargs):
+    """Stretch method based on the Reinhard algorithm, using luminance.
+
+    Args:
+        saturation: Saturation enhancement factor. Less is grayer. Neutral is 1.
+        white: the reflectance luminance to set to white (in %).
+
+
+    Reinhard, Erik & Stark, Michael & Shirley, Peter & Ferwerda, James. (2002).
+    Photographic Tone Reproduction For Digital Images. ACM Transactions on Graphics.
+    :doi: `21. 10.1145/566654.566575`
+    """
+    with xr.set_options(keep_attrs=True):
+        # scale the data to [0, 1] interval
+        rgb = img.data / 100
+        white /= 100
+
+        # extract color components
+        r = rgb.sel(bands='R').data
+        g = rgb.sel(bands='G').data
+        b = rgb.sel(bands='B').data
+
+        # saturate
+        luma = _compute_luminance_from_rgb(r, g, b)
+        rgb = (luma + (rgb - luma) * saturation).clip(0)
+
+        # reinhard
+        reinhard_luma = (luma / (1 + luma)) * (1 + luma/(white**2))
+        coef = reinhard_luma / luma
+        rgb = rgb * coef
+
+        # srgb gamma
+        rgb.data = _srgb_gamma(rgb.data)
+        img.data = rgb
+
+    return img.data
+
+
+def _compute_luminance_from_rgb(r, g, b):
+    """Compute the luminance of the image."""
+    return r * 0.2126 + g * 0.7152 + b * 0.0722
+
+
+def _srgb_gamma(arr):
+    """Apply the srgb gamma."""
+    return da.where(arr < 0.0031308, arr * 12.92, 1.055 * arr ** 0.41666 - 0.055)
+
+
 def _lookup_delayed(luts, band_data):
     # can't use luts.__getitem__ for some reason
     return luts[band_data]

satpy/tests/enhancement_tests/test_enhancements.py

@@ -66,6 +66,19 @@ def test_cira_stretch(self):
             [1.05181359, 1.11651012, 1.16635571, 1.20691137, 1.24110186]]])
         self._test_enhancement(cira_stretch, self.ch1, expected)
 
+    def test_reinhard(self):
+        """Test the reinhard algorithm."""
+        from satpy.enhancements import reinhard_to_srgb
+        expected = np.array([[[np.nan, 0., 0., 0.93333793, 1.29432402],
+                              [1.55428709, 1.76572249, 1.94738635, 2.10848544, 2.25432809]],
+
+                             [[np.nan, 0., 0., 0.93333793, 1.29432402],
+                              [1.55428709, 1.76572249, 1.94738635, 2.10848544, 2.25432809]],
+
+                             [[np.nan, 0., 0., 0.93333793, 1.29432402],
+                              [1.55428709, 1.76572249, 1.94738635, 2.10848544, 2.25432809]]])
+        self._test_enhancement(reinhard_to_srgb, self.rgb, expected)
+
     def test_lookup(self):
         """Test the lookup enhancement function."""
         from satpy.enhancements import lookup