(Multiband|)Tile Histogram Matching

raster/src/main/scala/geotrellis/raster/matching/HistogramMatching.scala

@@ -0,0 +1,94 @@
+package geotrellis.raster.matching
+
+import geotrellis.raster._
+import geotrellis.raster.histogram.Histogram
+import geotrellis.raster.histogram.StreamingHistogram
+import geotrellis.raster.equalization.HistogramEqualization._
+
+
+/**
+  * Uses the approach given here:
+  * http://fourier.eng.hmc.edu/e161/lectures/contrast_transform/node3.html
+  */
+object HistogramMatching {
+
+  /**
+    * Comparison class for sorting an array of (x, cdf(x))
+    * pairs by cdf(x).
+    */
+  private class BucketComparator extends java.util.Comparator[(Double, Double)] {
+    // Compare the (label, cdf(label)) pairs by their labels
+    def compare(left: (Double, Double), right: (Double, Double)): Int = {
+      if (left._2 < right._2) -1
+      else if (left._2 > right._2) 1
+      else 0
+    }
+  }
+
+  private val cmp = new BucketComparator
+
+  /**
+    * An implementation of the compound transformation referred to in
+    * the citation given above.  The idea is to transform from the
+    * source histogram to an equalized one, then transform from the
+    * equalized one to the target one via the inverse of the
+    * transformation that would equalize the target one.
+    */
+  @inline private def transform(
+    targetCdf: Array[(Double, Double)], fn: (Double => Double)
+  )(x: Double): Double = {
+    val cdfx = fn(x)
+    val i = java.util.Arrays.binarySearch(targetCdf, (0.0, cdfx), cmp)
+
+    if (i < 0) targetCdf(math.min(-(i+1), targetCdf.length-1))._1
+    else targetCdf(i)._1
+  }
+
+  def apply[T1 <: AnyVal, T2 <: AnyVal](
+    tile: Tile,
+    sourceHistogram: Histogram[T1],
+    targetHistogram: Histogram[T2]
+  ): Tile = {
+    val cellType = tile.cellType
+    val localIntensityToCdf = intensityToCdf(cellType, sourceHistogram.cdf)_
+    val localTransform = transform(targetHistogram.cdf, localIntensityToCdf)_
+
+    tile.mapDouble(localTransform)
+  }
+
+  def apply[T <: AnyVal](tile: Tile, targetHistogram: Histogram[T]): Tile =
+    HistogramMatching(
+      tile,
+      StreamingHistogram.fromTile(tile, 1<<17),
+      targetHistogram
+    )
+
+  def apply[T1 <: AnyVal, T2 <: AnyVal](
+    tile: MultibandTile,
+    sourceHistograms: Seq[Histogram[T1]],
+    targetHistograms: Seq[Histogram[T2]]
+  ): MultibandTile =
+    MultibandTile(
+      tile.bands
+        .zip(sourceHistograms.zip(targetHistograms))
+        .map({ case (tile, (source, target)) =>
+          HistogramMatching(tile, source, target)
+        })
+    )
+
+  def apply[T <: AnyVal](
+    tile: MultibandTile,
+    targetHistograms: Seq[Histogram[T]]
+  ): MultibandTile =
+    MultibandTile(
+      tile.bands
+        .zip(targetHistograms)
+        .map({ case (tile, target) =>
+          HistogramMatching(
+            tile,
+            StreamingHistogram.fromTile(tile, 1<<17),
+            target
+          )
+        })
+    )
+}

raster/src/main/scala/geotrellis/raster/matching/MultibandMatchingMethods.scala

@@ -0,0 +1,17 @@
+package geotrellis.raster.matching
+
+import geotrellis.raster.histogram.StreamingHistogram
+import geotrellis.raster.MultibandTile
+import geotrellis.util.MethodExtensions
+
+
+trait MultibandMatchingMethods extends MethodExtensions[MultibandTile] {
+
+  def matchHistogram(targetHistograms: Seq[StreamingHistogram]): MultibandTile =
+    HistogramMatching(self, targetHistograms)
+
+  def matchHistogram(
+    sourceHistograms: Seq[StreamingHistogram],
+    targetHistograms: Seq[StreamingHistogram]
+  ): MultibandTile = HistogramMatching(self, sourceHistograms, targetHistograms)
+}

raster/src/main/scala/geotrellis/raster/matching/SinglebandMatchingMethods.scala

@@ -0,0 +1,10 @@
+package geotrellis.raster.matching
+
+import geotrellis.raster.histogram.StreamingHistogram
+import geotrellis.raster.Tile
+import geotrellis.util.MethodExtensions
+
+
+trait SinglebandMatchingMethods extends MethodExtensions[Tile] {
+  def matchHistogram(histogram: StreamingHistogram): Tile = HistogramMatching(self, histogram)
+}

raster/src/main/scala/geotrellis/raster/package.scala

@@ -60,6 +60,7 @@ package object raster
       with mapalgebra.local.LocalMethods
       with mapalgebra.zonal.ZonalMethods
       with mask.SinglebandTileMaskMethods
+      with matching.SinglebandMatchingMethods
       with merge.SinglebandTileMergeMethods
       with prototype.SinglebandTilePrototypeMethods
       with regiongroup.RegionGroupMethods
@@ -79,6 +80,7 @@ package object raster
       with crop.MultibandTileCropMethods
       with equalization.MultibandEqualizationMethods
       with mask.MultibandTileMaskMethods
+      with matching.MultibandMatchingMethods
       with merge.MultibandTileMergeMethods
       with prototype.MultibandTilePrototypeMethods
       with render.MultibandColorMethods