feat: add optional dynamic pixel range operations to tile factory (#9)

The GDAL tile factory can now automatically rescale pixel values to match a desired output type. This is commonly used when converting panchromatic imagery which often has 11-16 bits per pixel into an 8-bits per pixel representation for visualization. These operations look at the histogram of the input image pixel values and then map them to the output range.

src/aws/osml/gdal/__init__.py

@@ -10,7 +10,7 @@
 from .gdal_utils import get_image_extension, get_type_and_scales, load_gdal_dataset
 from .nitf_des_accessor import NITFDESAccessor
 from .sensor_model_factory import ChippedImageInfoFacade, SensorModelFactory, SensorModelTypes
-from .typing import GDALCompressionOptions, GDALImageFormats
+from .typing import GDALCompressionOptions, GDALImageFormats, RangeAdjustmentType
 
 __all__ = [
     "set_gdal_default_configuration",
@@ -21,6 +21,7 @@
     "GDALConfigEnv",
     "GDALDigitalElevationModelTileFactory",
     "GDALImageFormats",
+    "RangeAdjustmentType",
     "NITFDESAccessor",
     "ChippedImageInfoFacade",
     "SensorModelFactory",

src/aws/osml/gdal/dynamic_range_adjustment.py

@@ -0,0 +1,86 @@
+from typing import List
+
+
+class DRAParameters:
+    """
+    This class manages a set of parameters used to perform a Dynamic Range Adjustment that is applied when
+    converting imagery pixel values (e.g. 11-bit per pixel panchromatic imagery to an 8-bit per pixel grayscale).
+    """
+
+    def __init__(
+        self, suggested_min_value: float, suggested_max_value: float, actual_min_value: float, actual_max_value: float
+    ):
+        """
+        Constructor for this class.
+
+        :param suggested_min_value: suggested minimum value of the relevant pixel range
+        :param suggested_max_value: suggested maximum value of the relevant pixel range
+        :param actual_min_value: actual minimum value of pixels in the image
+        :param actual_max_value: actual maximum value of pixels in the image
+        """
+        self.suggested_min_value = suggested_min_value
+        self.suggested_max_value = suggested_max_value
+        self.actual_min_value = actual_min_value
+        self.actual_max_value = actual_max_value
+
+    @staticmethod
+    def from_counts(
+        counts: List[float], min_percentage: float = 0.02, max_percentage: float = 0.98, a: float = 0.2, b: float = 0.4
+    ) -> "DRAParameters":
+        """
+        This static factory method computes a new set of DRA parameters given a histogram of pixel values.
+
+        :param counts: histogram of the pixel values
+        :param min_percentage: set point for low intensity pixels that may be outliers
+        :param max_percentage: set point for high intensity pixels that may be outliers
+        :param a: weighting factor for the low intensity range
+        :param b: weighting factor for the high intensity range
+        :return: a set of DRA parameters containing recommended and actual ranges of values
+        """
+        num_histogram_bins = len(counts)
+
+        # Find the first and last non-zero counts
+        actual_min_value = 0
+        while actual_min_value < num_histogram_bins and counts[actual_min_value] == 0:
+            actual_min_value += 1
+
+        actual_max_value = num_histogram_bins - 1
+        while actual_max_value > 0 and counts[actual_max_value] == 0:
+            actual_max_value -= 1
+
+        # Compute the cumulative distribution
+        cumulative_counts = counts.copy()
+        for i in range(1, len(cumulative_counts)):
+            cumulative_counts[i] = cumulative_counts[i] + cumulative_counts[i - 1]
+
+        # Find the values that exclude the lowest and highest percentages of the counts.
+        # This identifies the range that contains most of the pixels while excluding outliers.
+        max_counts = cumulative_counts[-1]
+        low_threshold = min_percentage * max_counts
+        e_min = 0
+        while cumulative_counts[e_min] < low_threshold:
+            e_min += 1
+
+        high_threshold = max_percentage * max_counts
+        e_max = num_histogram_bins - 1
+        while cumulative_counts[e_max] > high_threshold:
+            e_max -= 1
+
+        min_value = max([actual_min_value, e_min - a * (e_max - e_min)])
+        max_value = min([actual_max_value, e_max + b * (e_max - e_min)])
+
+        return DRAParameters(
+            suggested_min_value=min_value,
+            suggested_max_value=max_value,
+            actual_min_value=actual_min_value,
+            actual_max_value=actual_max_value,
+        )
+
+    def __repr__(self):
+        return (
+            f"DRAParameters(min_value={self.suggested_min_value}, "
+            f"max_value={self.suggested_max_value}, "
+            f"e_first={self.actual_min_value}, "
+            f"e_last={self.actual_max_value}, "
+            f")"
+        )

src/aws/osml/gdal/gdal_utils.py

@@ -7,7 +7,9 @@
 
 from aws.osml.photogrammetry import SensorModel
 
+from .dynamic_range_adjustment import DRAParameters
 from .sensor_model_factory import SensorModelFactory, SensorModelTypes
+from .typing import RangeAdjustmentType
 
 logger = logging.getLogger(__name__)
 
@@ -65,54 +67,100 @@ def load_gdal_dataset(image_path: str) -> Tuple[gdal.Dataset, Optional[SensorMod
     return ds, sensor_model
 
 
-def get_type_and_scales(raster_dataset: gdal.Dataset) -> Tuple[int, List[List[int]]]:
+def get_minmax_for_type(gdal_type: int) -> Tuple[float, float]:
+    """
+    This function computes the minimum and maximum values that can be stored in a given GDAL pixel type.
+
+    :param gdal_type: the pixel type
+    :return: tuple of min, max values
+    """
+    min_value = 0
+    max_value = 255
+    if gdal_type == gdalconst.GDT_Byte:
+        min_value = 0
+        max_value = 2**8 - 1
+    elif gdal_type == gdalconst.GDT_UInt16:
+        min_value = 0
+        max_value = 2**16 - 1
+    elif gdal_type == gdalconst.GDT_Int16:
+        min_value = -(2**15)
+        max_value = 2**15 - 1
+    elif gdal_type == gdalconst.GDT_UInt32:
+        min_value = 0
+        max_value = 2**32 - 1
+    elif gdal_type == gdalconst.GDT_Int32:
+        min_value = -(2**31)
+        max_value = 2**31 - 1
+    elif gdal_type == gdalconst.GDT_UInt64:
+        min_value = 0
+        max_value = 2**64 - 1
+    elif gdal_type == gdalconst.GDT_Int64:
+        min_value = -(2**63)
+        max_value = 2**63 - 1
+    elif gdal_type == gdalconst.GDT_Float32:
+        min_value = -3.4028235e38
+        max_value = 3.4028235e38
+    elif gdal_type == gdalconst.GDT_Float64:
+        min_value = -1.7976931348623157e308
+        max_value = 1.7976931348623157e308
+    else:
+        logger.warning("Image uses unsupported GDAL datatype {}. Defaulting to [0,255] range".format(gdal_type))
+
+    return min_value, max_value
+
+
+def get_type_and_scales(
+    raster_dataset: gdal.Dataset,
+    desired_output_type: Optional[int] = None,
+    range_adjustment: RangeAdjustmentType = RangeAdjustmentType.NONE,
+) -> Tuple[int, List[List[int]]]:
     """
     Get type and scales of a provided raster dataset
 
     :param raster_dataset: the raster dataset containing the region
+    :param desired_output_type: type to be output after dynamic range adjustments
+    :param range_adjustment: the type of pixel scaling effort to
 
     :return: a tuple containing type and scales
     """
     scale_params = []
-    num_bands = raster_dataset.RasterCount
     output_type = gdalconst.GDT_Byte
-    min = 0
-    max = 255
+    num_bands = raster_dataset.RasterCount
     for band_num in range(1, num_bands + 1):
         band = raster_dataset.GetRasterBand(band_num)
-        output_type = band.DataType
-        if output_type == gdalconst.GDT_Byte:
-            min = 0
-            max = 2**8 - 1
-        elif output_type == gdalconst.GDT_UInt16:
-            min = 0
-            max = 2**16 - 1
-        elif output_type == gdalconst.GDT_Int16:
-            min = -(2**15)
-            max = 2**15 - 1
-        elif output_type == gdalconst.GDT_UInt32:
-            min = 0
-            max = 2**32 - 1
-        elif output_type == gdalconst.GDT_Int32:
-            min = -(2**31)
-            max = 2**31 - 1
-        # TODO: Add these 64-bit cases in once GDAL is upgraded to a version that supports them
-        # elif output_type == gdalconst.GDT_UInt64:
-        #    min = 0
-        #    max = 2**64-1
-        # elif output_type == gdalconst.GDT_Int64:
-        #    min = -2**63
-        #    max = 2**63-1
-        elif output_type == gdalconst.GDT_Float32:
-            min = -3.4028235e38
-            max = 3.4028235e38
-        elif output_type == gdalconst.GDT_Float64:
-            min = -1.7976931348623157e308
-            max = 1.7976931348623157e308
-        else:
-            logger.warning("Image uses unsupported GDAL datatype {}. Defaulting to [0,255] range".format(output_type))
+        band_type = band.DataType
+        min_value, max_value = get_minmax_for_type(band_type)
 
-        scale_params.append([min, max, min, max])
+        if desired_output_type is None:
+            output_type = band_type
+            output_min = min_value
+            output_max = max_value
+        else:
+            output_type = desired_output_type
+            output_min, output_max = get_minmax_for_type(desired_output_type)
+
+        # If a range adjustment is requested compute the range of source pixel values that will be mapped to the full
+        # output range.
+        selected_min = min_value
+        selected_max = max_value
+        if range_adjustment is not RangeAdjustmentType.NONE:
+            num_buckets = int(max_value - min_value)
+            if band_type == gdalconst.GDT_Float32 or band_type == gdalconst.GDT_Float64:
+                num_buckets = 255
+            dra = DRAParameters.from_counts(
+                band.GetHistogram(buckets=num_buckets, max=max_value, min=min_value, include_out_of_range=1, approx_ok=1)
+            )
+            if range_adjustment == RangeAdjustmentType.DRA:
+                selected_min = dra.suggested_min_value
+                selected_max = dra.suggested_max_value
+            elif range_adjustment == RangeAdjustmentType.MINMAX:
+                selected_min = dra.actual_min_value
+                selected_max = dra.actual_max_value
+            else:
+                logger.warning(f"Unknown range adjustment selected {range_adjustment}. Skipping.")
+
+        band_scale_parameters = [selected_min, selected_max, output_min, output_max]
+        scale_params.append(band_scale_parameters)
 
     return output_type, scale_params
 

src/aws/osml/gdal/typing.py

@@ -22,3 +22,19 @@ class GDALImageFormats(str, Enum):
     JPEG = "JPEG"
     PNG = "PNG"
     GTIFF = "GTiff"
+
+
+class RangeAdjustmentType(str, Enum):
+    """
+    Enumeration defining ways to scale raw image pixels to an output value range.
+
+    - NONE indicates that the full range available to the input type will be used.
+    - MINMAX chooses the portion of the input range that actually contains values.
+    - DRA is a dynamic range adjustment that attempts to select the most important portion of the input range.
+      It differs from MINMAX in that it can exclude outliers to reduce the impact of unusually bright/dark
+      spots in an image.
+    """
+
+    NONE = "NONE"
+    MINMAX = "MINMAX"
+    DRA = "DRA"

src/aws/osml/image_processing/gdal_tile_factory.py

@@ -6,7 +6,7 @@
 from defusedxml import ElementTree
 from osgeo import gdal, gdalconst
 
-from aws.osml.gdal import GDALCompressionOptions, GDALImageFormats, NITFDESAccessor, get_type_and_scales
+from aws.osml.gdal import GDALCompressionOptions, GDALImageFormats, NITFDESAccessor, RangeAdjustmentType, get_type_and_scales
 from aws.osml.photogrammetry import ImageCoordinate, SensorModel
 
 from .sicd_updater import SICDUpdater
@@ -26,6 +26,8 @@ def __init__(
         sensor_model: Optional[SensorModel] = None,
         tile_format: GDALImageFormats = GDALImageFormats.NITF,
         tile_compression: GDALCompressionOptions = GDALCompressionOptions.NONE,
+        output_type: Optional[int] = None,
+        range_adjustment: RangeAdjustmentType = RangeAdjustmentType.NONE,
     ):
         """
         Constructs a new factory capable of producing tiles from a given GDAL raster dataset.
@@ -34,6 +36,8 @@ def __init__(
         :param sensor_model: the sensor model providing mensuration support for this image
         :param tile_format: the output tile format
         :param tile_compression: the output tile compression
+        :param output_type: the GDAL pixel type in the output tile
+        :param range_adjustment: the type of scaling used to convert raw pixel values to the output range
         """
         self.tile_format = tile_format
         self.tile_compression = tile_compression
@@ -42,6 +46,8 @@ def __init__(
         self.des_accessor = None
         self.sicd_updater = None
         self.sicd_des_header = None
+        self.range_adjustment = range_adjustment
+        self.output_type = output_type
 
         if self.raster_dataset.GetDriver().ShortName == "NITF":
             xml_des = self.raster_dataset.GetMetadata("xml:DES")
@@ -58,6 +64,8 @@ def __init__(
                 self.sicd_des_header = self.des_accessor.extract_des_header(sicd_des)
                 self.sicd_updater = SICDUpdater(sicd_metadata)
 
+        self.default_gdal_translate_kwargs = self._create_gdal_translate_kwargs()
+
     def create_encoded_tile(self, src_window: List[int]) -> Optional[bytearray]:
         """
         This method cuts a tile from the full image, updates the metadata as needed, and finally compresses/encodes
@@ -71,10 +79,10 @@ def create_encoded_tile(self, src_window: List[int]) -> Optional[bytearray]:
         # Use the request and metadata from the raster dataset to create a set of keyword
         # arguments for the gdal.Translate() function. This will configure that function to
         # create image tiles using the format, compression, etc. requested by the client.
-        gdal_translate_kwargs = self._create_gdal_translate_kwargs()
+        gdal_translate_kwargs = self.default_gdal_translate_kwargs.copy()
 
         # Create a new IGEOLO value based on the corner points of this tile
-        if self.sensor_model is not None:
+        if self.sensor_model is not None and self.tile_format == GDALImageFormats.NITF:
             gdal_translate_kwargs["creationOptions"].append("ICORDS=G")
             gdal_translate_kwargs["creationOptions"].append("IGEOLO=" + self.create_new_igeolo(src_window))
 
@@ -148,9 +156,10 @@ def _create_gdal_translate_kwargs(self) -> Dict[str, Any]:
 
         :return: Dict[str, any] = the dictionary of translate keyword arguments
         """
-        # Figure out what type of image this is and calculate a scale that does not force any range
-        # remapping
-        output_type, scale_params = get_type_and_scales(self.raster_dataset)
+        # Figure out what type of image this is and calculate a scale to map input pixels to the output type
+        output_type, scale_params = get_type_and_scales(
+            self.raster_dataset, desired_output_type=self.output_type, range_adjustment=self.range_adjustment
+        )
 
         gdal_translate_kwargs = {
             "scaleParams": scale_params,

test/aws/osml/gdal/test_dynamic_range_adjustment.py

@@ -0,0 +1,18 @@
+import unittest
+
+
+class TestDRAParameters(unittest.TestCase):
+    def test_from_counts(self):
+        from aws.osml.gdal.dynamic_range_adjustment import DRAParameters
+
+        counts = [0] * 1024
+        counts[1:99] = [1] * (99 - 1)
+        counts[100:400] = [200] * (400 - 100)
+        counts[1022] = 1
+
+        dra_parameters = DRAParameters.from_counts(counts=counts)
+
+        self.assertEquals(dra_parameters.actual_min_value, 1)
+        self.assertEquals(dra_parameters.actual_max_value, 1022)
+        self.assertAlmostEqual(dra_parameters.suggested_min_value, 47, delta=1)
+        self.assertAlmostEquals(dra_parameters.suggested_max_value, 506, delta=1)