Merge pull request #1668 from jamesmcclain/feature/histogram-equaliza…

…tion

Histogram Equalization

geowave/src/main/scala/geotrellis/spark/io/geowave/GeowaveLayerReader.scala

@@ -36,9 +36,9 @@ import scala.reflect._
 
 
 object GeowaveLayerReader {
-  val geometryFactory = new GeometryFactory
-  val tileClassTag = classTag[Tile]
-  val mbtClassTag = classTag[MultibandTile]
+  private val geometryFactory = new GeometryFactory
+  private val tileClassTag = classTag[Tile]
+  private val mbtClassTag = classTag[MultibandTile]
 
   /**
     * Given a map transform and a keybounds, produce a corresponding

raster-test/src/test/scala/geotrellis/raster/equalization/HistogramEqualizationSpec.scala

@@ -0,0 +1,77 @@
+/*
+ * Copyright (c) 2016 Azavea.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package geotrellis.raster.equalization
+
+import geotrellis.raster._
+
+import org.scalatest._
+
+
+class HistogramEqualizationSpec extends FunSpec with Matchers
+{
+  val _a = List(0)
+  val b  = List(1)
+  val c  = List(2)
+  val d  = List(4)
+  val e  = List(8)
+  val f  = List(16)
+  val g  = List(32)
+  val h  = List(64)
+
+  val data = (_a ++ b ++ c ++ d ++ e ++ f ++ g ++ h)
+
+  describe("Histogram Equalization") {
+
+    it("should work on floating-point rasters") {
+      val tile = DoubleArrayTile(data.map(_.toDouble).toArray, 1, 8).equalize
+      val array = tile.toArrayDouble
+
+      array.head should be (Double.MinValue)
+      array.last should be (Double.MaxValue)
+    }
+
+    it("should work on unsigned integral rasters") {
+      val tile = UShortArrayTile(data.map(_.toShort).toArray, 1, 8, UShortCellType).equalize
+      val array = tile.toArray
+
+      array.head should be (0)
+      array.last should be ((1<<16)-1)
+    }
+
+    it("should work on signed integral rasters") {
+      val tile = ShortArrayTile(data.map(_.toShort).toArray, 1, 8, ShortCellType).equalize
+      val array = tile.toArray
+
+      array.head should be (-(1<<15))
+      array.last should be ((1<<15)-1)
+    }
+
+    it("should work on multiband tiles") {
+      val data1 = (_a ++ b ++ c ++ d ++ List(0,0,0,0))
+      val data2 = (List(0,0,0,0) ++ e ++ f ++ g ++ h)
+
+      val tile1 = ShortArrayTile(data1.map(_.toShort).toArray, 1, 8, ShortCellType)
+      val tile2 = ShortArrayTile(data2.map(_.toShort).toArray, 1, 8, ShortCellType)
+      val tile = ArrayMultibandTile(tile1, tile2).equalize
+      val array = tile.bands.flatMap(_.toArray)
+
+      array.head should be (-(1<<15))
+      array.last should be ((1<<15)-1)
+    }
+
+  }
+}

raster/src/main/scala/geotrellis/raster/density/KernelDensity.scala

@@ -22,24 +22,26 @@ import geotrellis.raster.rasterize._
 import geotrellis.raster.mapalgebra.focal.Kernel
 
 /**
- * Object containing functions pertaining to kernel density estimation.
- */
+  * Object containing functions pertaining to kernel density
+  * estimation.
+  */
 object KernelDensity {
   /**
-    * Computes a Density raster based on the Kernel and set of points provided.
-    * Defaults to IntConstantNoDataCellType.
+    * Computes a Density raster based on the Kernel and set of points
+    * provided.  Defaults to IntConstantNoDataCellType.
     */
   def apply(points: Traversable[PointFeature[Int]],
             kernel: Kernel,
             rasterExtent: RasterExtent): Tile =
     apply(points, kernel, rasterExtent, IntConstantNoDataCellType)
 
   /**
-    * Computes a Density raster based on the Kernel and set of int point features provided.
+    * Computes a Density raster based on the Kernel and set of int
+    * point features provided.
     *
     * @param      points           Sequence of point features who's values will be used to
     *                              compute the density.
-    * @param      kernel           [[Kernel]] to be used in the computation.
+    * @param      kernel           [[geotrellis.raster.mapalgebra.focal.Kernel]] to be used in the computation.
     * @param      rasterExtent     Raster extent of the resulting raster.
     * @param      cellType         CellType of the resulting tile.
     *
@@ -72,11 +74,12 @@ object KernelDensity {
     apply(points, kernel, rasterExtent, DoubleConstantNoDataCellType)
 
   /**
-    * Computes a Density raster based on the Kernel and set of double point features provided.
+    * Computes a Density raster based on the Kernel and set of double
+    * point features provided.
     *
     * @param      points           Sequence of point features who's values will be used to
     *                              compute the density.
-    * @param      kernel           [[Kernel]] to be used in the computation.
+    * @param      kernel           [[geotrellis.raster.mapalgebra.focal.Kernel]] to be used in the computation.
     * @param      rasterExtent     Raster extent of the resulting raster.
     * @param      cellType         CellType of the resulting tile.
     *

raster/src/main/scala/geotrellis/raster/equalization/HistogramEqualization.scala

@@ -0,0 +1,119 @@
+package geotrellis.raster.equalization
+
+import geotrellis.raster._
+import geotrellis.raster.histogram.Histogram
+import geotrellis.raster.histogram.StreamingHistogram
+
+
+/**
+  * Uses the approach given here:
+  * http://www.math.uci.edu/icamp/courses/math77c/demos/hist_eq.pdf
+  */
+object HistogramEqualization {
+
+  /**
+    * Comparison class for sorting an array of (label, cdf(label))
+    * pairs by their label.
+    */
+  private class _cmp 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._1 < right._1) -1
+      else if (left._1 > right._1) 1
+      else 0
+    }
+  }
+
+  private val cmp = new _cmp()
+
+  /**
+    * An implementation of the transformation T referred to in the
+    * citation given above.
+    */
+  private def _T(cellType: CellType, cdf: Array[(Double, Double)])(x: Double): Double = {
+    val i = java.util.Arrays.binarySearch(cdf, (x, 0.0), cmp)
+    val bits = cellType.bits
+    val smallestCdf = cdf(0)._2
+    val rawCdf =
+      if (x < cdf(0)._1) { // x is smaller than any label in the array
+        smallestCdf // there is almost no mass here, so report 0.0
+      }
+      else if (-1 * i > cdf.length) { // x is larger than any label in the array
+        1.0 // there is almost no mass here, so report 1.0
+      }
+      else if (i >= 0) { // i is the location of the label x in the array
+        cdf(i)._2
+      }
+      else { // x is between two labels in the array
+        val j = (-1 * i - 2)
+        val cdf0 = cdf(j+0)._2
+        val cdf1 = cdf(j+1)._2
+        val t = (x - cdf0) / (cdf1 - cdf0)
+        (1.0-t)*cdf0 + t*cdf1
+      }
+    val normalizedCdf = math.max(0.0, math.min(1.0, (rawCdf - smallestCdf) / (1.0 - smallestCdf)))
+
+    cellType match {
+      case _: FloatCells => (Float.MaxValue * (2*normalizedCdf - 1.0))
+      case _: DoubleCells => (Double.MaxValue * (2*normalizedCdf - 1.0))
+      case _: BitCells | _: UByteCells | _: UShortCells =>
+        ((1<<bits) - 1) * normalizedCdf
+      case _: ByteCells | _: ShortCells | _: IntCells =>
+        (((1<<bits) - 1) * normalizedCdf) - (1<<(bits-1))
+    }
+  }
+
+  /**
+    * Given a [[Tile]], return a Tile with an equalized histogram.
+    *
+    * @param  tile  A singleband tile
+    * @return       A singleband tile with improved contrast
+    */
+  def apply(tile: Tile): Tile = {
+    val histogram = StreamingHistogram.fromTile(tile, 1<<17)
+    HistogramEqualization(tile, histogram)
+  }
+
+  /**
+    * Given a [[Tile]] and a [[Histogram]], return a Tile with an
+    * equalized histogram.
+    *
+    * @param  tile       A singleband tile
+    * @param  histogram  The histogram of the tile
+    * @return            A singleband tile with improved contrast
+    */
+  def apply[T <: AnyVal](tile: Tile, histogram: Histogram[T]): Tile = {
+    val T = _T(tile.cellType, histogram.cdf)_
+    tile.mapDouble(T)
+  }
+
+  /**
+    * Given a [[MultibandTile]], return a MultibandTile whose bands
+    * all have equalized histograms.
+    *
+    * @param  tile  A multiband tile
+    * @return       A multiband tile with improved contrast
+    */
+  def apply(tile: MultibandTile): MultibandTile = {
+    val histograms = tile.bands.map({ tile => StreamingHistogram.fromTile(tile, 1<<17) })
+    HistogramEqualization(tile, histograms)
+  }
+
+  /**
+    * Given a [[MultibandTile]] and a [[Histogram]] for each of its
+    * bands, return a MultibandTile whose bands all have equalized
+    * histograms.
+    *
+    * @param  tile        A multiband tile
+    * @param  histograms  A sequence of histograms, one for each band
+    * @return             A multiband tile with improved contrast
+    */
+  def apply[T <: AnyVal](tile: MultibandTile, histograms: Seq[Histogram[T]]): MultibandTile = {
+    MultibandTile(
+      tile.bands
+        .zip(histograms)
+        .map({ case (tile, histogram) => HistogramEqualization(tile, histogram) })
+    )
+  }
+
+}

raster/src/main/scala/geotrellis/raster/equalization/MultibandEqualizationMethods.scala

@@ -0,0 +1,24 @@
+package geotrellis.raster.equalization
+
+import geotrellis.raster.histogram.StreamingHistogram
+import geotrellis.raster.MultibandTile
+import geotrellis.util.MethodExtensions
+
+
+trait MultibandEqualizationMethods extends MethodExtensions[MultibandTile] {
+
+  /**
+    * Equalize the histograms of the bands of this [[MultibandTile]].
+    *
+    * @param  histogram  A sequence of [[StreamingHistogram]] objects, one for each band
+    * @return            A multiband tile whose bands have equalized histograms
+    */
+  def equalize(histograms: Array[StreamingHistogram]): MultibandTile = HistogramEqualization(self, histograms)
+
+  /**
+    * Equalize the histograms of the bands of this [[MultibandTile]].
+    *
+    * @return  A multiband tile whose bands have equalized histograms
+    */
+  def equalize(): MultibandTile = HistogramEqualization(self)
+}

raster/src/main/scala/geotrellis/raster/equalization/SinglebandEqualizationMethods.scala

@@ -0,0 +1,22 @@
+package geotrellis.raster.equalization
+
+import geotrellis.raster.histogram.StreamingHistogram
+import geotrellis.raster.Tile
+import geotrellis.util.MethodExtensions
+
+
+trait SinglebandEqualizationMethods extends MethodExtensions[Tile] {
+
+  /**
+    * Given a [[StreamingHistogram]] derived from this [[Tile]],
+    * equalize the histogram of this tile.
+    *
+    * @params  histogram  The histogram of this tile
+    */
+  def equalize(histogram: StreamingHistogram): Tile = HistogramEqualization(self, histogram)
+
+  /**
+    * Equalize the histogram of this [[Tile]].
+    */
+  def equalize(): Tile = HistogramEqualization(self)
+}

raster/src/main/scala/geotrellis/raster/histogram/FastMapHistogram.scala

@@ -310,6 +310,18 @@ class FastMapHistogram(_size: Int, _buckets: Array[Int], _counts: Array[Long], _
       None
   }
 
+  /**
+    * CDF of the distribution.
+    */
+  def cdf(): Array[(Double, Double)] = {
+    val pdf = counts.map({ n => (n / total.toDouble) })
+
+    buckets
+      .map({ label => label.toDouble })
+      .zip(pdf.scanLeft(0.0)(_ + _).drop(1))
+      .toArray
+  }
+
   /**
     * Return the smallest and largest values seen by the histogram, if
     * it has seen any values.

raster/src/main/scala/geotrellis/raster/histogram/Histogram.scala

@@ -44,6 +44,11 @@ abstract trait Histogram[@specialized (Int, Double) T <: AnyVal] extends Seriali
    */
   def maxValue(): Option[T]
 
+  /**
+    * CDF of the distribution.
+    */
+  def cdf(): Array[(Double, Double)]
+
   /**
    * Return the smallest and largest items seen as a tuple.
    */