From 63460f2c261c53182aedcaeb6f49fcc85803a33f Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Wed, 1 Feb 2017 16:28:57 -0500
Subject: [PATCH 01/28] Remove Excess Whitespace

Signed-off-by: James McClain <jmcclain@azavea.com>
---
 .../geotrellis/raster/costdistance/CostDistanceSpec.scala     | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
index 6b5db3febc..1a94fd612f 100644
--- a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
+++ b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
@@ -33,7 +33,7 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
 
   def asTile(a: Array[Int], cols: Int, rows: Int): Tile =
     IntArrayTile(a, cols, rows)
-  
+
   test("ESRI example") {
     val n = NODATA
     val N = Double.NaN
@@ -90,7 +90,7 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
     val cd = CostDistance(costTile, points)
 
     val d = cd.toArrayDouble()
-    
+
     val expected = Array(
       22,21,21,20,17,15,14,
       20,19,22,20,15,12,11,

From 2313427bc74e65e16f2e42fcc4252dd7ff1826b1 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Wed, 1 Feb 2017 16:57:27 -0500
Subject: [PATCH 02/28] Simplify CostDistance Implementation

---
 .../raster/costdistance/CostDistance.scala    | 288 ++++++------------
 1 file changed, 88 insertions(+), 200 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 867ab20313..a50638742b 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -27,225 +27,113 @@ import geotrellis.raster._
   */
 object CostDistance {
 
-  /**
-    * Generate a Cost-Distance raster based on a set of starting points
-    * and a cost raster
-    *
-    * @param cost     Cost Tile (Int)
-    * @param points   List of starting points as tuples
-    *
-    * @note    Operation will only work with integer typed Cost Tiles (BitCellType, ByteConstantNoDataCellType, ShortConstantNoDataCellType, IntConstantNoDataCellType).
-    *          If a double typed Cost Tile (FloatConstantNoDataCellType, DoubleConstantNoDataCellType) is passed in, those costs will be rounded
-    *          to their floor integer values.
-    *
-    */
-  def apply(cost: Tile, points: Seq[(Int, Int)]): Tile = {
-    val (cols, rows) = cost.dimensions
-    val output = DoubleArrayTile.empty(cols, rows)
+  type Cost = (Int, Int, Double, Double) // column, row, friction, cost
+  type Q = PriorityQueue[Cost]
 
-    for((c, r) <- points)
-      output.setDouble(c, r, 0.0)
+  val q: Q = new PriorityQueue(
+    (1<<11), new java.util.Comparator[Cost] {
+      override def equals(a: Any) = a.equals(this)
+      def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
+    })
 
-    val pqueue = new PriorityQueue(
-        1000, new java.util.Comparator[Cost] {
-          override def equals(a: Any) = a.equals(this)
-          def compare(a: Cost, b: Cost) = a._3.compareTo(b._3)
-        })
+  var cols: Int = -1
+  var rows: Int = -1
 
+  def inTile(col: Int, row: Int): Boolean =
+    ((0 <= col && col < cols) && (0 <= row && row < rows))
 
-    for((c, r) <- points) {
-      calcNeighbors(c, r, cost, output, pqueue)
-    }
+  def isPassable(f: Double): Boolean =
+    (isData(f) && 0 <= f)
 
-    while (!pqueue.isEmpty) {
-      val (c, r, v) = pqueue.poll
-      if (v == output.getDouble(c, r)) calcNeighbors(c, r, cost, output, pqueue)
-    }
-
-    output
-  }
 
   /**
-    * Predicate: are the given column and row within the bounds of the
-    * [[Tile]]?
+    * Generate a cost-distance raster based on a set of starting
+    * points and a friction raster.  This is an implementation of the
+    * standard algorithm from [1].
+    *
+    * 1. Tomlin, Dana.
+    *    "Propagating radial waves of travel cost in a grid."
+    *    International Journal of Geographical Information Science 24.9 (2010): 1391-1413.
+    *
+    * @param  friction  Friction tile
+    * @param  points    List of starting points as tuples
     *
-    * @param  c  The column
-    * @param  r  The row
-    */
-  private def isValid(c: Int, r: Int, cost: Tile): Boolean =
-    c >= 0 && r >= 0 && c < cost.cols && r < cost.rows
-
-  /**
-    * A class encoding directions.
     */
-  final class Dir(val dc: Int, val dr: Int) {
-    val diag = !(dc == 0 || dr == 0)
-
-    lazy val cornerOffsets = (dc, dr) match {
-      case (-1, -1) => Array((0, -1), (-1, 0))
-      case ( 1, -1) => Array((0, -1), ( 1, 0))
-      case (-1, 1) => Array((0, 1), (-1, 0))
-      case ( 1, 1) => Array((0, 1), ( 1, 0))
-      case _ => Array[(Int, Int)]()
-    }
-
-    def apply(c: Int, r: Int) = (c + dc, r + dr)
-
-    lazy val unitDistance = if (diag) 1.41421356237 else 1.0
+  def apply(frictionTile: Tile, points: Seq[(Int, Int)]): Tile = {
+    val dims = frictionTile.dimensions
+    cols = dims._1; rows = dims._2
+    val costTile = DoubleArrayTile.empty(cols, rows)
+
+    points.foreach({ case (col, row) =>
+      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
+      q.add(entry)
+    })
+    while (!q.isEmpty) processNext(frictionTile, costTile, q)
+
+    costTile
   }
 
-  val dirs: Array[Dir] = Array(
-    (-1, -1), ( 0, -1), ( 1, -1),
-    (-1, 0),          ( 1, 0),
-    (-1, 1), ( 0, 1), ( 1, 1)).map { case (c, r) => new Dir(c, r) }
-
-
   /**
-    * Input:
-    * (c, r) => Source cell
-    * (dc, dr) => Delta (direction)
-    * cost => Cost raster
-    * d => C - D output raster
+    * Given a location, an instantaneous cost at that neighboring
+    * location (friction), a friction tile, the cost to get to that
+    * location, and the distance from the neighboring pixel, enqueue a
+    * candidate path to this pixel.
     *
-    * Output:
-    * List((c, r)) <- list of cells set
+    * @param  col           The column of the given location
+    * @param  row           The row of the given location
+    * @param  friction1     The instantaneous cost (friction) at the neighboring location
+    * @param  frictionTile  The friction tile
+    * @param  cost          The length of the best-known path from a source to the neighboring location
+    * @param  distance      The distance from the neighboring location to this location
     */
-  private def calcCostCell(c: Int, r: Int, dir: Dir, cost: Tile, d: DoubleArrayTile) = {
-    val cr = dir(c, r)
-
-    if (isValid(cr._1, cr._2, cost)) {
-      val prev = d.getDouble(cr._1, cr._2)
-      if (prev == 0.0) { // This is a source cell, don't override and shortcircuit early
-        None
-      } else {
-        val source = d.getDouble(c, r)
-
-        // Previous value could be NODATA
-        val prevCost = if (isNoData(prev)) java.lang.Double.MAX_VALUE else prev
-
-        var curMinCost = Double.MaxValue
-
-        // Simple solution
-        val baseCostOpt = calcCost(c, r, dir, cost)
-        if (baseCostOpt.isDefined) {
-          curMinCost = source + baseCostOpt.get
-        }
-
-        // Alternative solutions (going around the corner)
-        // Generally you can check diags directly:
-        // +---+---+---+
-        // | a | b | c |
-        // +---+---+---+
-        // | d | e | f |
-        // +---+---+---+
-        // | g | h | i |
-        // +---+---+---+
-        //
-        // For instance, "eg" can be computed directly
-        // but it turns out "eg" can be more expensive
-        // than "edg" or "ehg" so we compute those right
-        // now just in case
-        val cornerOffsets = dir.cornerOffsets
-        val l = cornerOffsets.length
-        var z = 0
-        while(z < l) {
-          val p = cornerOffsets(z)
-          val c1 = p._1 + c
-          val r1 = p._2 + r
-          val cost1 = calcCost(c, r, c1, r1, cost)
-          if (cost1.isDefined) {
-            val cost2 = calcCost(c1, r1, dir(c, r), cost)
-            if (cost2.isDefined) {
-              curMinCost = math.min(curMinCost, source + cost1.get + cost2.get)
-            }
-          }
-          z += 1
-        }
-
-        if (curMinCost == Double.MaxValue) {
-          None // Possible all nodata values
-        } else {
-          if (curMinCost < prevCost) {
-            d.setDouble(cr._1, cr._2, curMinCost)
-
-            Some((cr._1, cr._2, curMinCost))
-          } else {
-            None
-          }
-        }
+  @inline private def enqueueNeighbor(
+    col: Int, row: Int, friction1: Double,
+    frictionTile: Tile, cost: Double,
+    distance: Double = 1.0
+  ): Unit = {
+    // If the location is inside of the tile ...
+    if (inTile(col, row)) {
+      val friction2 = frictionTile.getDouble(col, row)
+      // ... and if the location is passable, enqueue it for future processing
+      if (isPassable(friction2)) {
+        val entry = (col, row, friction2, cost + distance * (friction1 + friction2) / 2.0)
+        q.add(entry)
       }
-    } else {
-      None
     }
   }
 
-  type Cost = (Int, Int, Double)
-
-  private def calcNeighbors(c: Int, r: Int, cost: Tile, d: DoubleArrayTile, p: PriorityQueue[Cost]) {
-    val l = dirs.length
-    var z = 0
-
-    while(z < l) {
-      val opt = calcCostCell(c, r, dirs(z), cost, d)
-      if (opt.isDefined) {
-        p.add(opt.get)
+  /**
+    * Process the candidate path on the top of the queue.
+    *
+    * @param  frictionTile  The friction tile
+    * @param  costTile      The cost tile
+    * @param  q             The priority queue of candidate paths
+    */
+  private def processNext(frictionTile: Tile, costTile: DoubleArrayTile, q: Q): Unit = {
+    val (col, row, friction1, candidateCost) = q.poll
+    val currentCost =
+      if (inTile(col, row))
+        costTile.getDouble(col, row)
+      else
+        Double.NaN
+
+    // If the candidate path is a possible improvement ...
+    if (!isData(currentCost) || candidateCost <= currentCost) {
+      // ... over-write the current cost with the candidate cost
+      if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
+
+      // ... compute candidate costs for neighbors and enqueue them
+      if (isPassable(friction1)) {
+
+        enqueueNeighbor(col-1, row+0, friction1, frictionTile, candidateCost)
+        enqueueNeighbor(col+1, row+0, friction1, frictionTile, candidateCost)
+        enqueueNeighbor(col+0, row+1, friction1, frictionTile, candidateCost)
+        enqueueNeighbor(col+0, row-1, friction1, frictionTile, candidateCost)
+        enqueueNeighbor(col-1, row-1, friction1, frictionTile, candidateCost, math.sqrt(2))
+        enqueueNeighbor(col-1, row+1, friction1, frictionTile, candidateCost, math.sqrt(2))
+        enqueueNeighbor(col+1, row-1, friction1, frictionTile, candidateCost, math.sqrt(2))
+        enqueueNeighbor(col+1, row+1, friction1, frictionTile, candidateCost, math.sqrt(2))
       }
-      z += 1
     }
   }
-
-  private def factor(c: Int, r: Int, c1: Int, r1: Int) = if (c == c1 || r == r1) 1.0 else 1.41421356237
-
-  private def safeGet(c: Int, r: Int, cost: Tile): IOption = IOption(cost.get(c, r))
-
-  private def calcCost(c: Int, r: Int, c2: Int, r2: Int, cost: Tile): DOption = {
-    val cost1 = safeGet(c, r, cost)
-    val cost2 = safeGet(c2, r2, cost)
-
-    if (cost1.isDefined && cost2.isDefined) {
-      DOption(factor(c, r, c2, r2) * (cost1.get + cost2.get) / 2.0)
-    } else {
-      DOption.None
-    }
-  }
-
-  private def calcCost(c: Int, r: Int, cr2: (Int, Int), cost: Tile): DOption =
-    calcCost(c, r, cr2._1, cr2._2, cost)
-
-  private def calcCost(c: Int, r: Int, dir: Dir, cost: Tile): DOption =
-    calcCost(c, r, dir(c, r), cost)
-}
-
-/**
-  * Represents an optional integer using 'Tile NODATA' as a flag
-  */
-private [costdistance]
-class IOption(val v: Int) extends AnyVal {
-  def map(f: Int => Int) = if (isDefined) new IOption(f(v)) else this
-  def flatMap(f: Int => IOption) = if (isDefined) f(v) else this
-  def isDefined = isData(v)
-  def get = if (isDefined) v else sys.error("Get called on NODATA")
-}
-
-/**
-  * Represents an optional integer using 'Double.NaN' as a flag
-  */
-private [costdistance]
-class DOption(val v: Double) extends AnyVal {
-  def map(f: Double => Double) = if (isDefined) new DOption(f(v)) else this
-  def flatMap(f: Double => DOption) = if (isDefined) f(v) else this
-  def isDefined = isData(v)
-  def get = if (isDefined) v else sys.error("Get called on NaN")
-}
-
-private [costdistance]
-object IOption {
-  val None = new IOption(NODATA)
-  def apply(v: Int) = new IOption(v)
-}
-
-private [costdistance]
-object DOption {
-  val None = new DOption(Double.NaN)
-  def apply(v: Double) = new DOption(v)
 }

From 3f1320fdce40b018edf7cdfb2b4fed5625250b6a Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Wed, 1 Feb 2017 19:45:52 -0500
Subject: [PATCH 03/28] Bump Initial Priority Queue Size

Increase size by one binary order of magnitude.  Still proportional to
the square root of the area of the typical expected tile.
---
 .../geotrellis/raster/costdistance/CostDistance.scala  | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index a50638742b..16da9982db 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -31,7 +31,7 @@ object CostDistance {
   type Q = PriorityQueue[Cost]
 
   val q: Q = new PriorityQueue(
-    (1<<11), new java.util.Comparator[Cost] {
+    (1<<12), new java.util.Comparator[Cost] {
       override def equals(a: Any) = a.equals(this)
       def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
     })
@@ -75,9 +75,9 @@ object CostDistance {
 
   /**
     * Given a location, an instantaneous cost at that neighboring
-    * location (friction), a friction tile, the cost to get to that
-    * location, and the distance from the neighboring pixel, enqueue a
-    * candidate path to this pixel.
+    * location (friction), a friction tile, the cost to get to the
+    * neighboring location, and the distance from the neighboring
+    * pixel, enqueue a candidate path to the present pixel.
     *
     * @param  col           The column of the given location
     * @param  row           The row of the given location
@@ -119,12 +119,12 @@ object CostDistance {
 
     // If the candidate path is a possible improvement ...
     if (!isData(currentCost) || candidateCost <= currentCost) {
+
       // ... over-write the current cost with the candidate cost
       if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
 
       // ... compute candidate costs for neighbors and enqueue them
       if (isPassable(friction1)) {
-
         enqueueNeighbor(col-1, row+0, friction1, frictionTile, candidateCost)
         enqueueNeighbor(col+1, row+0, friction1, frictionTile, candidateCost)
         enqueueNeighbor(col+0, row+1, friction1, frictionTile, candidateCost)

From cb4a753efa76ca0c7bb0e51952eb163fd3fc7d49 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 2 Feb 2017 09:06:19 -0500
Subject: [PATCH 04/28] Make Thread Safe

---
 .../raster/costdistance/CostDistance.scala    | 160 +++++++++---------
 1 file changed, 77 insertions(+), 83 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 16da9982db..07debe8bb4 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -30,22 +30,6 @@ object CostDistance {
   type Cost = (Int, Int, Double, Double) // column, row, friction, cost
   type Q = PriorityQueue[Cost]
 
-  val q: Q = new PriorityQueue(
-    (1<<12), new java.util.Comparator[Cost] {
-      override def equals(a: Any) = a.equals(this)
-      def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
-    })
-
-  var cols: Int = -1
-  var rows: Int = -1
-
-  def inTile(col: Int, row: Int): Boolean =
-    ((0 <= col && col < cols) && (0 <= row && row < rows))
-
-  def isPassable(f: Double): Boolean =
-    (isData(f) && 0 <= f)
-
-
   /**
     * Generate a cost-distance raster based on a set of starting
     * points and a friction raster.  This is an implementation of the
@@ -60,80 +44,90 @@ object CostDistance {
     *
     */
   def apply(frictionTile: Tile, points: Seq[(Int, Int)]): Tile = {
-    val dims = frictionTile.dimensions
-    cols = dims._1; rows = dims._2
+    val cols = frictionTile.cols
+    val rows = frictionTile.rows
     val costTile = DoubleArrayTile.empty(cols, rows)
+    val q: Q = new PriorityQueue(
+      (cols*2 + rows*2), new java.util.Comparator[Cost] {
+        override def equals(a: Any) = a.equals(this)
+        def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
+      })
+
+    def inTile(col: Int, row: Int): Boolean =
+      ((0 <= col && col < cols) && (0 <= row && row < rows))
+
+    def isPassable(f: Double): Boolean =
+      (isData(f) && 0 <= f)
+
+    /**
+      * Given a location, an instantaneous cost at that neighboring
+      * location (friction), the cost to get to the neighboring
+      * location, and the distance from the neighboring pixel to this
+      * pixel, enqueue a candidate path to the present pixel.
+      *
+      * @param  col           The column of the given location
+      * @param  row           The row of the given location
+      * @param  friction1     The instantaneous cost (friction) at the neighboring location
+      * @param  cost          The length of the best-known path from a source to the neighboring location
+      * @param  distance      The distance from the neighboring location to this location
+      */
+    @inline def enqueueNeighbor(
+      col: Int, row: Int, friction1: Double,
+      cost: Double,
+      distance: Double = 1.0
+    ): Unit = {
+      // If the location is inside of the tile ...
+      if (inTile(col, row)) {
+        val friction2 = frictionTile.getDouble(col, row)
+        // ... and if the location is passable, enqueue it for future processing
+        if (isPassable(friction2)) {
+          val entry = (col, row, friction2, cost + distance * (friction1 + friction2) / 2.0)
+          q.add(entry)
+        }
+      }
+    }
+
+    /**
+      * Process the candidate path on the top of the queue.
+      *
+      * @param  frictionTile  The friction tile
+      * @param  costTile      The cost tile
+      * @param  q             The priority queue of candidate paths
+      */
+    def processNext(): Unit = {
+      val (col, row, friction1, candidateCost) = q.poll
+      val currentCost =
+        if (inTile(col, row))
+          costTile.getDouble(col, row)
+        else
+          Double.NaN
+
+      // If the candidate path is a possible improvement ...
+      if (!isData(currentCost) || candidateCost <= currentCost) {
+
+        // ... over-write the current cost with the candidate cost
+        if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
+
+        // ... compute candidate costs for neighbors and enqueue them
+        if (isPassable(friction1)) {
+          enqueueNeighbor(col-1, row+0, friction1, candidateCost)
+          enqueueNeighbor(col+1, row+0, friction1, candidateCost)
+          enqueueNeighbor(col+0, row+1, friction1, candidateCost)
+          enqueueNeighbor(col+0, row-1, friction1, candidateCost)
+          enqueueNeighbor(col-1, row-1, friction1, candidateCost, math.sqrt(2))
+          enqueueNeighbor(col-1, row+1, friction1, candidateCost, math.sqrt(2))
+          enqueueNeighbor(col+1, row-1, friction1, candidateCost, math.sqrt(2))
+          enqueueNeighbor(col+1, row+1, friction1, candidateCost, math.sqrt(2))
+        }
+      }
+    }
 
     points.foreach({ case (col, row) =>
       val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
       q.add(entry)
     })
-    while (!q.isEmpty) processNext(frictionTile, costTile, q)
+    while (!q.isEmpty) processNext
 
     costTile
   }
-
-  /**
-    * Given a location, an instantaneous cost at that neighboring
-    * location (friction), a friction tile, the cost to get to the
-    * neighboring location, and the distance from the neighboring
-    * pixel, enqueue a candidate path to the present pixel.
-    *
-    * @param  col           The column of the given location
-    * @param  row           The row of the given location
-    * @param  friction1     The instantaneous cost (friction) at the neighboring location
-    * @param  frictionTile  The friction tile
-    * @param  cost          The length of the best-known path from a source to the neighboring location
-    * @param  distance      The distance from the neighboring location to this location
-    */
-  @inline private def enqueueNeighbor(
-    col: Int, row: Int, friction1: Double,
-    frictionTile: Tile, cost: Double,
-    distance: Double = 1.0
-  ): Unit = {
-    // If the location is inside of the tile ...
-    if (inTile(col, row)) {
-      val friction2 = frictionTile.getDouble(col, row)
-      // ... and if the location is passable, enqueue it for future processing
-      if (isPassable(friction2)) {
-        val entry = (col, row, friction2, cost + distance * (friction1 + friction2) / 2.0)
-        q.add(entry)
-      }
-    }
-  }
-
-  /**
-    * Process the candidate path on the top of the queue.
-    *
-    * @param  frictionTile  The friction tile
-    * @param  costTile      The cost tile
-    * @param  q             The priority queue of candidate paths
-    */
-  private def processNext(frictionTile: Tile, costTile: DoubleArrayTile, q: Q): Unit = {
-    val (col, row, friction1, candidateCost) = q.poll
-    val currentCost =
-      if (inTile(col, row))
-        costTile.getDouble(col, row)
-      else
-        Double.NaN
-
-    // If the candidate path is a possible improvement ...
-    if (!isData(currentCost) || candidateCost <= currentCost) {
-
-      // ... over-write the current cost with the candidate cost
-      if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
-
-      // ... compute candidate costs for neighbors and enqueue them
-      if (isPassable(friction1)) {
-        enqueueNeighbor(col-1, row+0, friction1, frictionTile, candidateCost)
-        enqueueNeighbor(col+1, row+0, friction1, frictionTile, candidateCost)
-        enqueueNeighbor(col+0, row+1, friction1, frictionTile, candidateCost)
-        enqueueNeighbor(col+0, row-1, friction1, frictionTile, candidateCost)
-        enqueueNeighbor(col-1, row-1, friction1, frictionTile, candidateCost, math.sqrt(2))
-        enqueueNeighbor(col-1, row+1, friction1, frictionTile, candidateCost, math.sqrt(2))
-        enqueueNeighbor(col+1, row-1, friction1, frictionTile, candidateCost, math.sqrt(2))
-        enqueueNeighbor(col+1, row+1, friction1, frictionTile, candidateCost, math.sqrt(2))
-      }
-    }
-  }
 }

From 177bee2ae1f83e8c5b085bc09fb8d96b2cf90951 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 2 Feb 2017 09:49:32 -0500
Subject: [PATCH 05/28] Add Capability To Track Edge Changes

---
 .../raster/costdistance/CostDistance.scala    | 66 ++++++++++++++-----
 1 file changed, 51 insertions(+), 15 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 07debe8bb4..ddfb9da1dd 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -29,6 +29,27 @@ object CostDistance {
 
   type Cost = (Int, Int, Double, Double) // column, row, friction, cost
   type Q = PriorityQueue[Cost]
+  type EdgeCallback = (Cost => Unit)
+
+  def apply(frictionTile: Tile, points: Seq[(Int, Int)]): Tile = {
+    val cols = frictionTile.cols
+    val rows = frictionTile.rows
+    val costTile = DoubleArrayTile.empty(cols, rows)
+    val q: Q = new PriorityQueue(
+      (cols*2 + rows*2), new java.util.Comparator[Cost] {
+        override def equals(a: Any) = a.equals(this)
+        def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
+      })
+
+    def nop(cost: Cost): Unit = {}
+
+    points.foreach({ case (col, row) =>
+      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
+      q.add(entry)
+    })
+
+    compute(frictionTile, costTile, q, nop, nop, nop, nop)
+  }
 
   /**
     * Generate a cost-distance raster based on a set of starting
@@ -43,15 +64,17 @@ object CostDistance {
     * @param  points    List of starting points as tuples
     *
     */
-  def apply(frictionTile: Tile, points: Seq[(Int, Int)]): Tile = {
+  def compute(
+    frictionTile: Tile,
+    costTile: DoubleArrayTile,
+    q: Q,
+    leftCallback: EdgeCallback, rightCallback: EdgeCallback,
+    topCallback: EdgeCallback, bottomCallback: EdgeCallback
+  ): Tile = {
     val cols = frictionTile.cols
     val rows = frictionTile.rows
-    val costTile = DoubleArrayTile.empty(cols, rows)
-    val q: Q = new PriorityQueue(
-      (cols*2 + rows*2), new java.util.Comparator[Cost] {
-        override def equals(a: Any) = a.equals(this)
-        def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
-      })
+
+    require(frictionTile.dimensions == costTile.dimensions)
 
     def inTile(col: Int, row: Int): Boolean =
       ((0 <= col && col < cols) && (0 <= row && row < rows))
@@ -59,6 +82,14 @@ object CostDistance {
     def isPassable(f: Double): Boolean =
       (isData(f) && 0 <= f)
 
+    def leftEdge(col: Int, row: Int): Boolean = (col == 0)
+
+    def rightEdge(col: Int, row: Int): Boolean = (col == cols-1)
+
+    def topEdge(col: Int, row: Int): Boolean = (row == 0)
+
+    def bottomEdge(col: Int, row: Int): Boolean = (row == rows-1)
+
     /**
       * Given a location, an instantaneous cost at that neighboring
       * location (friction), the cost to get to the neighboring
@@ -95,7 +126,8 @@ object CostDistance {
       * @param  q             The priority queue of candidate paths
       */
     def processNext(): Unit = {
-      val (col, row, friction1, candidateCost) = q.poll
+      val cost: Cost = q.poll
+      val (col, row, friction1, candidateCost) = cost
       val currentCost =
         if (inTile(col, row))
           costTile.getDouble(col, row)
@@ -105,10 +137,18 @@ object CostDistance {
       // If the candidate path is a possible improvement ...
       if (!isData(currentCost) || candidateCost <= currentCost) {
 
-        // ... over-write the current cost with the candidate cost
-        if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
+        // Over-write the current cost with the candidate cost
+        if (inTile(col, row)) {
+          costTile.setDouble(col, row, candidateCost)
+
+          // Register changes on the boundary
+          if (leftEdge(col, row)) leftCallback(cost)
+          if (rightEdge(col, row)) rightCallback(cost)
+          if (topEdge(col, row)) topCallback(cost)
+          if (bottomEdge(col, row)) bottomCallback(cost)
+        }
 
-        // ... compute candidate costs for neighbors and enqueue them
+        // Compute candidate costs for neighbors and enqueue them
         if (isPassable(friction1)) {
           enqueueNeighbor(col-1, row+0, friction1, candidateCost)
           enqueueNeighbor(col+1, row+0, friction1, candidateCost)
@@ -122,10 +162,6 @@ object CostDistance {
       }
     }
 
-    points.foreach({ case (col, row) =>
-      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
-      q.add(entry)
-    })
     while (!q.isEmpty) processNext
 
     costTile

From fe930da2af3394d0dd0b0ccf207908cc0d4e0f85 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 2 Feb 2017 10:28:59 -0500
Subject: [PATCH 06/28] Add Maximum-Cost Capability

---
 .../raster/costdistance/CostDistance.scala    | 59 +++++++++++++------
 1 file changed, 42 insertions(+), 17 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index ddfb9da1dd..3abf528067 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -31,24 +31,16 @@ object CostDistance {
   type Q = PriorityQueue[Cost]
   type EdgeCallback = (Cost => Unit)
 
-  def apply(frictionTile: Tile, points: Seq[(Int, Int)]): Tile = {
-    val cols = frictionTile.cols
-    val rows = frictionTile.rows
-    val costTile = DoubleArrayTile.empty(cols, rows)
-    val q: Q = new PriorityQueue(
-      (cols*2 + rows*2), new java.util.Comparator[Cost] {
+  /**
+    * Generate a Queue suitable for working with a tile of the given
+    * dimensions.
+    */
+  def generateQueue(cols: Int, rows: Int): Q = {
+    new PriorityQueue(
+      (cols*16 + rows*16), new java.util.Comparator[Cost] {
         override def equals(a: Any) = a.equals(this)
         def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
       })
-
-    def nop(cost: Cost): Unit = {}
-
-    points.foreach({ case (col, row) =>
-      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
-      q.add(entry)
-    })
-
-    compute(frictionTile, costTile, q, nop, nop, nop, nop)
   }
 
   /**
@@ -64,13 +56,46 @@ object CostDistance {
     * @param  points    List of starting points as tuples
     *
     */
+  def apply(
+    frictionTile: Tile,
+    points: Seq[(Int, Int)],
+    maxCost: Double = Double.PositiveInfinity
+  ): DoubleArrayTile = {
+    val cols = frictionTile.cols
+    val rows = frictionTile.rows
+    val costTile = DoubleArrayTile.empty(cols, rows)
+    val q: Q = generateQueue(cols, rows)
+
+    def nop(cost: Cost): Unit = {}
+
+    points.foreach({ case (col, row) =>
+      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
+      q.add(entry)
+    })
+
+    compute(frictionTile, costTile, maxCost, q, nop, nop, nop, nop)
+  }
+
+  /**
+    * Compute a cost tile.
+    *
+    * @param  frictionTile    The friction tile
+    * @param  costTile        The tile that will contain the costs
+    * @param  maxCost         The maximum cost of any path (truncates to limit computational cost)
+    * @param  q               A priority queue of Cost objects (a.k.a. candidate paths)
+    * @param  leftCallback    Called when a pixel in the left-most column is updated
+    * @param  rightCallback   Called when a pixel in the right-most column is updated
+    * @param  topCallbck      Called when a pixel in the top-most row is updated
+    * @param  bottomCallback  Called when a pixel in the bottom-most row is updated
+    */
   def compute(
     frictionTile: Tile,
     costTile: DoubleArrayTile,
+    maxCost: Double,
     q: Q,
     leftCallback: EdgeCallback, rightCallback: EdgeCallback,
     topCallback: EdgeCallback, bottomCallback: EdgeCallback
-  ): Tile = {
+  ): DoubleArrayTile = {
     val cols = frictionTile.cols
     val rows = frictionTile.rows
 
@@ -135,7 +160,7 @@ object CostDistance {
           Double.NaN
 
       // If the candidate path is a possible improvement ...
-      if (!isData(currentCost) || candidateCost <= currentCost) {
+      if ((candidateCost <= maxCost) && (!isData(currentCost) || candidateCost <= currentCost)) {
 
         // Over-write the current cost with the candidate cost
         if (inTile(col, row)) {

From a406228dae6779a93d65ca34d104724e4cf5d637 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 2 Feb 2017 12:58:33 -0500
Subject: [PATCH 07/28] Add Tests for Boundary Callbacks

---
 .../costdistance/CostDistanceSpec.scala       | 55 +++++++++++++++++++
 .../raster/costdistance/CostDistance.scala    | 17 ++++--
 2 files changed, 66 insertions(+), 6 deletions(-)

diff --git a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
index 1a94fd612f..02878d4803 100644
--- a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
+++ b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
@@ -108,6 +108,61 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
     }
   }
 
+  test("Boundary Callbacks") {
+    val n = NODATA
+    val N = Double.NaN
+    val size = 6
+
+    // Example from ESRI
+    // http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Cost_Distance_algorithm
+    val frictionTile = Array(
+      1,3,4,4,3,2,
+      4,6,2,3,7,6,
+      5,8,7,5,6,6,
+      1,4,5,n,5,1,
+      4,7,5,n,2,6,
+      1,2,2,1,3,4)
+
+    val actualLeft = Array.ofDim[Double](size)
+    val actualRight = Array.ofDim[Double](size)
+    val actualTop = Array.ofDim[Double](size)
+    val actualBottom = Array.ofDim[Double](size)
+
+    val expectedLeft = Array(2.0, 4.5, 8.0, 5.0, 2.5, 0.0)
+    val expectedRight = Array(9.2, 13.1, 12.7, 9.2, 11.1, 10.5)
+    val expectedTop = Array(2.0, 0.0, 0.0, 4.0, 6.7, 9.2)
+    val expectedBottom = Array(0.0, 1.5, 3.5, 5.0, 7.0, 10.5)
+
+    // Produce priority queue
+    val q = CostDistance.generateEmptyQueue(size, size)
+
+    // Insert starting points
+    val points = Seq((1,0), (2,0), (2,1), (0,5))
+    points.foreach({ case (col, row) =>
+      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
+      q.add(entry)
+    })
+
+    // Generate initial (empty) cost tile
+    val costTile = CostDistance.generateEmptyCostTile(size, size)
+
+    // Various callbacks to fill in the "actual" arrays
+    val leftCb: CostDistance.EdgeCallback = { case (_, row: Int, _, cost: Double) => actualLeft(row) = math.floor(cost*10)/10 }
+    val topCb: CostDistance.EdgeCallback = { case (col: Int, _, _, cost: Double) => actualTop(col) = math.floor(cost*10)/10 }
+    val rightCb: CostDistance.EdgeCallback = { case (_, row: Int, _, cost: Double) => actualRight(row) = math.floor(cost*10)/10 }
+    val bottomCb: CostDistance.EdgeCallback = { case (col: Int, _, _, cost: Double) => actualBottom(col) = math.floor(cost*10)/10 }
+
+    CostDistance.compute(
+      frictionTile, costTile,
+      Double.PositiveInfinity, q,
+      leftCb, topCb, rightCb, bottomCb)
+
+    actualLeft should be (expectedLeft)
+    actualTop should be (expectedTop)
+    actualRight should be (expectedRight)
+    actualBottom should be (expectedBottom)
+  }
+
   def print(d: DoubleArrayTile):Unit = println(d.array.toList.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i)).grouped(d.cols).map(_.mkString(",")).mkString("\n"))
 
 }
diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 3abf528067..8d3f6a60cf 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -35,7 +35,7 @@ object CostDistance {
     * Generate a Queue suitable for working with a tile of the given
     * dimensions.
     */
-  def generateQueue(cols: Int, rows: Int): Q = {
+  def generateEmptyQueue(cols: Int, rows: Int): Q = {
     new PriorityQueue(
       (cols*16 + rows*16), new java.util.Comparator[Cost] {
         override def equals(a: Any) = a.equals(this)
@@ -43,6 +43,9 @@ object CostDistance {
       })
   }
 
+  def generateEmptyCostTile(cols: Int, rows: Int): DoubleArrayTile =
+    DoubleArrayTile.empty(cols, rows)
+
   /**
     * Generate a cost-distance raster based on a set of starting
     * points and a friction raster.  This is an implementation of the
@@ -63,8 +66,8 @@ object CostDistance {
   ): DoubleArrayTile = {
     val cols = frictionTile.cols
     val rows = frictionTile.rows
-    val costTile = DoubleArrayTile.empty(cols, rows)
-    val q: Q = generateQueue(cols, rows)
+    val costTile = generateEmptyCostTile(cols, rows)
+    val q: Q = generateEmptyQueue(cols, rows)
 
     def nop(cost: Cost): Unit = {}
 
@@ -84,8 +87,8 @@ object CostDistance {
     * @param  maxCost         The maximum cost of any path (truncates to limit computational cost)
     * @param  q               A priority queue of Cost objects (a.k.a. candidate paths)
     * @param  leftCallback    Called when a pixel in the left-most column is updated
-    * @param  rightCallback   Called when a pixel in the right-most column is updated
     * @param  topCallbck      Called when a pixel in the top-most row is updated
+    * @param  rightCallback   Called when a pixel in the right-most column is updated
     * @param  bottomCallback  Called when a pixel in the bottom-most row is updated
     */
   def compute(
@@ -93,8 +96,10 @@ object CostDistance {
     costTile: DoubleArrayTile,
     maxCost: Double,
     q: Q,
-    leftCallback: EdgeCallback, rightCallback: EdgeCallback,
-    topCallback: EdgeCallback, bottomCallback: EdgeCallback
+    leftCallback: EdgeCallback,
+    topCallback: EdgeCallback,
+    rightCallback: EdgeCallback,
+    bottomCallback: EdgeCallback
   ): DoubleArrayTile = {
     val cols = frictionTile.cols
     val rows = frictionTile.rows

From 921902446a33f9e2daf990406d0ae4f61e9007bd Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 2 Feb 2017 13:12:44 -0500
Subject: [PATCH 08/28] Add Max Distance Tests

---
 .../costdistance/CostDistanceSpec.scala       | 44 ++++++++++++++++++-
 1 file changed, 42 insertions(+), 2 deletions(-)

diff --git a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
index 02878d4803..c08cb24352 100644
--- a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
+++ b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
@@ -16,14 +16,15 @@
 
 package geotrellis.raster.costdistance
 
-import java.util.Locale
-
 import geotrellis.raster._
 import geotrellis.raster.testkit._
+
 import org.scalatest._
 
+import java.util.Locale
 import scala.language.implicitConversions
 
+
 class CostDistanceSpec extends FunSuite with RasterMatchers {
   implicit def array2Tile(a: Array[Int]): Tile = {
     val size = math.sqrt(a.length).toInt
@@ -163,6 +164,45 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
     actualBottom should be (expectedBottom)
   }
 
+  test("Max Distance") {
+    val n = NODATA
+    val N = Double.NaN
+
+    // Example from ESRI
+    // http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Cost_Distance_algorithm
+    val costTile = Array(
+      1,3,4,4,3,2,
+      4,6,2,3,7,6,
+      5,8,7,5,6,6,
+      1,4,5,n,5,1,
+      4,7,5,n,2,6,
+      1,2,2,1,3,4)
+
+    val points = Seq(
+      (1,0),
+      (2,0),
+      (2,1),
+      (0,5))
+
+    val cd = CostDistance(costTile, points, 2.6)
+
+    val d = cd.toArrayDouble()
+
+    val expected = Array(
+      2.0, 0.0, 0.0,   N,   N,   N,
+        N,   N, 0.0, 2.5,   N,   N,
+        N,   N,   N,   N,   N,   N,
+        N,   N,   N,   N,   N,   N,
+      2.5,   N,   N,   N,   N,   N,
+      0.0, 1.5,   N,   N,   N,   N).map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
+
+    val strings = d.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
+
+    for(i <- 0 until strings.length) {
+      strings(i) should be (expected(i))
+    }
+  }
+
   def print(d: DoubleArrayTile):Unit = println(d.array.toList.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i)).grouped(d.cols).map(_.mkString(",")).mkString("\n"))
 
 }

From c44f4d50a3d75997ccbea8c2cebb62b4d48b5456 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Fri, 3 Feb 2017 08:55:06 -0500
Subject: [PATCH 09/28] Reduce Number of Callbacks

---
 .../costdistance/CostDistanceSpec.scala       | 24 ++++++++-------
 .../raster/costdistance/CostDistance.scala    | 30 +++++--------------
 2 files changed, 21 insertions(+), 33 deletions(-)

diff --git a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
index c08cb24352..48c3e31a43 100644
--- a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
+++ b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
@@ -125,14 +125,14 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
       1,2,2,1,3,4)
 
     val actualLeft = Array.ofDim[Double](size)
+    val actualUp = Array.ofDim[Double](size)
     val actualRight = Array.ofDim[Double](size)
-    val actualTop = Array.ofDim[Double](size)
-    val actualBottom = Array.ofDim[Double](size)
+    val actualDown = Array.ofDim[Double](size)
 
     val expectedLeft = Array(2.0, 4.5, 8.0, 5.0, 2.5, 0.0)
+    val expectedUp = Array(2.0, 0.0, 0.0, 4.0, 6.7, 9.2)
     val expectedRight = Array(9.2, 13.1, 12.7, 9.2, 11.1, 10.5)
-    val expectedTop = Array(2.0, 0.0, 0.0, 4.0, 6.7, 9.2)
-    val expectedBottom = Array(0.0, 1.5, 3.5, 5.0, 7.0, 10.5)
+    val expectedDown = Array(0.0, 1.5, 3.5, 5.0, 7.0, 10.5)
 
     // Produce priority queue
     val q = CostDistance.generateEmptyQueue(size, size)
@@ -148,20 +148,22 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
     val costTile = CostDistance.generateEmptyCostTile(size, size)
 
     // Various callbacks to fill in the "actual" arrays
-    val leftCb: CostDistance.EdgeCallback = { case (_, row: Int, _, cost: Double) => actualLeft(row) = math.floor(cost*10)/10 }
-    val topCb: CostDistance.EdgeCallback = { case (col: Int, _, _, cost: Double) => actualTop(col) = math.floor(cost*10)/10 }
-    val rightCb: CostDistance.EdgeCallback = { case (_, row: Int, _, cost: Double) => actualRight(row) = math.floor(cost*10)/10 }
-    val bottomCb: CostDistance.EdgeCallback = { case (col: Int, _, _, cost: Double) => actualBottom(col) = math.floor(cost*10)/10 }
+    val edgeCb: CostDistance.EdgeCallback = { case (col: Int, row: Int, _: Double, cost: Double) =>
+      if (col == 0) actualLeft(row) = math.floor(cost*10)/10
+      if (col == size-1) actualRight(row) = math.floor(cost*10)/10
+      if (row == 0) actualUp(col) = math.floor(cost*10)/10
+      if (row == size-1) actualDown(col) = math.floor(cost*10)/10
+    }
 
     CostDistance.compute(
       frictionTile, costTile,
       Double.PositiveInfinity, q,
-      leftCb, topCb, rightCb, bottomCb)
+      edgeCb)
 
     actualLeft should be (expectedLeft)
-    actualTop should be (expectedTop)
+    actualUp should be (expectedUp)
     actualRight should be (expectedRight)
-    actualBottom should be (expectedBottom)
+    actualDown should be (expectedDown)
   }
 
   test("Max Distance") {
diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 8d3f6a60cf..44aeaaa788 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -16,10 +16,10 @@
 
 package geotrellis.raster.costdistance
 
-import java.util.PriorityQueue
-
 import geotrellis.raster._
 
+import java.util.PriorityQueue
+
 
 /**
   * Object housing various functions related to Cost-Distance
@@ -76,7 +76,7 @@ object CostDistance {
       q.add(entry)
     })
 
-    compute(frictionTile, costTile, maxCost, q, nop, nop, nop, nop)
+    compute(frictionTile, costTile, maxCost, q, nop)
   }
 
   /**
@@ -86,20 +86,14 @@ object CostDistance {
     * @param  costTile        The tile that will contain the costs
     * @param  maxCost         The maximum cost of any path (truncates to limit computational cost)
     * @param  q               A priority queue of Cost objects (a.k.a. candidate paths)
-    * @param  leftCallback    Called when a pixel in the left-most column is updated
-    * @param  topCallbck      Called when a pixel in the top-most row is updated
-    * @param  rightCallback   Called when a pixel in the right-most column is updated
-    * @param  bottomCallback  Called when a pixel in the bottom-most row is updated
+    * @param  edgeCallback    Called when a pixel on the edge of the tile is updated
     */
   def compute(
     frictionTile: Tile,
     costTile: DoubleArrayTile,
     maxCost: Double,
     q: Q,
-    leftCallback: EdgeCallback,
-    topCallback: EdgeCallback,
-    rightCallback: EdgeCallback,
-    bottomCallback: EdgeCallback
+    edgeCallback: EdgeCallback
   ): DoubleArrayTile = {
     val cols = frictionTile.cols
     val rows = frictionTile.rows
@@ -112,13 +106,8 @@ object CostDistance {
     def isPassable(f: Double): Boolean =
       (isData(f) && 0 <= f)
 
-    def leftEdge(col: Int, row: Int): Boolean = (col == 0)
-
-    def rightEdge(col: Int, row: Int): Boolean = (col == cols-1)
-
-    def topEdge(col: Int, row: Int): Boolean = (row == 0)
-
-    def bottomEdge(col: Int, row: Int): Boolean = (row == rows-1)
+    def onEdge(col: Int, row: Int): Boolean =
+      ((col == 0) || (row == 0) || (col == cols-1) || (row == rows-1))
 
     /**
       * Given a location, an instantaneous cost at that neighboring
@@ -172,10 +161,7 @@ object CostDistance {
           costTile.setDouble(col, row, candidateCost)
 
           // Register changes on the boundary
-          if (leftEdge(col, row)) leftCallback(cost)
-          if (rightEdge(col, row)) rightCallback(cost)
-          if (topEdge(col, row)) topCallback(cost)
-          if (bottomEdge(col, row)) bottomCallback(cost)
+          if (onEdge(col, row)) edgeCallback(cost)
         }
 
         // Compute candidate costs for neighbors and enqueue them

From 56c353c1759fa46de50a33a689973428b386de64 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 2 Feb 2017 14:23:37 -0500
Subject: [PATCH 10/28] Initial Foundation of MrGeo-inspired Cost-Distance

---
 .../costdistance/MrGeoCostDistance.scala      | 125 ++++++++++++++++++
 1 file changed, 125 insertions(+)
 create mode 100644 spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
new file mode 100644
index 0000000000..980f77d694
--- /dev/null
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
@@ -0,0 +1,125 @@
+package geotrellis.spark.costdistance
+
+import geotrellis.raster._
+import geotrellis.spark._
+import geotrellis.vector._
+import geotrellis.raster.costdistance.CostDistance
+
+import org.apache.spark.rdd.RDD
+import org.apache.spark.SparkContext
+import org.apache.spark.util.AccumulatorV2
+
+import scala.collection.mutable
+
+
+/**
+  * This Spark-enabled implementation of the standard cost-distance
+  * algorithm [1] is "heavily inspired" by the MrGeo implementation
+  * [2] but does not share any code with it.
+  *
+  * 1. Tomlin, Dana.
+  *    "Propagating radial waves of travel cost in a grid."
+  *    International Journal of Geographical Information Science 24.9 (2010): 1391-1413.
+  *
+  * 2. https://github.com/ngageoint/mrgeo/blob/0c6ed4a7e66bb0923ec5c570b102862aee9e885e/mrgeo-mapalgebra/mrgeo-mapalgebra-costdistance/src/main/scala/org/mrgeo/mapalgebra/CostDistanceMapOp.scala
+  */
+object MrGeoCostDistance {
+
+  type CostList = List[CostDistance.Cost]
+  type KeyListPair = (SpatialKey, CostList)
+  type PixelMap = mutable.Map[SpatialKey, CostList]
+
+  /**
+    * An accumulator to hold lists of edge changes.
+    */
+  class PixelMapAccumulator extends AccumulatorV2[KeyListPair, PixelMap] {
+    private var map: PixelMap = mutable.Map.empty
+
+    def copy: PixelMapAccumulator = {
+      val other = new PixelMapAccumulator
+      other.merge(this)
+      other
+    }
+    def add(kv: (SpatialKey, CostList)): Unit = {
+      val key = kv._1
+      val list = kv._2
+
+      if (!map.contains(key)) map.put(key, list)
+      else map.put(key, map.get(key).get ++ list)
+    }
+    def isZero: Boolean = map.isEmpty
+    def merge(other: AccumulatorV2[(SpatialKey, CostList), PixelMap]): Unit = map ++= other.value
+    def reset: Unit = map = mutable.Map.empty
+    def value: PixelMap = map
+  }
+
+  /**
+    * Perform the cost-distance computation.
+    */
+  def apply[K: (? => SpatialKey), V: (? => Tile)](
+    friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]],
+    points: List[Point],
+    maxCost: Double = Double.PositiveInfinity
+  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] = {
+
+    val accumulator = new PixelMapAccumulator
+    sc.register(accumulator, "Pixel Map Accumulator")
+
+    val mt = friction.metadata.mapTransform
+
+    // Load the accumulator with the starting values
+    friction.foreach({ case (k, v) =>
+      val key = implicitly[SpatialKey](k)
+      val tile = implicitly[Tile](v)
+      val cols = tile.cols
+      val rows = tile.rows
+      val extent = mt(key)
+      val rasterExtent = RasterExtent(extent, cols, rows)
+      val costs: List[CostDistance.Cost] = points
+        .filter({ point => extent.contains(point) })
+        .map({ point =>
+          val col = rasterExtent.mapXToGrid(point.x)
+          val row = rasterExtent.mapYToGrid(point.y)
+          val friction = tile.getDouble(col, row)
+          (col, row, friction, 0.0)
+        })
+
+      accumulator.add((key, costs))
+    })
+
+    // Move starting values from accumulator to local variable
+    var changes: PixelMap = accumulator.value
+    accumulator.reset
+
+    // Create RDD of initial (empty) cost tiles
+    var costs: RDD[(K, V, DoubleArrayTile)] = friction.map({ case (k, v) =>
+      val tile = implicitly[Tile](v)
+      val cols = tile.cols
+      val rows = tile.rows
+
+      (k, v, CostDistance.generateEmptyCostTile(cols, rows))
+    })
+
+    do {
+      costs = costs.map({ case (k, v, cost) =>
+        val key = implicitly[SpatialKey](k)
+        val friction = implicitly[Tile](v)
+        val cols = friction.cols
+        val rows = friction.rows
+        val q: CostDistance.Q = {
+          val q = CostDistance.generateEmptyQueue(cols, rows)
+
+          changes
+            .getOrElse(key, List.empty[CostDistance.Cost])
+            .foreach({ (entry: CostDistance.Cost) => q.add(entry) })
+          q
+        }
+
+        (k, v, CostDistance.compute(friction, cost, maxCost, q, CostDistance.nop)) })
+      changes = accumulator.value
+      accumulator.reset
+    } while (changes.size > 0)
+
+    costs.map({ case (k, _, cost) => (k, cost) })
+  }
+}

From 881f16a31120a9466860177aecc1525fa9c2d4dc Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Fri, 3 Feb 2017 20:26:10 -0500
Subject: [PATCH 11/28] Allow Serial Cost-Distance To Terminate

Force cost values up as quickly as possible to allow paths to be pruned.
---
 .../raster/costdistance/CostDistance.scala    | 26 ++++++++++++++-----
 1 file changed, 20 insertions(+), 6 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 44aeaaa788..36f1f8cb36 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -31,6 +31,11 @@ object CostDistance {
   type Q = PriorityQueue[Cost]
   type EdgeCallback = (Cost => Unit)
 
+  /**
+    * NOP EdgeCallback
+    */
+  def nop(cost: Cost): Unit = {}
+
   /**
     * Generate a Queue suitable for working with a tile of the given
     * dimensions.
@@ -69,8 +74,6 @@ object CostDistance {
     val costTile = generateEmptyCostTile(cols, rows)
     val q: Q = generateEmptyQueue(cols, rows)
 
-    def nop(cost: Cost): Unit = {}
-
     points.foreach({ case (col, row) =>
       val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
       q.add(entry)
@@ -129,10 +132,21 @@ object CostDistance {
       // If the location is inside of the tile ...
       if (inTile(col, row)) {
         val friction2 = frictionTile.getDouble(col, row)
-        // ... and if the location is passable, enqueue it for future processing
+        val currentCost = costTile.getDouble(col, row)
+
+        // ... and if the location is passable ...
         if (isPassable(friction2)) {
           val entry = (col, row, friction2, cost + distance * (friction1 + friction2) / 2.0)
-          q.add(entry)
+          val candidateCost = entry._4
+
+          // ... and the candidate cost is less than the maximum cost ...
+          if (candidateCost <= maxCost) {
+            // ... and the candidate is a possible improvement ...
+            if ((isData(currentCost) && candidateCost < currentCost) || !isData(currentCost)) {
+              costTile.setDouble(col, row, candidateCost) // then increase lower bound on pixel,
+              q.add(entry) // and enqueue candidate for future processing
+            }
+          }
         }
       }
     }
@@ -154,11 +168,11 @@ object CostDistance {
           Double.NaN
 
       // If the candidate path is a possible improvement ...
-      if ((candidateCost <= maxCost) && (!isData(currentCost) || candidateCost <= currentCost)) {
+      if (!isData(currentCost) || candidateCost <= currentCost) {
 
         // Over-write the current cost with the candidate cost
         if (inTile(col, row)) {
-          costTile.setDouble(col, row, candidateCost)
+          costTile.setDouble(col, row, candidateCost) // XXX
 
           // Register changes on the boundary
           if (onEdge(col, row)) edgeCallback(cost)

From 3d0f99625278a4c7b9edc11252c1bef58f276790 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Fri, 3 Feb 2017 15:39:59 -0500
Subject: [PATCH 12/28] First Iteration

This code implements the first iteration of the n-iteration distributed
cost-distance algorithm.
---
 .../costdistance/MrGeoCostDistance.scala      | 49 ++++++++++---------
 1 file changed, 25 insertions(+), 24 deletions(-)

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
index 980f77d694..52b08a03fd 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
@@ -60,7 +60,7 @@ object MrGeoCostDistance {
     friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]],
     points: List[Point],
     maxCost: Double = Double.PositiveInfinity
-  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] = {
+  )(implicit sc: SparkContext): RDD[(K, Tile)] = {
 
     val accumulator = new PixelMapAccumulator
     sc.register(accumulator, "Pixel Map Accumulator")
@@ -87,10 +87,6 @@ object MrGeoCostDistance {
       accumulator.add((key, costs))
     })
 
-    // Move starting values from accumulator to local variable
-    var changes: PixelMap = accumulator.value
-    accumulator.reset
-
     // Create RDD of initial (empty) cost tiles
     var costs: RDD[(K, V, DoubleArrayTile)] = friction.map({ case (k, v) =>
       val tile = implicitly[Tile](v)
@@ -100,25 +96,30 @@ object MrGeoCostDistance {
       (k, v, CostDistance.generateEmptyCostTile(cols, rows))
     })
 
-    do {
-      costs = costs.map({ case (k, v, cost) =>
-        val key = implicitly[SpatialKey](k)
-        val friction = implicitly[Tile](v)
-        val cols = friction.cols
-        val rows = friction.rows
-        val q: CostDistance.Q = {
-          val q = CostDistance.generateEmptyQueue(cols, rows)
-
-          changes
-            .getOrElse(key, List.empty[CostDistance.Cost])
-            .foreach({ (entry: CostDistance.Cost) => q.add(entry) })
-          q
-        }
-
-        (k, v, CostDistance.compute(friction, cost, maxCost, q, CostDistance.nop)) })
-      changes = accumulator.value
-      accumulator.reset
-    } while (changes.size > 0)
+    // XXX Should be shared variable
+    val changes = accumulator.value
+
+    // do {
+    costs = costs.map({ case (k, v, cost) =>
+      val key = implicitly[SpatialKey](k)
+      val friction = implicitly[Tile](v)
+      val cols = friction.cols
+      val rows = friction.rows
+      val q: CostDistance.Q = {
+        val q = CostDistance.generateEmptyQueue(cols, rows)
+
+        changes
+          .getOrElse(key, List.empty[CostDistance.Cost])
+          .foreach({ (entry: CostDistance.Cost) => q.add(entry) })
+        q
+      }
+
+      (k, v, CostDistance.compute(friction, cost, maxCost, q, CostDistance.nop))
+      // (k, v, CostDistance.compute(CostDistance.generateEmptyCostTile(cols, rows), cost, maxCost, q, CostDistance.nop))
+    })
+
+    accumulator.reset
+    // } while (changes.size > 0)
 
     costs.map({ case (k, _, cost) => (k, cost) })
   }

From f3328d0a2be91fd1a8c6106f44bfff60cb3a1e35 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Sun, 5 Feb 2017 14:37:54 -0500
Subject: [PATCH 13/28] nth Iteration

05 Feb 19:44:24 INFO [cdistance.CostDistance$] - MILLIS: 155272
05 Feb 19:48:29 INFO [cdistance.CostDistance$] - MILLIS: 153500
---
 .../raster/costdistance/CostDistance.scala    |  16 +--
 .../costdistance/MrGeoCostDistance.scala      | 135 +++++++++++++-----
 2 files changed, 103 insertions(+), 48 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 36f1f8cb36..130a3b6d9c 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -159,24 +159,18 @@ object CostDistance {
       * @param  q             The priority queue of candidate paths
       */
     def processNext(): Unit = {
-      val cost: Cost = q.poll
-      val (col, row, friction1, candidateCost) = cost
+      val entry: Cost = q.poll
+      val (col, row, friction1, candidateCost) = entry
       val currentCost =
         if (inTile(col, row))
           costTile.getDouble(col, row)
         else
           Double.NaN
 
-      // If the candidate path is a possible improvement ...
+      // If the candidate path is an improvement ...
       if (!isData(currentCost) || candidateCost <= currentCost) {
-
-        // Over-write the current cost with the candidate cost
-        if (inTile(col, row)) {
-          costTile.setDouble(col, row, candidateCost) // XXX
-
-          // Register changes on the boundary
-          if (onEdge(col, row)) edgeCallback(cost)
-        }
+        if (inTile(col, row)) costTile.setDouble(col, row, candidateCost) // XXX
+        if (onEdge(col, row)) edgeCallback(entry) // XXX
 
         // Compute candidate costs for neighbors and enqueue them
         if (isPassable(friction1)) {
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
index 52b08a03fd..747c437d4e 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
@@ -8,6 +8,7 @@ import geotrellis.raster.costdistance.CostDistance
 import org.apache.spark.rdd.RDD
 import org.apache.spark.SparkContext
 import org.apache.spark.util.AccumulatorV2
+import org.apache.spark.storage.StorageLevel
 
 import scala.collection.mutable
 
@@ -27,16 +28,16 @@ object MrGeoCostDistance {
 
   type CostList = List[CostDistance.Cost]
   type KeyListPair = (SpatialKey, CostList)
-  type PixelMap = mutable.Map[SpatialKey, CostList]
+  type Changes = mutable.Map[SpatialKey, CostList]
 
   /**
     * An accumulator to hold lists of edge changes.
     */
-  class PixelMapAccumulator extends AccumulatorV2[KeyListPair, PixelMap] {
-    private var map: PixelMap = mutable.Map.empty
+  class ChangesAccumulator extends AccumulatorV2[KeyListPair, Changes] {
+    private var map: Changes = mutable.Map.empty
 
-    def copy: PixelMapAccumulator = {
-      val other = new PixelMapAccumulator
+    def copy: ChangesAccumulator = {
+      val other = new ChangesAccumulator
       other.merge(this)
       other
     }
@@ -48,9 +49,9 @@ object MrGeoCostDistance {
       else map.put(key, map.get(key).get ++ list)
     }
     def isZero: Boolean = map.isEmpty
-    def merge(other: AccumulatorV2[(SpatialKey, CostList), PixelMap]): Unit = map ++= other.value
-    def reset: Unit = map = mutable.Map.empty
-    def value: PixelMap = map
+    def merge(other: AccumulatorV2[(SpatialKey, CostList), Changes]): Unit = map ++= other.value
+    def reset: Unit = map.clear
+    def value: Changes = map
   }
 
   /**
@@ -62,10 +63,13 @@ object MrGeoCostDistance {
     maxCost: Double = Double.PositiveInfinity
   )(implicit sc: SparkContext): RDD[(K, Tile)] = {
 
-    val accumulator = new PixelMapAccumulator
-    sc.register(accumulator, "Pixel Map Accumulator")
-
     val mt = friction.metadata.mapTransform
+    val bounds = friction.metadata.bounds.asInstanceOf[KeyBounds[K]]
+    val minKey = implicitly[SpatialKey](bounds.minKey)
+    val maxKey = implicitly[SpatialKey](bounds.maxKey)
+
+    val accumulator = new ChangesAccumulator
+    sc.register(accumulator, "Changes")
 
     // Load the accumulator with the starting values
     friction.foreach({ case (k, v) =>
@@ -83,7 +87,6 @@ object MrGeoCostDistance {
           val friction = tile.getDouble(col, row)
           (col, row, friction, 0.0)
         })
-
       accumulator.add((key, costs))
     })
 
@@ -92,34 +95,92 @@ object MrGeoCostDistance {
       val tile = implicitly[Tile](v)
       val cols = tile.cols
       val rows = tile.rows
-
       (k, v, CostDistance.generateEmptyCostTile(cols, rows))
     })
 
-    // XXX Should be shared variable
-    val changes = accumulator.value
-
-    // do {
-    costs = costs.map({ case (k, v, cost) =>
-      val key = implicitly[SpatialKey](k)
-      val friction = implicitly[Tile](v)
-      val cols = friction.cols
-      val rows = friction.rows
-      val q: CostDistance.Q = {
-        val q = CostDistance.generateEmptyQueue(cols, rows)
-
-        changes
-          .getOrElse(key, List.empty[CostDistance.Cost])
-          .foreach({ (entry: CostDistance.Cost) => q.add(entry) })
-        q
-      }
-
-      (k, v, CostDistance.compute(friction, cost, maxCost, q, CostDistance.nop))
-      // (k, v, CostDistance.compute(CostDistance.generateEmptyCostTile(cols, rows), cost, maxCost, q, CostDistance.nop))
-    })
-
-    accumulator.reset
-    // } while (changes.size > 0)
+    // Repeatedly map over the RDD of cost tiles until no more changes
+    // occur on the periphery of any tile.
+    do {
+      val changes = sc.broadcast(accumulator.value.toMap)
+      println(s"${changes.value.size}") // XXX
+
+      accumulator.reset
+
+      val previous = costs
+
+      costs = previous.map({ case (k, v, oldCostTile) =>
+        val key = implicitly[SpatialKey](k)
+        val frictionTile = implicitly[Tile](v)
+        val keyCol = key._1
+        val keyRow = key._2
+        val frictionTileCols = frictionTile.cols
+        val frictionTileRows = frictionTile.rows
+        val q: CostDistance.Q = {
+          val q = CostDistance.generateEmptyQueue(frictionTileCols, frictionTileRows)
+          changes.value
+            .getOrElse(key, List.empty[CostDistance.Cost])
+            .foreach({ (entry: CostDistance.Cost) => q.add(entry) })
+          q
+        }
+        val buffer: mutable.ArrayBuffer[CostDistance.Cost] = mutable.ArrayBuffer.empty
+        val newCostTile = CostDistance.compute(
+          frictionTile, oldCostTile,
+          maxCost, q,
+          { (entry: CostDistance.Cost) => buffer.append(entry) }
+        )
+
+        // Register changes on left periphery
+        if (minKey._1 <= keyCol-1) {
+          val leftKey = SpatialKey(keyCol-1, keyRow)
+          val leftList = buffer
+            .filter({ (entry: CostDistance.Cost) => entry._1 == 0 })
+            .map({ case (_, row: Int, f: Double, c: Double) =>
+              (frictionTileCols, row, f, c) })
+            .toList
+          if (leftList.size > 0)
+            accumulator.add((leftKey, leftList))
+        }
+
+        // Register changes on upper periphery
+        if (keyRow+1 <= maxKey._2) {
+          val upKey = SpatialKey(keyCol, keyRow+1)
+          val upList = buffer
+            .filter({ (entry: CostDistance.Cost) => entry._2 == frictionTileRows-1 })
+            .map({ case (col: Int, _, f: Double, c: Double) => (col, -1, f, c) })
+            .toList
+          if (upList.size > 0)
+            accumulator.add((upKey, upList))
+        }
+
+        // Register changes on right periphery
+        if (keyCol+1 <= maxKey._1) {
+          val rightKey = SpatialKey(keyCol+1, keyRow)
+          val rightList = buffer
+            .filter({ (entry: CostDistance.Cost) => entry._1 == frictionTileCols-1 })
+            .map({ case (_, row: Int, f: Double, c: Double) => (-1, row, f, c) })
+            .toList
+          if (rightList.size > 0)
+            accumulator.add((rightKey, rightList))
+        }
+
+        // Register changes on lower periphery
+        if (minKey._2 <= keyRow-1) {
+          val downKey = SpatialKey(keyCol, keyRow-1)
+          val downList = buffer
+            .filter({ (entry: CostDistance.Cost) => entry._2 == 0 })
+            .map({ case (col: Int, _, f: Double, c: Double) =>
+              (col, frictionTileRows, f, c) })
+            .toList
+          if (downList.size > 0)
+            accumulator.add((downKey, downList))
+        }
+
+        (k, v, newCostTile)
+      }).persist(StorageLevel.MEMORY_AND_DISK_SER)
+
+      costs.count
+      previous.unpersist()
+    } while (accumulator.value.size > 0)
 
     costs.map({ case (k, _, cost) => (k, cost) })
   }

From 225d8b27a7ff55ff23291081583cc8c95add4946 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Mon, 6 Feb 2017 11:41:15 -0500
Subject: [PATCH 14/28] Loop Fusion

---
 .../raster/costdistance/CostDistance.scala    |   9 +-
 .../costdistance/MrGeoCostDistance.scala      | 150 +++++++++---------
 2 files changed, 85 insertions(+), 74 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 130a3b6d9c..359aae3b9e 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -48,6 +48,10 @@ object CostDistance {
       })
   }
 
+  /**
+    * Generate an empty double-valued array tile of the correct
+    * dimensions.
+    */
   def generateEmptyCostTile(cols: Int, rows: Int): DoubleArrayTile =
     DoubleArrayTile.empty(cols, rows)
 
@@ -62,7 +66,6 @@ object CostDistance {
     *
     * @param  friction  Friction tile
     * @param  points    List of starting points as tuples
-    *
     */
   def apply(
     frictionTile: Tile,
@@ -169,8 +172,8 @@ object CostDistance {
 
       // If the candidate path is an improvement ...
       if (!isData(currentCost) || candidateCost <= currentCost) {
-        if (inTile(col, row)) costTile.setDouble(col, row, candidateCost) // XXX
-        if (onEdge(col, row)) edgeCallback(entry) // XXX
+        if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
+        if (onEdge(col, row)) edgeCallback(entry)
 
         // Compute candidate costs for neighbors and enqueue them
         if (isPassable(friction1)) {
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
index 747c437d4e..87bfcd212c 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
@@ -5,10 +5,11 @@ import geotrellis.spark._
 import geotrellis.vector._
 import geotrellis.raster.costdistance.CostDistance
 
+import org.apache.log4j.Logger
 import org.apache.spark.rdd.RDD
 import org.apache.spark.SparkContext
-import org.apache.spark.util.AccumulatorV2
 import org.apache.spark.storage.StorageLevel
+import org.apache.spark.util.AccumulatorV2
 
 import scala.collection.mutable
 
@@ -30,6 +31,8 @@ object MrGeoCostDistance {
   type KeyListPair = (SpatialKey, CostList)
   type Changes = mutable.Map[SpatialKey, CostList]
 
+  val logger = Logger.getLogger(MrGeoCostDistance.getClass)
+
   /**
     * An accumulator to hold lists of edge changes.
     */
@@ -71,8 +74,9 @@ object MrGeoCostDistance {
     val accumulator = new ChangesAccumulator
     sc.register(accumulator, "Changes")
 
-    // Load the accumulator with the starting values
-    friction.foreach({ case (k, v) =>
+    // Create RDD of initial (empty) cost tiles and load the
+    // accumulator with the starting values.
+    var costs: RDD[(K, V, DoubleArrayTile)] = friction.map({ case (k, v) =>
       val key = implicitly[SpatialKey](k)
       val tile = implicitly[Tile](v)
       val cols = tile.cols
@@ -87,14 +91,10 @@ object MrGeoCostDistance {
           val friction = tile.getDouble(col, row)
           (col, row, friction, 0.0)
         })
-      accumulator.add((key, costs))
-    })
 
-    // Create RDD of initial (empty) cost tiles
-    var costs: RDD[(K, V, DoubleArrayTile)] = friction.map({ case (k, v) =>
-      val tile = implicitly[Tile](v)
-      val cols = tile.cols
-      val rows = tile.rows
+      if (costs.length > 0)
+        accumulator.add((key, costs))
+
       (k, v, CostDistance.generateEmptyCostTile(cols, rows))
     })
 
@@ -102,7 +102,7 @@ object MrGeoCostDistance {
     // occur on the periphery of any tile.
     do {
       val changes = sc.broadcast(accumulator.value.toMap)
-      println(s"${changes.value.size}") // XXX
+      logger.debug(s"At least ${changes.value.size} changed tiles")
 
       accumulator.reset
 
@@ -115,67 +115,75 @@ object MrGeoCostDistance {
         val keyRow = key._2
         val frictionTileCols = frictionTile.cols
         val frictionTileRows = frictionTile.rows
-        val q: CostDistance.Q = {
-          val q = CostDistance.generateEmptyQueue(frictionTileCols, frictionTileRows)
-          changes.value
-            .getOrElse(key, List.empty[CostDistance.Cost])
-            .foreach({ (entry: CostDistance.Cost) => q.add(entry) })
-          q
-        }
-        val buffer: mutable.ArrayBuffer[CostDistance.Cost] = mutable.ArrayBuffer.empty
-        val newCostTile = CostDistance.compute(
-          frictionTile, oldCostTile,
-          maxCost, q,
-          { (entry: CostDistance.Cost) => buffer.append(entry) }
-        )
-
-        // Register changes on left periphery
-        if (minKey._1 <= keyCol-1) {
-          val leftKey = SpatialKey(keyCol-1, keyRow)
-          val leftList = buffer
-            .filter({ (entry: CostDistance.Cost) => entry._1 == 0 })
-            .map({ case (_, row: Int, f: Double, c: Double) =>
-              (frictionTileCols, row, f, c) })
-            .toList
-          if (leftList.size > 0)
-            accumulator.add((leftKey, leftList))
+        val localChanges = changes.value.getOrElse(key, List.empty[CostDistance.Cost])
+
+        if (localChanges.length > 0) {
+          val q: CostDistance.Q = CostDistance.generateEmptyQueue(frictionTileCols, frictionTileRows)
+          localChanges.foreach({ (entry: CostDistance.Cost) => q.add(entry) })
+
+          val buffer: mutable.ArrayBuffer[CostDistance.Cost] = mutable.ArrayBuffer.empty
+
+          val newCostTile = CostDistance.compute(
+            frictionTile, oldCostTile,
+            maxCost, q,
+            { (entry: CostDistance.Cost) => buffer.append(entry) }
+          )
+
+          // Register changes on left periphery
+          if (minKey._1 <= keyCol-1) {
+            val leftKey = SpatialKey(keyCol-1, keyRow)
+            val leftList = buffer
+              .filter({ (entry: CostDistance.Cost) => entry._1 == 0 })
+              .map({ case (_, row: Int, f: Double, c: Double) =>
+                (frictionTileCols, row, f, c) })
+              .toList
+            if (leftList.size > 0)
+              accumulator.add((leftKey, leftList))
+          }
+
+          // Register changes on upper periphery
+          if (keyRow+1 <= maxKey._2) {
+            val upKey = SpatialKey(keyCol, keyRow+1)
+            val upList = buffer
+              .filter({ (entry: CostDistance.Cost) => entry._2 == frictionTileRows-1 })
+              .map({ case (col: Int, _, f: Double, c: Double) => (col, -1, f, c) })
+              .toList
+            if (upList.size > 0)
+              accumulator.add((upKey, upList))
+          }
+
+          // Register changes on right periphery
+          if (keyCol+1 <= maxKey._1) {
+            val rightKey = SpatialKey(keyCol+1, keyRow)
+            val rightList = buffer
+              .filter({ (entry: CostDistance.Cost) => entry._1 == frictionTileCols-1 })
+              .map({ case (_, row: Int, f: Double, c: Double) => (-1, row, f, c) })
+              .toList
+            if (rightList.size > 0)
+              accumulator.add((rightKey, rightList))
+          }
+
+          // Register changes on lower periphery
+          if (minKey._2 <= keyRow-1) {
+            val downKey = SpatialKey(keyCol, keyRow-1)
+            val downList = buffer
+              .filter({ (entry: CostDistance.Cost) => entry._2 == 0 })
+              .map({ case (col: Int, _, f: Double, c: Double) =>
+                (col, frictionTileRows, f, c) })
+              .toList
+            if (downList.size > 0)
+              accumulator.add((downKey, downList))
+          }
+
+          // XXX It would be slightly more correct to include the four
+          // diagonal tiles as well, but there would be at most a one
+          // pixel contribution each, so it probably is not worth the
+          // expense.
+
+          (k, v, newCostTile)
         }
-
-        // Register changes on upper periphery
-        if (keyRow+1 <= maxKey._2) {
-          val upKey = SpatialKey(keyCol, keyRow+1)
-          val upList = buffer
-            .filter({ (entry: CostDistance.Cost) => entry._2 == frictionTileRows-1 })
-            .map({ case (col: Int, _, f: Double, c: Double) => (col, -1, f, c) })
-            .toList
-          if (upList.size > 0)
-            accumulator.add((upKey, upList))
-        }
-
-        // Register changes on right periphery
-        if (keyCol+1 <= maxKey._1) {
-          val rightKey = SpatialKey(keyCol+1, keyRow)
-          val rightList = buffer
-            .filter({ (entry: CostDistance.Cost) => entry._1 == frictionTileCols-1 })
-            .map({ case (_, row: Int, f: Double, c: Double) => (-1, row, f, c) })
-            .toList
-          if (rightList.size > 0)
-            accumulator.add((rightKey, rightList))
-        }
-
-        // Register changes on lower periphery
-        if (minKey._2 <= keyRow-1) {
-          val downKey = SpatialKey(keyCol, keyRow-1)
-          val downList = buffer
-            .filter({ (entry: CostDistance.Cost) => entry._2 == 0 })
-            .map({ case (col: Int, _, f: Double, c: Double) =>
-              (col, frictionTileRows, f, c) })
-            .toList
-          if (downList.size > 0)
-            accumulator.add((downKey, downList))
-        }
-
-        (k, v, newCostTile)
+        else
+          (k, v, oldCostTile)
       }).persist(StorageLevel.MEMORY_AND_DISK_SER)
 
       costs.count

From 610ab937efd4ae3945e2f267e4929d6839dde33b Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Mon, 6 Feb 2017 14:08:59 -0500
Subject: [PATCH 15/28] Rename Object

---
 .../{MrGeoCostDistance.scala => IterativeCostDistance.scala}  | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)
 rename spark/src/main/scala/geotrellis/spark/costdistance/{MrGeoCostDistance.scala => IterativeCostDistance.scala} (98%)

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
similarity index 98%
rename from spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
rename to spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 87bfcd212c..d500448eb3 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/MrGeoCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -25,13 +25,13 @@ import scala.collection.mutable
   *
   * 2. https://github.com/ngageoint/mrgeo/blob/0c6ed4a7e66bb0923ec5c570b102862aee9e885e/mrgeo-mapalgebra/mrgeo-mapalgebra-costdistance/src/main/scala/org/mrgeo/mapalgebra/CostDistanceMapOp.scala
   */
-object MrGeoCostDistance {
+object IterativeCostDistance {
 
   type CostList = List[CostDistance.Cost]
   type KeyListPair = (SpatialKey, CostList)
   type Changes = mutable.Map[SpatialKey, CostList]
 
-  val logger = Logger.getLogger(MrGeoCostDistance.getClass)
+  val logger = Logger.getLogger(IterativeCostDistance.getClass)
 
   /**
     * An accumulator to hold lists of edge changes.

From 5c4c9e88d70f008ea73ccb971875e02009fd2fc9 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Mon, 6 Feb 2017 16:52:26 -0500
Subject: [PATCH 16/28] Compute Layer Resolution

---
 .../costdistance/CostDistanceSpec.scala       |  4 +--
 .../raster/costdistance/CostDistance.scala    | 31 ++++++++++++-------
 .../costdistance/IterativeCostDistance.scala  | 28 ++++++++++++++---
 3 files changed, 46 insertions(+), 17 deletions(-)

diff --git a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
index 48c3e31a43..01fa082f83 100644
--- a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
+++ b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
@@ -157,8 +157,8 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
 
     CostDistance.compute(
       frictionTile, costTile,
-      Double.PositiveInfinity, q,
-      edgeCb)
+      Double.PositiveInfinity, 1,
+      q, edgeCb)
 
     actualLeft should be (expectedLeft)
     actualUp should be (expectedUp)
diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 359aae3b9e..87440d4bdd 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -39,6 +39,9 @@ object CostDistance {
   /**
     * Generate a Queue suitable for working with a tile of the given
     * dimensions.
+    *
+    * @param  cols  The number of columns of the friction tile
+    * @param  rows  The number of rows of the friction tile
     */
   def generateEmptyQueue(cols: Int, rows: Int): Q = {
     new PriorityQueue(
@@ -51,6 +54,9 @@ object CostDistance {
   /**
     * Generate an empty double-valued array tile of the correct
     * dimensions.
+    *
+    * @param  cols  The number of cols of the friction tile (and therefore the cost tile)
+    * @param  rows  The number of rows of the frition tile and cost tiles
     */
   def generateEmptyCostTile(cols: Int, rows: Int): DoubleArrayTile =
     DoubleArrayTile.empty(cols, rows)
@@ -64,13 +70,16 @@ object CostDistance {
     *    "Propagating radial waves of travel cost in a grid."
     *    International Journal of Geographical Information Science 24.9 (2010): 1391-1413.
     *
-    * @param  friction  Friction tile
-    * @param  points    List of starting points as tuples
+    * @param  friction    Friction tile; pixels are interpreted as "second per meter"
+    * @param  points      List of starting points as tuples
+    * @param  maxCost     The maximum cost before pruning a path (in units of "seconds")
+    * @param  resolution  The resolution of the tiles (in units of "meters per pixel")
     */
   def apply(
     frictionTile: Tile,
     points: Seq[(Int, Int)],
-    maxCost: Double = Double.PositiveInfinity
+    maxCost: Double = Double.PositiveInfinity,
+    resolution: Double = 1
   ): DoubleArrayTile = {
     val cols = frictionTile.cols
     val rows = frictionTile.rows
@@ -82,7 +91,7 @@ object CostDistance {
       q.add(entry)
     })
 
-    compute(frictionTile, costTile, maxCost, q, nop)
+    compute(frictionTile, costTile, maxCost, resolution, q, nop)
   }
 
   /**
@@ -90,16 +99,16 @@ object CostDistance {
     *
     * @param  frictionTile    The friction tile
     * @param  costTile        The tile that will contain the costs
-    * @param  maxCost         The maximum cost of any path (truncates to limit computational cost)
+    * @param  maxCost         The maximum cost before pruning a path (in units of "seconds")
+    * @param  resolution      The resolution of the tiles (in units of "meters per pixel")
     * @param  q               A priority queue of Cost objects (a.k.a. candidate paths)
     * @param  edgeCallback    Called when a pixel on the edge of the tile is updated
     */
   def compute(
     frictionTile: Tile,
     costTile: DoubleArrayTile,
-    maxCost: Double,
-    q: Q,
-    edgeCallback: EdgeCallback
+    maxCost: Double, resolution: Double,
+    q: Q, edgeCallback: EdgeCallback
   ): DoubleArrayTile = {
     val cols = frictionTile.cols
     val rows = frictionTile.rows
@@ -128,8 +137,7 @@ object CostDistance {
       * @param  distance      The distance from the neighboring location to this location
       */
     @inline def enqueueNeighbor(
-      col: Int, row: Int, friction1: Double,
-      cost: Double,
+      col: Int, row: Int, friction1: Double, cost: Double,
       distance: Double = 1.0
     ): Unit = {
       // If the location is inside of the tile ...
@@ -139,7 +147,8 @@ object CostDistance {
 
         // ... and if the location is passable ...
         if (isPassable(friction2)) {
-          val entry = (col, row, friction2, cost + distance * (friction1 + friction2) / 2.0)
+          val step = resolution * distance * (friction1 + friction2) / 2.0
+          val entry = (col, row, friction2, cost + step)
           val candidateCost = entry._4
 
           // ... and the candidate cost is less than the maximum cost ...
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index d500448eb3..9db5be5936 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -1,9 +1,10 @@
 package geotrellis.spark.costdistance
 
+import geotrellis.proj4.LatLng
 import geotrellis.raster._
+import geotrellis.raster.costdistance.CostDistance
 import geotrellis.spark._
 import geotrellis.vector._
-import geotrellis.raster.costdistance.CostDistance
 
 import org.apache.log4j.Logger
 import org.apache.spark.rdd.RDD
@@ -59,6 +60,10 @@ object IterativeCostDistance {
 
   /**
     * Perform the cost-distance computation.
+    *
+    * @param  friction  The friction layer; pixels are in units of "seconds per meter"
+    * @param  points    The starting locations from-which to compute the cost of traveling
+    * @param  maxCost   The maximum cost before pruning a path (in units of "seconds")
     */
   def apply[K: (? => SpatialKey), V: (? => Tile)](
     friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]],
@@ -66,7 +71,22 @@ object IterativeCostDistance {
     maxCost: Double = Double.PositiveInfinity
   )(implicit sc: SparkContext): RDD[(K, Tile)] = {
 
-    val mt = friction.metadata.mapTransform
+    val md = friction.metadata
+
+    val mt = md.mapTransform
+
+    val resolution = {
+      val kv = friction.first
+      val key = implicitly[SpatialKey](kv._1)
+      val tile = implicitly[Tile](kv._2)
+      val extent = mt(key).reproject(md.crs, LatLng)
+      val degrees = extent.xmax - extent.xmin
+      val meters = degrees * (6378137 * 2.0 * math.Pi) / 360.0
+      val pixels = tile.cols
+      math.abs(meters / pixels)
+    }
+    logger.debug(s"Computed resolution: $resolution meters/pixel")
+
     val bounds = friction.metadata.bounds.asInstanceOf[KeyBounds[K]]
     val minKey = implicitly[SpatialKey](bounds.minKey)
     val maxKey = implicitly[SpatialKey](bounds.maxKey)
@@ -125,8 +145,8 @@ object IterativeCostDistance {
 
           val newCostTile = CostDistance.compute(
             frictionTile, oldCostTile,
-            maxCost, q,
-            { (entry: CostDistance.Cost) => buffer.append(entry) }
+            maxCost, resolution,
+            q, { (entry: CostDistance.Cost) => buffer.append(entry) }
           )
 
           // Register changes on left periphery

From 5c716f3451626eae85c65e7285c26aa51fa16ebb Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Tue, 7 Feb 2017 15:14:44 -0500
Subject: [PATCH 17/28] Do Not Name Accumulator

Naming this accumulator adds a lot of clutter to the Spark UI.
---
 .../geotrellis/spark/costdistance/IterativeCostDistance.scala   | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 9db5be5936..fcc52c6d97 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -92,7 +92,7 @@ object IterativeCostDistance {
     val maxKey = implicitly[SpatialKey](bounds.maxKey)
 
     val accumulator = new ChangesAccumulator
-    sc.register(accumulator, "Changes")
+    sc.register(accumulator)
 
     // Create RDD of initial (empty) cost tiles and load the
     // accumulator with the starting values.

From 0e4de08d5a834e9d704dce0cbe4d4f3bba4c30aa Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 9 Feb 2017 09:08:17 -0500
Subject: [PATCH 18/28] Iterative Cost-Distance Unit Tests

---
 .../costdistance/IterativeCostDistance.scala  |  28 ++--
 .../costdistance/IterativeCostDistance.scala  | 120 ++++++++++++++++++
 .../RDDHistogramEqualizationSpec.scala        |   2 +-
 3 files changed, 137 insertions(+), 13 deletions(-)
 create mode 100644 spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index fcc52c6d97..7516de7489 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -58,6 +58,21 @@ object IterativeCostDistance {
     def value: Changes = map
   }
 
+  def computeResolution[K: (? => SpatialKey), V: (? => Tile)](
+    friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]]
+  ) = {
+    val md = friction.metadata
+    val mt = md.mapTransform
+    val kv = friction.first
+    val key = implicitly[SpatialKey](kv._1)
+    val tile = implicitly[Tile](kv._2)
+    val extent = mt(key).reproject(md.crs, LatLng)
+    val degrees = extent.xmax - extent.xmin
+    val meters = degrees * (6378137 * 2.0 * math.Pi) / 360.0
+    val pixels = tile.cols
+    math.abs(meters / pixels)
+  }
+
   /**
     * Perform the cost-distance computation.
     *
@@ -72,19 +87,8 @@ object IterativeCostDistance {
   )(implicit sc: SparkContext): RDD[(K, Tile)] = {
 
     val md = friction.metadata
-
     val mt = md.mapTransform
-
-    val resolution = {
-      val kv = friction.first
-      val key = implicitly[SpatialKey](kv._1)
-      val tile = implicitly[Tile](kv._2)
-      val extent = mt(key).reproject(md.crs, LatLng)
-      val degrees = extent.xmax - extent.xmin
-      val meters = degrees * (6378137 * 2.0 * math.Pi) / 360.0
-      val pixels = tile.cols
-      math.abs(meters / pixels)
-    }
+    val resolution = computeResolution(friction)
     logger.debug(s"Computed resolution: $resolution meters/pixel")
 
     val bounds = friction.metadata.bounds.asInstanceOf[KeyBounds[K]]
diff --git a/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
new file mode 100644
index 0000000000..523ac44617
--- /dev/null
+++ b/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -0,0 +1,120 @@
+/*
+ * Copyright 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.spark.costdistance
+
+import geotrellis.proj4.LatLng
+import geotrellis.raster._
+import geotrellis.spark._
+import geotrellis.spark.tiling.LayoutDefinition
+import geotrellis.vector._
+
+import org.scalatest._
+
+
+class IterativeCostDistanceSpec extends FunSpec
+    with Matchers
+    with TestEnvironment {
+
+  val rdd1 = {
+    val tile = IntArrayTile(Array.fill[Int](25)(1), 5, 5)
+    val skey = SpatialKey(0, 0)
+    val extent = Extent(0, 0, 5, 5)
+    val gridExtent = GridExtent(extent, 1, 1) // 5×5 pixels
+    val layoutDefinition = LayoutDefinition(gridExtent, 5)
+    val bounds = Bounds(skey, skey)
+    val tileLayerMetadata = TileLayerMetadata(IntCellType, layoutDefinition, extent, LatLng, bounds)
+    val list = List((skey, tile))
+    ContextRDD(sc.parallelize(list), tileLayerMetadata)
+  }
+
+  val rdd2 = {
+    val tile = IntArrayTile(Array.fill[Int](25)(1), 5, 5)
+    val extent = Extent(0, 0, 10, 5)
+    val gridExtent = GridExtent(extent, 1, 1) // 10×5 pixels
+    val layoutDefinition = LayoutDefinition(gridExtent, 10, 5)
+    val bounds = Bounds(SpatialKey(0,0), SpatialKey(1,0))
+    val tileLayerMetadata = TileLayerMetadata(IntCellType, layoutDefinition, extent, LatLng, bounds)
+    val list = List((SpatialKey(0,0), tile), (SpatialKey(1,0), tile))
+    ContextRDD(sc.parallelize(list), tileLayerMetadata)
+  }
+
+  val rdd3 = {
+    val tile = IntArrayTile(Array.fill[Int](25)(1), 5, 5)
+    val extent = Extent(0, 0, 5, 10)
+    val gridExtent = GridExtent(extent, 1, 1) // 5×10 pixels
+    val layoutDefinition = LayoutDefinition(gridExtent, 5, 10)
+    val bounds = Bounds(SpatialKey(0,0), SpatialKey(0,1))
+    val tileLayerMetadata = TileLayerMetadata(IntCellType, layoutDefinition, extent, LatLng, bounds)
+    val list = List((SpatialKey(0,0), tile), (SpatialKey(0,1), tile))
+    ContextRDD(sc.parallelize(list), tileLayerMetadata)
+  }
+
+  describe("Iterative Cost Distance") {
+
+    it("Should correctly compute resolution") {
+      val resolution = IterativeCostDistance.computeResolution(rdd1)
+      math.floor(resolution/1000) should be (111) // ≈ 111 kilometers per degree
+    }
+
+    it("Should correctly project input points") {
+      val costs = IterativeCostDistance(rdd1, List(Point(2.5, 2.5)))
+      val cost = costs.first._2
+      cost.getDouble(2,2) should be (0.0)
+    }
+
+    it("Should propogate left") {
+      val costs = IterativeCostDistance(rdd2, List(Point(2.5+5.0, 2.5)))
+      val right = costs.filter({ case (k, _) => k == SpatialKey(1, 0) }).first._2
+      val left = costs.filter({ case (k, _) => k == SpatialKey(0, 0) }).first._2
+      val resolution = IterativeCostDistance.computeResolution(rdd2)
+      val hops = (right.getDouble(3,2) - left.getDouble(3,2)) / resolution
+
+      hops should be (5.0)
+    }
+
+    it("Should propogate right") {
+      val costs = IterativeCostDistance(rdd2, List(Point(2.5, 2.5)))
+      val right = costs.filter({ case (k, _) => k == SpatialKey(1, 0) }).first._2
+      val left = costs.filter({ case (k, _) => k == SpatialKey(0, 0) }).first._2
+      val resolution = IterativeCostDistance.computeResolution(rdd2)
+      val hops = (right.getDouble(1,2) - left.getDouble(1,2)) / resolution
+
+      hops should be (5.0)
+    }
+
+    it("Should propogate up") {
+      val costs = IterativeCostDistance(rdd3, List(Point(2.5, 2.5)))
+      val up = costs.filter({ case (k, _) => k == SpatialKey(0, 1) }).first._2
+      val down = costs.filter({ case (k, _) => k == SpatialKey(0, 0) }).first._2
+      val resolution = IterativeCostDistance.computeResolution(rdd3)
+      val hops = (up.getDouble(2,3) - down.getDouble(2,3)) / resolution
+
+      hops should be (5.0)
+    }
+
+    it("Should propogate down") {
+      val costs = IterativeCostDistance(rdd3, List(Point(2.5, 2.5+5.0)))
+      val up = costs.filter({ case (k, _) => k == SpatialKey(0, 1) }).first._2
+      val down = costs.filter({ case (k, _) => k == SpatialKey(0, 0) }).first._2
+      val resolution = IterativeCostDistance.computeResolution(rdd3)
+      val hops = (up.getDouble(2,1) - down.getDouble(2,1)) / resolution
+
+      hops should be (5.0)
+    }
+
+  }
+}
diff --git a/spark/src/test/scala/geotrellis/spark/equalization/RDDHistogramEqualizationSpec.scala b/spark/src/test/scala/geotrellis/spark/equalization/RDDHistogramEqualizationSpec.scala
index d101caa0d5..4fa4f926aa 100644
--- a/spark/src/test/scala/geotrellis/spark/equalization/RDDHistogramEqualizationSpec.scala
+++ b/spark/src/test/scala/geotrellis/spark/equalization/RDDHistogramEqualizationSpec.scala
@@ -40,7 +40,7 @@ class RDDHistogramEqualizationSpec extends FunSpec
 
   describe("RDD Histogram Equalization") {
 
-    it("should work with float-point rasters") {
+    it("should work with floating-point rasters") {
       val tile1: Tile = DoubleArrayTile(data1.map(_.toDouble).toArray, 1, 8)
       val tile2: Tile = DoubleArrayTile(data2.map(_.toDouble).toArray, 1, 8)
       val rdd = ContextRDD(sc.parallelize(List((0, tile1), (1, tile2))), 33).equalize

From 3ec6ff1b821a8b86d11d4eb6586ccc69f8371cc6 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 9 Feb 2017 14:13:53 -0500
Subject: [PATCH 19/28] Address Not-Fun, Not-Funny Concurrency Bug

---
 .../raster/costdistance/CostDistance.scala    |  10 +-
 .../costdistance/IterativeCostDistance.scala  | 149 +++++++-----------
 2 files changed, 65 insertions(+), 94 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 87440d4bdd..28f99d909e 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -119,7 +119,7 @@ object CostDistance {
       ((0 <= col && col < cols) && (0 <= row && row < rows))
 
     def isPassable(f: Double): Boolean =
-      (isData(f) && 0 <= f)
+      (isData(f) && 0.0 <= f)
 
     def onEdge(col: Int, row: Int): Boolean =
       ((col == 0) || (row == 0) || (col == cols-1) || (row == rows-1))
@@ -137,7 +137,7 @@ object CostDistance {
       * @param  distance      The distance from the neighboring location to this location
       */
     @inline def enqueueNeighbor(
-      col: Int, row: Int, friction1: Double, cost: Double,
+      col: Int, row: Int, friction1: Double, neighborCost: Double,
       distance: Double = 1.0
     ): Unit = {
       // If the location is inside of the tile ...
@@ -148,13 +148,13 @@ object CostDistance {
         // ... and if the location is passable ...
         if (isPassable(friction2)) {
           val step = resolution * distance * (friction1 + friction2) / 2.0
-          val entry = (col, row, friction2, cost + step)
-          val candidateCost = entry._4
+          val candidateCost = neighborCost + step
+          val entry = (col, row, friction2, candidateCost)
 
           // ... and the candidate cost is less than the maximum cost ...
           if (candidateCost <= maxCost) {
             // ... and the candidate is a possible improvement ...
-            if ((isData(currentCost) && candidateCost < currentCost) || !isData(currentCost)) {
+            if (!isData(currentCost) || candidateCost < currentCost) {
               costTile.setDouble(col, row, candidateCost) // then increase lower bound on pixel,
               q.add(entry) // and enqueue candidate for future processing
             }
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 7516de7489..059a3a8ea0 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -28,34 +28,30 @@ import scala.collection.mutable
   */
 object IterativeCostDistance {
 
-  type CostList = List[CostDistance.Cost]
-  type KeyListPair = (SpatialKey, CostList)
-  type Changes = mutable.Map[SpatialKey, CostList]
+  type KeyCostPair = (SpatialKey, CostDistance.Cost)
+  type Changes = mutable.ArrayBuffer[KeyCostPair]
 
   val logger = Logger.getLogger(IterativeCostDistance.getClass)
 
   /**
     * An accumulator to hold lists of edge changes.
     */
-  class ChangesAccumulator extends AccumulatorV2[KeyListPair, Changes] {
-    private var map: Changes = mutable.Map.empty
+  class ChangesAccumulator extends AccumulatorV2[KeyCostPair, Changes] {
+    private val list: Changes = mutable.ArrayBuffer.empty
 
     def copy: ChangesAccumulator = {
       val other = new ChangesAccumulator
       other.merge(this)
       other
     }
-    def add(kv: (SpatialKey, CostList)): Unit = {
-      val key = kv._1
-      val list = kv._2
-
-      if (!map.contains(key)) map.put(key, list)
-      else map.put(key, map.get(key).get ++ list)
+    def add(pair: KeyCostPair): Unit = {
+      this.synchronized { list.append(pair) }
     }
-    def isZero: Boolean = map.isEmpty
-    def merge(other: AccumulatorV2[(SpatialKey, CostList), Changes]): Unit = map ++= other.value
-    def reset: Unit = map.clear
-    def value: Changes = map
+    def isZero: Boolean = list.isEmpty
+    def merge(other: AccumulatorV2[KeyCostPair, Changes]): Unit =
+      this.synchronized { list ++= other.value }
+    def reset: Unit = this.synchronized { list.clear }
+    def value: Changes = list
   }
 
   def computeResolution[K: (? => SpatialKey), V: (? => Tile)](
@@ -93,7 +89,11 @@ object IterativeCostDistance {
 
     val bounds = friction.metadata.bounds.asInstanceOf[KeyBounds[K]]
     val minKey = implicitly[SpatialKey](bounds.minKey)
+    val minKeyCol = minKey._1
+    val minKeyRow = minKey._2
     val maxKey = implicitly[SpatialKey](bounds.maxKey)
+    val maxKeyCol = maxKey._1
+    val maxKeyRow = maxKey._2
 
     val accumulator = new ChangesAccumulator
     sc.register(accumulator)
@@ -107,25 +107,31 @@ object IterativeCostDistance {
       val rows = tile.rows
       val extent = mt(key)
       val rasterExtent = RasterExtent(extent, cols, rows)
-      val costs: List[CostDistance.Cost] = points
+
+      points
         .filter({ point => extent.contains(point) })
         .map({ point =>
           val col = rasterExtent.mapXToGrid(point.x)
           val row = rasterExtent.mapYToGrid(point.y)
           val friction = tile.getDouble(col, row)
-          (col, row, friction, 0.0)
+          val cost = (col, row, friction, 0.0)
+          accumulator.add((key, cost))
         })
 
-      if (costs.length > 0)
-        accumulator.add((key, costs))
-
       (k, v, CostDistance.generateEmptyCostTile(cols, rows))
-    })
+    }).persist(StorageLevel.MEMORY_AND_DISK_SER)
+
+    costs.count
 
     // Repeatedly map over the RDD of cost tiles until no more changes
     // occur on the periphery of any tile.
     do {
-      val changes = sc.broadcast(accumulator.value.toMap)
+      val _changes: Map[SpatialKey, Seq[CostDistance.Cost]] =
+        accumulator.value
+          .groupBy(_._1)
+          .map({ case (k, list) => (k, list.map({ case (_, v) => v })) })
+          .toMap
+      val changes = sc.broadcast(_changes)
       logger.debug(s"At least ${changes.value.size} changed tiles")
 
       accumulator.reset
@@ -139,75 +145,40 @@ object IterativeCostDistance {
         val keyRow = key._2
         val frictionTileCols = frictionTile.cols
         val frictionTileRows = frictionTile.rows
-        val localChanges = changes.value.getOrElse(key, List.empty[CostDistance.Cost])
-
-        if (localChanges.length > 0) {
-          val q: CostDistance.Q = CostDistance.generateEmptyQueue(frictionTileCols, frictionTileRows)
-          localChanges.foreach({ (entry: CostDistance.Cost) => q.add(entry) })
-
-          val buffer: mutable.ArrayBuffer[CostDistance.Cost] = mutable.ArrayBuffer.empty
-
-          val newCostTile = CostDistance.compute(
-            frictionTile, oldCostTile,
-            maxCost, resolution,
-            q, { (entry: CostDistance.Cost) => buffer.append(entry) }
-          )
-
-          // Register changes on left periphery
-          if (minKey._1 <= keyCol-1) {
-            val leftKey = SpatialKey(keyCol-1, keyRow)
-            val leftList = buffer
-              .filter({ (entry: CostDistance.Cost) => entry._1 == 0 })
-              .map({ case (_, row: Int, f: Double, c: Double) =>
-                (frictionTileCols, row, f, c) })
-              .toList
-            if (leftList.size > 0)
-              accumulator.add((leftKey, leftList))
-          }
-
-          // Register changes on upper periphery
-          if (keyRow+1 <= maxKey._2) {
-            val upKey = SpatialKey(keyCol, keyRow+1)
-            val upList = buffer
-              .filter({ (entry: CostDistance.Cost) => entry._2 == frictionTileRows-1 })
-              .map({ case (col: Int, _, f: Double, c: Double) => (col, -1, f, c) })
-              .toList
-            if (upList.size > 0)
-              accumulator.add((upKey, upList))
+        val localChanges: Option[Seq[CostDistance.Cost]] = changes.value.get(key)
+
+        localChanges match {
+          case Some(localChanges) => {
+            val q: CostDistance.Q = {
+              val q = CostDistance.generateEmptyQueue(frictionTileCols, frictionTileRows)
+              localChanges.foreach({ (entry: CostDistance.Cost) => q.add(entry) })
+              q
+            }
+
+            val newCostTile = CostDistance.compute(
+              frictionTile, oldCostTile,
+              maxCost, resolution,
+              q, { (entry: CostDistance.Cost) =>
+                val (col, row, f, c) = entry
+                if (col == 0 && (minKeyCol <= keyCol-1))
+                  accumulator.add((SpatialKey(keyCol-1, keyRow), (frictionTileCols, row, f, c)))
+                if (row == frictionTileRows-1 && (keyRow+1 <= maxKeyRow))
+                  accumulator.add((SpatialKey(keyCol, keyRow+1), (col, -1, f, c)))
+                if (col == frictionTileCols-1 && (keyCol+1 <= maxKeyCol))
+                  accumulator.add((SpatialKey(keyCol+1, keyRow), (-1, row, f, c)))
+                if (row == 0 && (minKeyRow <= keyRow-1))
+                  accumulator.add((SpatialKey(keyCol, keyRow-1)), (col, frictionTileRows, f, c))
+              })
+
+            // XXX It would be slightly more correct to include the four
+            // diagonal tiles as well, but there would be at most a one
+            // pixel contribution each, so it probably is not worth the
+            // expense.
+
+            (k, v, newCostTile)
           }
-
-          // Register changes on right periphery
-          if (keyCol+1 <= maxKey._1) {
-            val rightKey = SpatialKey(keyCol+1, keyRow)
-            val rightList = buffer
-              .filter({ (entry: CostDistance.Cost) => entry._1 == frictionTileCols-1 })
-              .map({ case (_, row: Int, f: Double, c: Double) => (-1, row, f, c) })
-              .toList
-            if (rightList.size > 0)
-              accumulator.add((rightKey, rightList))
-          }
-
-          // Register changes on lower periphery
-          if (minKey._2 <= keyRow-1) {
-            val downKey = SpatialKey(keyCol, keyRow-1)
-            val downList = buffer
-              .filter({ (entry: CostDistance.Cost) => entry._2 == 0 })
-              .map({ case (col: Int, _, f: Double, c: Double) =>
-                (col, frictionTileRows, f, c) })
-              .toList
-            if (downList.size > 0)
-              accumulator.add((downKey, downList))
-          }
-
-          // XXX It would be slightly more correct to include the four
-          // diagonal tiles as well, but there would be at most a one
-          // pixel contribution each, so it probably is not worth the
-          // expense.
-
-          (k, v, newCostTile)
+          case None => (k, v, oldCostTile)
         }
-        else
-          (k, v, oldCostTile)
       }).persist(StorageLevel.MEMORY_AND_DISK_SER)
 
       costs.count

From 48a8b362a31ff9c9fdcda81e764cec1ef9c504d0 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Thu, 9 Feb 2017 18:10:57 -0500
Subject: [PATCH 20/28] Add Diagonals

---
 .../costdistance/IterativeCostDistance.scala  | 30 ++++++++++++-------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 059a3a8ea0..9e2378276a 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -160,20 +160,30 @@ object IterativeCostDistance {
               maxCost, resolution,
               q, { (entry: CostDistance.Cost) =>
                 val (col, row, f, c) = entry
-                if (col == 0 && (minKeyCol <= keyCol-1))
+                if (col == 0 && (minKeyCol <= keyCol-1)) // left
                   accumulator.add((SpatialKey(keyCol-1, keyRow), (frictionTileCols, row, f, c)))
-                if (row == frictionTileRows-1 && (keyRow+1 <= maxKeyRow))
+
+                if (row == frictionTileRows-1 && (keyRow+1 <= maxKeyRow)) // up
                   accumulator.add((SpatialKey(keyCol, keyRow+1), (col, -1, f, c)))
-                if (col == frictionTileCols-1 && (keyCol+1 <= maxKeyCol))
+
+                if (col == frictionTileCols-1 && (keyCol+1 <= maxKeyCol)) // right
                   accumulator.add((SpatialKey(keyCol+1, keyRow), (-1, row, f, c)))
-                if (row == 0 && (minKeyRow <= keyRow-1))
-                  accumulator.add((SpatialKey(keyCol, keyRow-1)), (col, frictionTileRows, f, c))
-              })
 
-            // XXX It would be slightly more correct to include the four
-            // diagonal tiles as well, but there would be at most a one
-            // pixel contribution each, so it probably is not worth the
-            // expense.
+                if (row == 0 && (minKeyRow <= keyRow-1)) // down
+                  accumulator.add((SpatialKey(keyCol, keyRow-1), (col, frictionTileRows, f, c)))
+
+                if (col == 0 && row == 0 && (minKeyCol <= keyCol-1) && (minKeyRow <= keyRow-1)) // upper-left
+                  accumulator.add((SpatialKey(keyCol-1,keyRow-1), (frictionTileCols, frictionTileRows, f, c)))
+
+                if (col == frictionTileCols-1 && row == 0 && (keyCol+1 <= maxKeyCol) && (minKeyRow <= keyRow-1)) // upper-right
+                  accumulator.add((SpatialKey(keyCol+1,keyRow-1), (-1, frictionTileRows, f, c)))
+
+                if (col == frictionTileCols-1 && row == frictionTileRows-1 && (keyCol+1 <= maxKeyCol) && (keyRow+1 <= maxKeyCol)) // lower-right
+                  accumulator.add((SpatialKey(keyCol+1,keyRow+1), (-1, -1, f, c)))
+
+                if (col == 0 && row == frictionTileRows-1 && (minKeyCol <= keyCol-1) && (keyRow+1 <= maxKeyRow)) // lower-left
+                  accumulator.add((SpatialKey(keyCol-1,keyRow+1), (frictionTileCols, -1, f, c)))
+              })
 
             (k, v, newCostTile)
           }

From 44984238807b977512b184b63b97beed625a453e Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Fri, 10 Feb 2017 12:24:09 -0500
Subject: [PATCH 21/28] Update Comments

---
 .../scala/geotrellis/raster/costdistance/CostDistance.scala  | 2 +-
 .../spark/costdistance/IterativeCostDistance.scala           | 5 +++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 28f99d909e..69e1398174 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -64,7 +64,7 @@ object CostDistance {
   /**
     * Generate a cost-distance raster based on a set of starting
     * points and a friction raster.  This is an implementation of the
-    * standard algorithm from [1].
+    * standard algorithm mentioned in the "previous work" section of [1].
     *
     * 1. Tomlin, Dana.
     *    "Propagating radial waves of travel cost in a grid."
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 9e2378276a..07320a6615 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -17,8 +17,9 @@ import scala.collection.mutable
 
 /**
   * This Spark-enabled implementation of the standard cost-distance
-  * algorithm [1] is "heavily inspired" by the MrGeo implementation
-  * [2] but does not share any code with it.
+  * algorithm mentioned in the "previous work" section of [1] is
+  * "heavily inspired" by the MrGeo implementation [2] but does not
+  * share any code with it.
   *
   * 1. Tomlin, Dana.
   *    "Propagating radial waves of travel cost in a grid."

From ec8bc2a78244eb12a262ae9a2ffaf904610faca2 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Sat, 18 Feb 2017 19:38:37 -0500
Subject: [PATCH 22/28] Add Extension Methods

---
 .../spark/costdistance/Implicits.scala        | 30 +++++++++++++
 .../costdistance/IterativeCostDistance.scala  | 26 ++++++++++--
 .../costdistance/RDDCostDistanceMethods.scala | 42 +++++++++++++++++++
 .../main/scala/geotrellis/spark/package.scala |  1 +
 .../costdistance/IterativeCostDistance.scala  |  2 +-
 5 files changed, 97 insertions(+), 4 deletions(-)
 create mode 100644 spark/src/main/scala/geotrellis/spark/costdistance/Implicits.scala
 create mode 100644 spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/Implicits.scala b/spark/src/main/scala/geotrellis/spark/costdistance/Implicits.scala
new file mode 100644
index 0000000000..4b78e16917
--- /dev/null
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/Implicits.scala
@@ -0,0 +1,30 @@
+/*
+ * Copyright 2017 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.spark.costdistance
+
+import geotrellis.spark._
+import geotrellis.raster.Tile
+
+import org.apache.spark.rdd.RDD
+
+
+object Implicits extends Implicits
+
+trait Implicits {
+  implicit class withRDDCostDistanceMethods[K: (? => SpatialKey), V: (? => Tile)](val self: RDD[(K, V)] with Metadata[TileLayerMetadata[K]])
+      extends RDDCostDistanceMethods[K, V]
+}
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 07320a6615..5fa2080c5f 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -1,3 +1,19 @@
+/*
+ * Copyright 2017 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.spark.costdistance
 
 import geotrellis.proj4.LatLng
@@ -79,9 +95,9 @@ object IterativeCostDistance {
     */
   def apply[K: (? => SpatialKey), V: (? => Tile)](
     friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]],
-    points: List[Point],
+    points: Seq[Point],
     maxCost: Double = Double.PositiveInfinity
-  )(implicit sc: SparkContext): RDD[(K, Tile)] = {
+  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] with Metadata[TileLayerMetadata[K]]= {
 
     val md = friction.metadata
     val mt = md.mapTransform
@@ -196,6 +212,10 @@ object IterativeCostDistance {
       previous.unpersist()
     } while (accumulator.value.size > 0)
 
-    costs.map({ case (k, _, cost) => (k, cost) })
+    // Construct return value and return it
+    val metadata = TileLayerMetadata(DoubleCellType, md.layout, md.extent, md.crs, md.bounds)
+    val rdd = costs.map({ case (k, _, cost) => (k, cost) })
+    ContextRDD(rdd, metadata)
   }
+
 }
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala b/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
new file mode 100644
index 0000000000..ad3f5f2001
--- /dev/null
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
@@ -0,0 +1,42 @@
+/*
+ * Copyright 2017 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.spark.costdistance
+
+import geotrellis.raster.{Tile, DoubleArrayTile}
+import geotrellis.spark._
+import geotrellis.util.MethodExtensions
+import geotrellis.vector.Point
+
+import org.apache.spark.rdd.RDD
+import org.apache.spark.SparkContext
+
+
+abstract class RDDCostDistanceMethods[K: (? => SpatialKey), V: (? => Tile)]
+    extends MethodExtensions[RDD[(K, V)] with Metadata[TileLayerMetadata[K]]] {
+
+  def costdistance(
+    points: Seq[Point]
+  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] with Metadata[TileLayerMetadata[K]] =
+    IterativeCostDistance(self, points)
+
+  def costdistance(
+    points: Seq[Point],
+    maxCost: Double
+  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] with Metadata[TileLayerMetadata[K]] =
+    IterativeCostDistance(self, points, maxCost)
+
+}
diff --git a/spark/src/main/scala/geotrellis/spark/package.scala b/spark/src/main/scala/geotrellis/spark/package.scala
index aa1d47a561..22614f9bb7 100644
--- a/spark/src/main/scala/geotrellis/spark/package.scala
+++ b/spark/src/main/scala/geotrellis/spark/package.scala
@@ -36,6 +36,7 @@ import java.time.Instant
 
 package object spark
     extends buffer.Implicits
+    with costdistance.Implicits
     with crop.Implicits
     with density.Implicits
     with equalization.Implicits
diff --git a/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 523ac44617..1ab55cfb0f 100644
--- a/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -1,5 +1,5 @@
 /*
- * Copyright 2016 Azavea
+ * Copyright 2017 Azavea
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.

From 2e26151da08e523adabeced10fd9ee8f629dd61e Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Sat, 18 Feb 2017 19:52:23 -0500
Subject: [PATCH 23/28] Cost-Distance Extension Method Tests

---
 .../RDDCostDistanceMethodsSpec.scala          | 62 +++++++++++++++++++
 1 file changed, 62 insertions(+)
 create mode 100644 spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala

diff --git a/spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala b/spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala
new file mode 100644
index 0000000000..190d471018
--- /dev/null
+++ b/spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala
@@ -0,0 +1,62 @@
+/*
+ * Copyright 2017 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.spark.costdistance
+
+import geotrellis.proj4.LatLng
+import geotrellis.raster._
+import geotrellis.spark._
+import geotrellis.spark.tiling.LayoutDefinition
+import geotrellis.vector._
+
+import org.apache.spark.rdd.RDD
+
+import org.scalatest._
+
+
+class RDDCostDistanceMethodsSpec extends FunSpec
+    with Matchers
+    with TestEnvironment {
+
+  val rdd: RDD[(SpatialKey, Tile)] with Metadata[TileLayerMetadata[SpatialKey]] = {
+    val tile = IntArrayTile(Array.fill[Int](25)(1), 5, 5)
+    val extent = Extent(0, 0, 10, 5)
+    val gridExtent = GridExtent(extent, 1, 1) // 10×5 pixels
+    val layoutDefinition = LayoutDefinition(gridExtent, 10, 5)
+    val bounds = Bounds(SpatialKey(0,0), SpatialKey(1,0))
+    val tileLayerMetadata = TileLayerMetadata(IntCellType, layoutDefinition, extent, LatLng, bounds)
+    val list = List((SpatialKey(0,0), tile), (SpatialKey(1,0), tile))
+    ContextRDD(sc.parallelize(list), tileLayerMetadata)
+  }
+
+  val points = List(Point(2.5+5.0, 2.5))
+
+  it("The costdistance Method Should Work (1/2)") {
+    val expected = IterativeCostDistance(rdd, points).collect.toList
+    val actual = rdd.costdistance(points).collect.toList
+
+    actual should be (expected)
+  }
+
+  it("The costdistance Method Should Work (2/2)") {
+    val resolution = IterativeCostDistance.computeResolution(rdd)
+    val expected = IterativeCostDistance(rdd, points, resolution).collect.toList
+    val actual = rdd.costdistance(points, resolution).collect.toList
+
+    actual should be (expected)
+  }
+
+}

From e26657ceda1035bfc6ab3b13dfbc087bd0b47cac Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Mon, 20 Feb 2017 15:02:06 -0500
Subject: [PATCH 24/28] Loosen RDD Type

---
 .../geotrellis/spark/costdistance/IterativeCostDistance.scala | 4 ++--
 .../spark/costdistance/RDDCostDistanceMethods.scala           | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 5fa2080c5f..a4c8ec6538 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -97,7 +97,7 @@ object IterativeCostDistance {
     friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]],
     points: Seq[Point],
     maxCost: Double = Double.PositiveInfinity
-  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] with Metadata[TileLayerMetadata[K]]= {
+  )(implicit sc: SparkContext): RDD[(K, Tile)] with Metadata[TileLayerMetadata[K]]= {
 
     val md = friction.metadata
     val mt = md.mapTransform
@@ -214,7 +214,7 @@ object IterativeCostDistance {
 
     // Construct return value and return it
     val metadata = TileLayerMetadata(DoubleCellType, md.layout, md.extent, md.crs, md.bounds)
-    val rdd = costs.map({ case (k, _, cost) => (k, cost) })
+    val rdd = costs.map({ case (k, _, cost) => (k, cost.asInstanceOf[Tile]) })
     ContextRDD(rdd, metadata)
   }
 
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala b/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
index ad3f5f2001..4bc792c34c 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
@@ -30,13 +30,13 @@ abstract class RDDCostDistanceMethods[K: (? => SpatialKey), V: (? => Tile)]
 
   def costdistance(
     points: Seq[Point]
-  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] with Metadata[TileLayerMetadata[K]] =
+  )(implicit sc: SparkContext): RDD[(K, Tile)] with Metadata[TileLayerMetadata[K]] =
     IterativeCostDistance(self, points)
 
   def costdistance(
     points: Seq[Point],
     maxCost: Double
-  )(implicit sc: SparkContext): RDD[(K, DoubleArrayTile)] with Metadata[TileLayerMetadata[K]] =
+  )(implicit sc: SparkContext): RDD[(K, Tile)] with Metadata[TileLayerMetadata[K]] =
     IterativeCostDistance(self, points, maxCost)
 
 }

From 640efef11592f21c9e2cefee447e8fe050d5c2ac Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Fri, 24 Feb 2017 09:06:24 -0500
Subject: [PATCH 25/28] Perform Cost-Distance Over Geometries

Previously, only points were accepted.
---
 .../costdistance/IterativeCostDistance.scala  | 63 +++++++++++++++----
 .../costdistance/RDDCostDistanceMethods.scala | 10 +--
 2 files changed, 55 insertions(+), 18 deletions(-)

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index a4c8ec6538..9acb25a99b 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -19,7 +19,10 @@ package geotrellis.spark.costdistance
 import geotrellis.proj4.LatLng
 import geotrellis.raster._
 import geotrellis.raster.costdistance.CostDistance
+import geotrellis.raster.rasterize.Rasterizer
 import geotrellis.spark._
+import geotrellis.spark.tiling._
+import geotrellis.util._
 import geotrellis.vector._
 
 import org.apache.log4j.Logger
@@ -86,16 +89,45 @@ object IterativeCostDistance {
     math.abs(meters / pixels)
   }
 
+  private def geometryToKeys[K: (? => SpatialKey)](
+    md: TileLayerMetadata[K],
+    g: Geometry
+  ) = {
+    val keys = mutable.ArrayBuffer.empty[SpatialKey]
+    val bounds = md.layout.mapTransform(g.envelope)
+
+    var row = bounds.rowMin; while (row <= bounds.rowMax) {
+      var col = bounds.colMin; while (col <= bounds.colMax) {
+        keys.append(SpatialKey(col, row))
+        col += 1
+      }
+      row += 1
+    }
+
+    keys.toList
+  }
+
+  private def geometryMap[K: (? => SpatialKey)](
+    md: TileLayerMetadata[K],
+    gs: Seq[Geometry]
+  ): Map[SpatialKey, Seq[Geometry]] = {
+    gs
+      .flatMap({ g => geometryToKeys(md, g).map({ k => (k, g) }) })
+      .groupBy(_._1)
+      .mapValues({ list => list.map({ case (_, v) => v }) })
+      .toMap
+  }
+
   /**
     * Perform the cost-distance computation.
     *
-    * @param  friction  The friction layer; pixels are in units of "seconds per meter"
-    * @param  points    The starting locations from-which to compute the cost of traveling
-    * @param  maxCost   The maximum cost before pruning a path (in units of "seconds")
+    * @param  friction    The friction layer; pixels are in units of "seconds per meter"
+    * @param  geometries  The starting locations from-which to compute the cost of traveling
+    * @param  maxCost     The maximum cost before pruning a path (in units of "seconds")
     */
   def apply[K: (? => SpatialKey), V: (? => Tile)](
     friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]],
-    points: Seq[Point],
+    geometries: Seq[Geometry],
     maxCost: Double = Double.PositiveInfinity
   )(implicit sc: SparkContext): RDD[(K, Tile)] with Metadata[TileLayerMetadata[K]]= {
 
@@ -115,6 +147,9 @@ object IterativeCostDistance {
     val accumulator = new ChangesAccumulator
     sc.register(accumulator)
 
+    // Index the input geometry by SpatialKey
+    val gs = sc.broadcast(geometryMap(md, geometries))
+
     // Create RDD of initial (empty) cost tiles and load the
     // accumulator with the starting values.
     var costs: RDD[(K, V, DoubleArrayTile)] = friction.map({ case (k, v) =>
@@ -124,15 +159,17 @@ object IterativeCostDistance {
       val rows = tile.rows
       val extent = mt(key)
       val rasterExtent = RasterExtent(extent, cols, rows)
-
-      points
-        .filter({ point => extent.contains(point) })
-        .map({ point =>
-          val col = rasterExtent.mapXToGrid(point.x)
-          val row = rasterExtent.mapYToGrid(point.y)
-          val friction = tile.getDouble(col, row)
-          val cost = (col, row, friction, 0.0)
-          accumulator.add((key, cost))
+      val options = Rasterizer.Options.DEFAULT
+
+      gs.value.getOrElse(key, List.empty[Geometry])
+        .filter({ geometry => extent.intersects(geometry) })
+        .foreach({ geometry =>
+          Rasterizer
+            .foreachCellByGeometry(geometry, rasterExtent, options)({ (col, row) =>
+              val friction = tile.getDouble(col, row)
+              val entry = (col, row, friction, 0.0)
+              accumulator.add((key, entry))
+            })
         })
 
       (k, v, CostDistance.generateEmptyCostTile(cols, rows))
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala b/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
index 4bc792c34c..bf8d6861c8 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/RDDCostDistanceMethods.scala
@@ -19,7 +19,7 @@ package geotrellis.spark.costdistance
 import geotrellis.raster.{Tile, DoubleArrayTile}
 import geotrellis.spark._
 import geotrellis.util.MethodExtensions
-import geotrellis.vector.Point
+import geotrellis.vector.Geometry
 
 import org.apache.spark.rdd.RDD
 import org.apache.spark.SparkContext
@@ -29,14 +29,14 @@ abstract class RDDCostDistanceMethods[K: (? => SpatialKey), V: (? => Tile)]
     extends MethodExtensions[RDD[(K, V)] with Metadata[TileLayerMetadata[K]]] {
 
   def costdistance(
-    points: Seq[Point]
+    geometries: Seq[Geometry]
   )(implicit sc: SparkContext): RDD[(K, Tile)] with Metadata[TileLayerMetadata[K]] =
-    IterativeCostDistance(self, points)
+    IterativeCostDistance(self, geometries)
 
   def costdistance(
-    points: Seq[Point],
+    geometries: Seq[Geometry],
     maxCost: Double
   )(implicit sc: SparkContext): RDD[(K, Tile)] with Metadata[TileLayerMetadata[K]] =
-    IterativeCostDistance(self, points, maxCost)
+    IterativeCostDistance(self, geometries, maxCost)
 
 }

From 94daa42a4c5e97d1bb9b3027ff4496f139363879 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Fri, 3 Mar 2017 10:14:00 -0500
Subject: [PATCH 26/28] TestEnvironment

---
 .../geotrellis/spark/costdistance/IterativeCostDistance.scala    | 1 +
 .../spark/costdistance/RDDCostDistanceMethodsSpec.scala          | 1 +
 2 files changed, 2 insertions(+)

diff --git a/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 1ab55cfb0f..0f6bf0bcc5 100644
--- a/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -19,6 +19,7 @@ package geotrellis.spark.costdistance
 import geotrellis.proj4.LatLng
 import geotrellis.raster._
 import geotrellis.spark._
+import geotrellis.spark.testkit.TestEnvironment
 import geotrellis.spark.tiling.LayoutDefinition
 import geotrellis.vector._
 
diff --git a/spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala b/spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala
index 190d471018..9a0a94e194 100644
--- a/spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala
+++ b/spark/src/test/scala/geotrellis/spark/costdistance/RDDCostDistanceMethodsSpec.scala
@@ -19,6 +19,7 @@ package geotrellis.spark.costdistance
 import geotrellis.proj4.LatLng
 import geotrellis.raster._
 import geotrellis.spark._
+import geotrellis.spark.testkit.TestEnvironment
 import geotrellis.spark.tiling.LayoutDefinition
 import geotrellis.vector._
 

From 4ed92a2e392f5482a76972df62116411e6278094 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Mon, 27 Mar 2017 07:24:08 -0400
Subject: [PATCH 27/28] Restore Original CostDistance.scala

---
 .../costdistance/CostDistanceSpec.scala       |  96 -----
 .../costdistance/SimpleCostDistanceSpec.scala | 210 +++++++++++
 .../raster/costdistance/CostDistance.scala    | 346 ++++++++++--------
 .../costdistance/SimpleCostDistance.scala     | 205 +++++++++++
 .../costdistance/IterativeCostDistance.scala  |  20 +-
 5 files changed, 621 insertions(+), 256 deletions(-)
 create mode 100644 raster-test/src/test/scala/geotrellis/raster/costdistance/SimpleCostDistanceSpec.scala
 create mode 100644 raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala

diff --git a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
index 01fa082f83..0dfc97224f 100644
--- a/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
+++ b/raster-test/src/test/scala/geotrellis/raster/costdistance/CostDistanceSpec.scala
@@ -109,102 +109,6 @@ class CostDistanceSpec extends FunSuite with RasterMatchers {
     }
   }
 
-  test("Boundary Callbacks") {
-    val n = NODATA
-    val N = Double.NaN
-    val size = 6
-
-    // Example from ESRI
-    // http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Cost_Distance_algorithm
-    val frictionTile = Array(
-      1,3,4,4,3,2,
-      4,6,2,3,7,6,
-      5,8,7,5,6,6,
-      1,4,5,n,5,1,
-      4,7,5,n,2,6,
-      1,2,2,1,3,4)
-
-    val actualLeft = Array.ofDim[Double](size)
-    val actualUp = Array.ofDim[Double](size)
-    val actualRight = Array.ofDim[Double](size)
-    val actualDown = Array.ofDim[Double](size)
-
-    val expectedLeft = Array(2.0, 4.5, 8.0, 5.0, 2.5, 0.0)
-    val expectedUp = Array(2.0, 0.0, 0.0, 4.0, 6.7, 9.2)
-    val expectedRight = Array(9.2, 13.1, 12.7, 9.2, 11.1, 10.5)
-    val expectedDown = Array(0.0, 1.5, 3.5, 5.0, 7.0, 10.5)
-
-    // Produce priority queue
-    val q = CostDistance.generateEmptyQueue(size, size)
-
-    // Insert starting points
-    val points = Seq((1,0), (2,0), (2,1), (0,5))
-    points.foreach({ case (col, row) =>
-      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
-      q.add(entry)
-    })
-
-    // Generate initial (empty) cost tile
-    val costTile = CostDistance.generateEmptyCostTile(size, size)
-
-    // Various callbacks to fill in the "actual" arrays
-    val edgeCb: CostDistance.EdgeCallback = { case (col: Int, row: Int, _: Double, cost: Double) =>
-      if (col == 0) actualLeft(row) = math.floor(cost*10)/10
-      if (col == size-1) actualRight(row) = math.floor(cost*10)/10
-      if (row == 0) actualUp(col) = math.floor(cost*10)/10
-      if (row == size-1) actualDown(col) = math.floor(cost*10)/10
-    }
-
-    CostDistance.compute(
-      frictionTile, costTile,
-      Double.PositiveInfinity, 1,
-      q, edgeCb)
-
-    actualLeft should be (expectedLeft)
-    actualUp should be (expectedUp)
-    actualRight should be (expectedRight)
-    actualDown should be (expectedDown)
-  }
-
-  test("Max Distance") {
-    val n = NODATA
-    val N = Double.NaN
-
-    // Example from ESRI
-    // http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Cost_Distance_algorithm
-    val costTile = Array(
-      1,3,4,4,3,2,
-      4,6,2,3,7,6,
-      5,8,7,5,6,6,
-      1,4,5,n,5,1,
-      4,7,5,n,2,6,
-      1,2,2,1,3,4)
-
-    val points = Seq(
-      (1,0),
-      (2,0),
-      (2,1),
-      (0,5))
-
-    val cd = CostDistance(costTile, points, 2.6)
-
-    val d = cd.toArrayDouble()
-
-    val expected = Array(
-      2.0, 0.0, 0.0,   N,   N,   N,
-        N,   N, 0.0, 2.5,   N,   N,
-        N,   N,   N,   N,   N,   N,
-        N,   N,   N,   N,   N,   N,
-      2.5,   N,   N,   N,   N,   N,
-      0.0, 1.5,   N,   N,   N,   N).map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
-
-    val strings = d.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
-
-    for(i <- 0 until strings.length) {
-      strings(i) should be (expected(i))
-    }
-  }
-
   def print(d: DoubleArrayTile):Unit = println(d.array.toList.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i)).grouped(d.cols).map(_.mkString(",")).mkString("\n"))
 
 }
diff --git a/raster-test/src/test/scala/geotrellis/raster/costdistance/SimpleCostDistanceSpec.scala b/raster-test/src/test/scala/geotrellis/raster/costdistance/SimpleCostDistanceSpec.scala
new file mode 100644
index 0000000000..8d8ff97ad6
--- /dev/null
+++ b/raster-test/src/test/scala/geotrellis/raster/costdistance/SimpleCostDistanceSpec.scala
@@ -0,0 +1,210 @@
+/*
+ * Copyright 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.costdistance
+
+import geotrellis.raster._
+import geotrellis.raster.testkit._
+
+import org.scalatest._
+
+import java.util.Locale
+import scala.language.implicitConversions
+
+
+class SimpleCostDistanceSpec extends FunSuite with RasterMatchers {
+  implicit def array2Tile(a: Array[Int]): Tile = {
+    val size = math.sqrt(a.length).toInt
+
+    IntArrayTile(a, size, size)
+  }
+
+  def asTile(a: Array[Int], cols: Int, rows: Int): Tile =
+    IntArrayTile(a, cols, rows)
+
+  test("ESRI example") {
+    val n = NODATA
+    val N = Double.NaN
+
+    // Example from ESRI
+    // http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Cost_Distance_algorithm
+    val costTile = Array(
+      1,3,4,4,3,2,
+      4,6,2,3,7,6,
+      5,8,7,5,6,6,
+      1,4,5,n,5,1,
+      4,7,5,n,2,6,
+      1,2,2,1,3,4)
+
+    val points = Seq(
+      (1,0),
+      (2,0),
+      (2,1),
+      (0,5))
+
+    val cd = SimpleCostDistance(costTile, points)
+
+    val d = cd.toArrayDouble()
+
+    val expected = Array(
+      2.0, 0.0, 0.0, 4.0, 6.7, 9.2,
+      4.5, 4.0, 0.0, 2.5, 7.5, 13.1,
+      8.0, 7.1, 4.5, 4.9, 8.9, 12.7,
+      5.0, 7.5, 10.5,  N, 10.6, 9.2,
+      2.5, 5.7, 6.4,   N, 7.1, 11.1,
+      0.0, 1.5, 3.5, 5.0, 7.0, 10.5).map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
+
+    val strings = d.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
+
+    for(i <- 0 until strings.length) {
+      strings(i) should be (expected(i))
+    }
+  }
+
+  test("GRASS example") {
+
+    // Example from GRASS
+    // http://grass.osgeo.org/grass64/manuals/r.cost.html
+    val costTile = asTile(Array(
+       2 , 2 , 1 , 1 , 5 , 5 , 5 ,
+       2 , 2 , 8 , 8 , 5 , 2 , 1 ,
+       7 , 1 , 1 , 8 , 2 , 2 , 2 ,
+       8 , 7 , 8 , 8 , 8 , 8 , 5 ,
+       8 , 8 , 1 , 1 , 5 , 3 , 9 ,
+       8 , 1 , 1 , 2 , 5 , 3 , 9), 7, 6)
+
+    val points = Seq((5,4))
+
+    val cd = SimpleCostDistance(costTile, points)
+
+    val d = cd.toArrayDouble()
+
+    val expected = Array(
+      22,21,21,20,17,15,14,
+      20,19,22,20,15,12,11,
+      22,18,17,18,13,11, 9,
+      21,14,13,12, 8, 6, 6,
+      16,13, 8, 7, 4, 0, 6,
+      14, 9, 8, 9, 6, 3, 8)
+
+    val values = d.toSeq.map(i => (i + 0.5).toInt)
+    // println(values.toSeq)
+    // println(expected.toSeq)
+
+    for(i <- 0 until values.length) {
+      values(i) should be (expected(i))
+    }
+  }
+
+  test("Boundary Callbacks") {
+    val n = NODATA
+    val N = Double.NaN
+    val size = 6
+
+    // Example from ESRI
+    // http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Cost_Distance_algorithm
+    val frictionTile = Array(
+      1,3,4,4,3,2,
+      4,6,2,3,7,6,
+      5,8,7,5,6,6,
+      1,4,5,n,5,1,
+      4,7,5,n,2,6,
+      1,2,2,1,3,4)
+
+    val actualLeft = Array.ofDim[Double](size)
+    val actualUp = Array.ofDim[Double](size)
+    val actualRight = Array.ofDim[Double](size)
+    val actualDown = Array.ofDim[Double](size)
+
+    val expectedLeft = Array(2.0, 4.5, 8.0, 5.0, 2.5, 0.0)
+    val expectedUp = Array(2.0, 0.0, 0.0, 4.0, 6.7, 9.2)
+    val expectedRight = Array(9.2, 13.1, 12.7, 9.2, 11.1, 10.5)
+    val expectedDown = Array(0.0, 1.5, 3.5, 5.0, 7.0, 10.5)
+
+    // Produce priority queue
+    val q = SimpleCostDistance.generateEmptyQueue(size, size)
+
+    // Insert starting points
+    val points = Seq((1,0), (2,0), (2,1), (0,5))
+    points.foreach({ case (col, row) =>
+      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
+      q.add(entry)
+    })
+
+    // Generate initial (empty) cost tile
+    val costTile = SimpleCostDistance.generateEmptyCostTile(size, size)
+
+    // Various callbacks to fill in the "actual" arrays
+    val edgeCb: SimpleCostDistance.EdgeCallback = { case (col: Int, row: Int, _: Double, cost: Double) =>
+      if (col == 0) actualLeft(row) = math.floor(cost*10)/10
+      if (col == size-1) actualRight(row) = math.floor(cost*10)/10
+      if (row == 0) actualUp(col) = math.floor(cost*10)/10
+      if (row == size-1) actualDown(col) = math.floor(cost*10)/10
+    }
+
+    SimpleCostDistance.compute(
+      frictionTile, costTile,
+      Double.PositiveInfinity, 1,
+      q, edgeCb)
+
+    actualLeft should be (expectedLeft)
+    actualUp should be (expectedUp)
+    actualRight should be (expectedRight)
+    actualDown should be (expectedDown)
+  }
+
+  test("Max Distance") {
+    val n = NODATA
+    val N = Double.NaN
+
+    // Example from ESRI
+    // http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Cost_Distance_algorithm
+    val costTile = Array(
+      1,3,4,4,3,2,
+      4,6,2,3,7,6,
+      5,8,7,5,6,6,
+      1,4,5,n,5,1,
+      4,7,5,n,2,6,
+      1,2,2,1,3,4)
+
+    val points = Seq(
+      (1,0),
+      (2,0),
+      (2,1),
+      (0,5))
+
+    val cd = SimpleCostDistance(costTile, points, 2.6)
+
+    val d = cd.toArrayDouble()
+
+    val expected = Array(
+      2.0, 0.0, 0.0,   N,   N,   N,
+        N,   N, 0.0, 2.5,   N,   N,
+        N,   N,   N,   N,   N,   N,
+        N,   N,   N,   N,   N,   N,
+      2.5,   N,   N,   N,   N,   N,
+      0.0, 1.5,   N,   N,   N,   N).map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
+
+    val strings = d.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i))
+
+    for(i <- 0 until strings.length) {
+      strings(i) should be (expected(i))
+    }
+  }
+
+  def print(d: DoubleArrayTile):Unit = println(d.array.toList.map(i => " %04.1f ".formatLocal(Locale.ENGLISH, i)).grouped(d.cols).map(_.mkString(",")).mkString("\n"))
+
+}
diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
index 69e1398174..867ab20313 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/CostDistance.scala
@@ -16,10 +16,10 @@
 
 package geotrellis.raster.costdistance
 
-import geotrellis.raster._
-
 import java.util.PriorityQueue
 
+import geotrellis.raster._
+
 
 /**
   * Object housing various functions related to Cost-Distance
@@ -27,179 +27,225 @@ import java.util.PriorityQueue
   */
 object CostDistance {
 
-  type Cost = (Int, Int, Double, Double) // column, row, friction, cost
-  type Q = PriorityQueue[Cost]
-  type EdgeCallback = (Cost => Unit)
-
-  /**
-    * NOP EdgeCallback
-    */
-  def nop(cost: Cost): Unit = {}
-
   /**
-    * Generate a Queue suitable for working with a tile of the given
-    * dimensions.
+    * Generate a Cost-Distance raster based on a set of starting points
+    * and a cost raster
+    *
+    * @param cost     Cost Tile (Int)
+    * @param points   List of starting points as tuples
+    *
+    * @note    Operation will only work with integer typed Cost Tiles (BitCellType, ByteConstantNoDataCellType, ShortConstantNoDataCellType, IntConstantNoDataCellType).
+    *          If a double typed Cost Tile (FloatConstantNoDataCellType, DoubleConstantNoDataCellType) is passed in, those costs will be rounded
+    *          to their floor integer values.
     *
-    * @param  cols  The number of columns of the friction tile
-    * @param  rows  The number of rows of the friction tile
     */
-  def generateEmptyQueue(cols: Int, rows: Int): Q = {
-    new PriorityQueue(
-      (cols*16 + rows*16), new java.util.Comparator[Cost] {
-        override def equals(a: Any) = a.equals(this)
-        def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
-      })
+  def apply(cost: Tile, points: Seq[(Int, Int)]): Tile = {
+    val (cols, rows) = cost.dimensions
+    val output = DoubleArrayTile.empty(cols, rows)
+
+    for((c, r) <- points)
+      output.setDouble(c, r, 0.0)
+
+    val pqueue = new PriorityQueue(
+        1000, new java.util.Comparator[Cost] {
+          override def equals(a: Any) = a.equals(this)
+          def compare(a: Cost, b: Cost) = a._3.compareTo(b._3)
+        })
+
+
+    for((c, r) <- points) {
+      calcNeighbors(c, r, cost, output, pqueue)
+    }
+
+    while (!pqueue.isEmpty) {
+      val (c, r, v) = pqueue.poll
+      if (v == output.getDouble(c, r)) calcNeighbors(c, r, cost, output, pqueue)
+    }
+
+    output
   }
 
   /**
-    * Generate an empty double-valued array tile of the correct
-    * dimensions.
+    * Predicate: are the given column and row within the bounds of the
+    * [[Tile]]?
     *
-    * @param  cols  The number of cols of the friction tile (and therefore the cost tile)
-    * @param  rows  The number of rows of the frition tile and cost tiles
+    * @param  c  The column
+    * @param  r  The row
     */
-  def generateEmptyCostTile(cols: Int, rows: Int): DoubleArrayTile =
-    DoubleArrayTile.empty(cols, rows)
+  private def isValid(c: Int, r: Int, cost: Tile): Boolean =
+    c >= 0 && r >= 0 && c < cost.cols && r < cost.rows
 
   /**
-    * Generate a cost-distance raster based on a set of starting
-    * points and a friction raster.  This is an implementation of the
-    * standard algorithm mentioned in the "previous work" section of [1].
-    *
-    * 1. Tomlin, Dana.
-    *    "Propagating radial waves of travel cost in a grid."
-    *    International Journal of Geographical Information Science 24.9 (2010): 1391-1413.
-    *
-    * @param  friction    Friction tile; pixels are interpreted as "second per meter"
-    * @param  points      List of starting points as tuples
-    * @param  maxCost     The maximum cost before pruning a path (in units of "seconds")
-    * @param  resolution  The resolution of the tiles (in units of "meters per pixel")
+    * A class encoding directions.
     */
-  def apply(
-    frictionTile: Tile,
-    points: Seq[(Int, Int)],
-    maxCost: Double = Double.PositiveInfinity,
-    resolution: Double = 1
-  ): DoubleArrayTile = {
-    val cols = frictionTile.cols
-    val rows = frictionTile.rows
-    val costTile = generateEmptyCostTile(cols, rows)
-    val q: Q = generateEmptyQueue(cols, rows)
-
-    points.foreach({ case (col, row) =>
-      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
-      q.add(entry)
-    })
-
-    compute(frictionTile, costTile, maxCost, resolution, q, nop)
+  final class Dir(val dc: Int, val dr: Int) {
+    val diag = !(dc == 0 || dr == 0)
+
+    lazy val cornerOffsets = (dc, dr) match {
+      case (-1, -1) => Array((0, -1), (-1, 0))
+      case ( 1, -1) => Array((0, -1), ( 1, 0))
+      case (-1, 1) => Array((0, 1), (-1, 0))
+      case ( 1, 1) => Array((0, 1), ( 1, 0))
+      case _ => Array[(Int, Int)]()
+    }
+
+    def apply(c: Int, r: Int) = (c + dc, r + dr)
+
+    lazy val unitDistance = if (diag) 1.41421356237 else 1.0
   }
 
+  val dirs: Array[Dir] = Array(
+    (-1, -1), ( 0, -1), ( 1, -1),
+    (-1, 0),          ( 1, 0),
+    (-1, 1), ( 0, 1), ( 1, 1)).map { case (c, r) => new Dir(c, r) }
+
+
   /**
-    * Compute a cost tile.
+    * Input:
+    * (c, r) => Source cell
+    * (dc, dr) => Delta (direction)
+    * cost => Cost raster
+    * d => C - D output raster
     *
-    * @param  frictionTile    The friction tile
-    * @param  costTile        The tile that will contain the costs
-    * @param  maxCost         The maximum cost before pruning a path (in units of "seconds")
-    * @param  resolution      The resolution of the tiles (in units of "meters per pixel")
-    * @param  q               A priority queue of Cost objects (a.k.a. candidate paths)
-    * @param  edgeCallback    Called when a pixel on the edge of the tile is updated
+    * Output:
+    * List((c, r)) <- list of cells set
     */
-  def compute(
-    frictionTile: Tile,
-    costTile: DoubleArrayTile,
-    maxCost: Double, resolution: Double,
-    q: Q, edgeCallback: EdgeCallback
-  ): DoubleArrayTile = {
-    val cols = frictionTile.cols
-    val rows = frictionTile.rows
-
-    require(frictionTile.dimensions == costTile.dimensions)
-
-    def inTile(col: Int, row: Int): Boolean =
-      ((0 <= col && col < cols) && (0 <= row && row < rows))
-
-    def isPassable(f: Double): Boolean =
-      (isData(f) && 0.0 <= f)
-
-    def onEdge(col: Int, row: Int): Boolean =
-      ((col == 0) || (row == 0) || (col == cols-1) || (row == rows-1))
-
-    /**
-      * Given a location, an instantaneous cost at that neighboring
-      * location (friction), the cost to get to the neighboring
-      * location, and the distance from the neighboring pixel to this
-      * pixel, enqueue a candidate path to the present pixel.
-      *
-      * @param  col           The column of the given location
-      * @param  row           The row of the given location
-      * @param  friction1     The instantaneous cost (friction) at the neighboring location
-      * @param  cost          The length of the best-known path from a source to the neighboring location
-      * @param  distance      The distance from the neighboring location to this location
-      */
-    @inline def enqueueNeighbor(
-      col: Int, row: Int, friction1: Double, neighborCost: Double,
-      distance: Double = 1.0
-    ): Unit = {
-      // If the location is inside of the tile ...
-      if (inTile(col, row)) {
-        val friction2 = frictionTile.getDouble(col, row)
-        val currentCost = costTile.getDouble(col, row)
-
-        // ... and if the location is passable ...
-        if (isPassable(friction2)) {
-          val step = resolution * distance * (friction1 + friction2) / 2.0
-          val candidateCost = neighborCost + step
-          val entry = (col, row, friction2, candidateCost)
-
-          // ... and the candidate cost is less than the maximum cost ...
-          if (candidateCost <= maxCost) {
-            // ... and the candidate is a possible improvement ...
-            if (!isData(currentCost) || candidateCost < currentCost) {
-              costTile.setDouble(col, row, candidateCost) // then increase lower bound on pixel,
-              q.add(entry) // and enqueue candidate for future processing
+  private def calcCostCell(c: Int, r: Int, dir: Dir, cost: Tile, d: DoubleArrayTile) = {
+    val cr = dir(c, r)
+
+    if (isValid(cr._1, cr._2, cost)) {
+      val prev = d.getDouble(cr._1, cr._2)
+      if (prev == 0.0) { // This is a source cell, don't override and shortcircuit early
+        None
+      } else {
+        val source = d.getDouble(c, r)
+
+        // Previous value could be NODATA
+        val prevCost = if (isNoData(prev)) java.lang.Double.MAX_VALUE else prev
+
+        var curMinCost = Double.MaxValue
+
+        // Simple solution
+        val baseCostOpt = calcCost(c, r, dir, cost)
+        if (baseCostOpt.isDefined) {
+          curMinCost = source + baseCostOpt.get
+        }
+
+        // Alternative solutions (going around the corner)
+        // Generally you can check diags directly:
+        // +---+---+---+
+        // | a | b | c |
+        // +---+---+---+
+        // | d | e | f |
+        // +---+---+---+
+        // | g | h | i |
+        // +---+---+---+
+        //
+        // For instance, "eg" can be computed directly
+        // but it turns out "eg" can be more expensive
+        // than "edg" or "ehg" so we compute those right
+        // now just in case
+        val cornerOffsets = dir.cornerOffsets
+        val l = cornerOffsets.length
+        var z = 0
+        while(z < l) {
+          val p = cornerOffsets(z)
+          val c1 = p._1 + c
+          val r1 = p._2 + r
+          val cost1 = calcCost(c, r, c1, r1, cost)
+          if (cost1.isDefined) {
+            val cost2 = calcCost(c1, r1, dir(c, r), cost)
+            if (cost2.isDefined) {
+              curMinCost = math.min(curMinCost, source + cost1.get + cost2.get)
             }
           }
+          z += 1
+        }
+
+        if (curMinCost == Double.MaxValue) {
+          None // Possible all nodata values
+        } else {
+          if (curMinCost < prevCost) {
+            d.setDouble(cr._1, cr._2, curMinCost)
+
+            Some((cr._1, cr._2, curMinCost))
+          } else {
+            None
+          }
         }
       }
+    } else {
+      None
     }
+  }
 
-    /**
-      * Process the candidate path on the top of the queue.
-      *
-      * @param  frictionTile  The friction tile
-      * @param  costTile      The cost tile
-      * @param  q             The priority queue of candidate paths
-      */
-    def processNext(): Unit = {
-      val entry: Cost = q.poll
-      val (col, row, friction1, candidateCost) = entry
-      val currentCost =
-        if (inTile(col, row))
-          costTile.getDouble(col, row)
-        else
-          Double.NaN
-
-      // If the candidate path is an improvement ...
-      if (!isData(currentCost) || candidateCost <= currentCost) {
-        if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
-        if (onEdge(col, row)) edgeCallback(entry)
-
-        // Compute candidate costs for neighbors and enqueue them
-        if (isPassable(friction1)) {
-          enqueueNeighbor(col-1, row+0, friction1, candidateCost)
-          enqueueNeighbor(col+1, row+0, friction1, candidateCost)
-          enqueueNeighbor(col+0, row+1, friction1, candidateCost)
-          enqueueNeighbor(col+0, row-1, friction1, candidateCost)
-          enqueueNeighbor(col-1, row-1, friction1, candidateCost, math.sqrt(2))
-          enqueueNeighbor(col-1, row+1, friction1, candidateCost, math.sqrt(2))
-          enqueueNeighbor(col+1, row-1, friction1, candidateCost, math.sqrt(2))
-          enqueueNeighbor(col+1, row+1, friction1, candidateCost, math.sqrt(2))
-        }
+  type Cost = (Int, Int, Double)
+
+  private def calcNeighbors(c: Int, r: Int, cost: Tile, d: DoubleArrayTile, p: PriorityQueue[Cost]) {
+    val l = dirs.length
+    var z = 0
+
+    while(z < l) {
+      val opt = calcCostCell(c, r, dirs(z), cost, d)
+      if (opt.isDefined) {
+        p.add(opt.get)
       }
+      z += 1
     }
+  }
 
-    while (!q.isEmpty) processNext
+  private def factor(c: Int, r: Int, c1: Int, r1: Int) = if (c == c1 || r == r1) 1.0 else 1.41421356237
 
-    costTile
+  private def safeGet(c: Int, r: Int, cost: Tile): IOption = IOption(cost.get(c, r))
+
+  private def calcCost(c: Int, r: Int, c2: Int, r2: Int, cost: Tile): DOption = {
+    val cost1 = safeGet(c, r, cost)
+    val cost2 = safeGet(c2, r2, cost)
+
+    if (cost1.isDefined && cost2.isDefined) {
+      DOption(factor(c, r, c2, r2) * (cost1.get + cost2.get) / 2.0)
+    } else {
+      DOption.None
+    }
   }
+
+  private def calcCost(c: Int, r: Int, cr2: (Int, Int), cost: Tile): DOption =
+    calcCost(c, r, cr2._1, cr2._2, cost)
+
+  private def calcCost(c: Int, r: Int, dir: Dir, cost: Tile): DOption =
+    calcCost(c, r, dir(c, r), cost)
+}
+
+/**
+  * Represents an optional integer using 'Tile NODATA' as a flag
+  */
+private [costdistance]
+class IOption(val v: Int) extends AnyVal {
+  def map(f: Int => Int) = if (isDefined) new IOption(f(v)) else this
+  def flatMap(f: Int => IOption) = if (isDefined) f(v) else this
+  def isDefined = isData(v)
+  def get = if (isDefined) v else sys.error("Get called on NODATA")
+}
+
+/**
+  * Represents an optional integer using 'Double.NaN' as a flag
+  */
+private [costdistance]
+class DOption(val v: Double) extends AnyVal {
+  def map(f: Double => Double) = if (isDefined) new DOption(f(v)) else this
+  def flatMap(f: Double => DOption) = if (isDefined) f(v) else this
+  def isDefined = isData(v)
+  def get = if (isDefined) v else sys.error("Get called on NaN")
+}
+
+private [costdistance]
+object IOption {
+  val None = new IOption(NODATA)
+  def apply(v: Int) = new IOption(v)
+}
+
+private [costdistance]
+object DOption {
+  val None = new DOption(Double.NaN)
+  def apply(v: Double) = new DOption(v)
 }
diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala
new file mode 100644
index 0000000000..800c70916a
--- /dev/null
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala
@@ -0,0 +1,205 @@
+/*
+ * Copyright 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.costdistance
+
+import geotrellis.raster._
+
+import java.util.PriorityQueue
+
+
+/**
+  * Object housing various functions related to Cost-Distance
+  * computations.
+  */
+object SimpleCostDistance {
+
+  type Cost = (Int, Int, Double, Double) // column, row, friction, cost
+  type Q = PriorityQueue[Cost]
+  type EdgeCallback = (Cost => Unit)
+
+  /**
+    * NOP EdgeCallback
+    */
+  def nop(cost: Cost): Unit = {}
+
+  /**
+    * Generate a Queue suitable for working with a tile of the given
+    * dimensions.
+    *
+    * @param  cols  The number of columns of the friction tile
+    * @param  rows  The number of rows of the friction tile
+    */
+  def generateEmptyQueue(cols: Int, rows: Int): Q = {
+    new PriorityQueue(
+      (cols*16 + rows*16), new java.util.Comparator[Cost] {
+        override def equals(a: Any) = a.equals(this)
+        def compare(a: Cost, b: Cost) = a._4.compareTo(b._4)
+      })
+  }
+
+  /**
+    * Generate an empty double-valued array tile of the correct
+    * dimensions.
+    *
+    * @param  cols  The number of cols of the friction tile (and therefore the cost tile)
+    * @param  rows  The number of rows of the frition tile and cost tiles
+    */
+  def generateEmptyCostTile(cols: Int, rows: Int): DoubleArrayTile =
+    DoubleArrayTile.empty(cols, rows)
+
+  /**
+    * Generate a cost-distance raster based on a set of starting
+    * points and a friction raster.  This is an implementation of the
+    * standard algorithm mentioned in the "previous work" section of [1].
+    *
+    * 1. Tomlin, Dana.
+    *    "Propagating radial waves of travel cost in a grid."
+    *    International Journal of Geographical Information Science 24.9 (2010): 1391-1413.
+    *
+    * @param  friction    Friction tile; pixels are interpreted as "second per meter"
+    * @param  points      List of starting points as tuples
+    * @param  maxCost     The maximum cost before pruning a path (in units of "seconds")
+    * @param  resolution  The resolution of the tiles (in units of "meters per pixel")
+    */
+  def apply(
+    frictionTile: Tile,
+    points: Seq[(Int, Int)],
+    maxCost: Double = Double.PositiveInfinity,
+    resolution: Double = 1
+  ): DoubleArrayTile = {
+    val cols = frictionTile.cols
+    val rows = frictionTile.rows
+    val costTile = generateEmptyCostTile(cols, rows)
+    val q: Q = generateEmptyQueue(cols, rows)
+
+    points.foreach({ case (col, row) =>
+      val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
+      q.add(entry)
+    })
+
+    compute(frictionTile, costTile, maxCost, resolution, q, nop)
+  }
+
+  /**
+    * Compute a cost tile.
+    *
+    * @param  frictionTile    The friction tile
+    * @param  costTile        The tile that will contain the costs
+    * @param  maxCost         The maximum cost before pruning a path (in units of "seconds")
+    * @param  resolution      The resolution of the tiles (in units of "meters per pixel")
+    * @param  q               A priority queue of Cost objects (a.k.a. candidate paths)
+    * @param  edgeCallback    Called when a pixel on the edge of the tile is updated
+    */
+  def compute(
+    frictionTile: Tile,
+    costTile: DoubleArrayTile,
+    maxCost: Double, resolution: Double,
+    q: Q, edgeCallback: EdgeCallback
+  ): DoubleArrayTile = {
+    val cols = frictionTile.cols
+    val rows = frictionTile.rows
+
+    require(frictionTile.dimensions == costTile.dimensions)
+
+    def inTile(col: Int, row: Int): Boolean =
+      ((0 <= col && col < cols) && (0 <= row && row < rows))
+
+    def isPassable(f: Double): Boolean =
+      (isData(f) && 0.0 <= f)
+
+    def onEdge(col: Int, row: Int): Boolean =
+      ((col == 0) || (row == 0) || (col == cols-1) || (row == rows-1))
+
+    /**
+      * Given a location, an instantaneous cost at that neighboring
+      * location (friction), the cost to get to the neighboring
+      * location, and the distance from the neighboring pixel to this
+      * pixel, enqueue a candidate path to the present pixel.
+      *
+      * @param  col           The column of the given location
+      * @param  row           The row of the given location
+      * @param  friction1     The instantaneous cost (friction) at the neighboring location
+      * @param  cost          The length of the best-known path from a source to the neighboring location
+      * @param  distance      The distance from the neighboring location to this location
+      */
+    @inline def enqueueNeighbor(
+      col: Int, row: Int, friction1: Double, neighborCost: Double,
+      distance: Double = 1.0
+    ): Unit = {
+      // If the location is inside of the tile ...
+      if (inTile(col, row)) {
+        val friction2 = frictionTile.getDouble(col, row)
+        val currentCost = costTile.getDouble(col, row)
+
+        // ... and if the location is passable ...
+        if (isPassable(friction2)) {
+          val step = resolution * distance * (friction1 + friction2) / 2.0
+          val candidateCost = neighborCost + step
+          val entry = (col, row, friction2, candidateCost)
+
+          // ... and the candidate cost is less than the maximum cost ...
+          if (candidateCost <= maxCost) {
+            // ... and the candidate is a possible improvement ...
+            if (!isData(currentCost) || candidateCost < currentCost) {
+              costTile.setDouble(col, row, candidateCost) // then increase lower bound on pixel,
+              q.add(entry) // and enqueue candidate for future processing
+            }
+          }
+        }
+      }
+    }
+
+    /**
+      * Process the candidate path on the top of the queue.
+      *
+      * @param  frictionTile  The friction tile
+      * @param  costTile      The cost tile
+      * @param  q             The priority queue of candidate paths
+      */
+    def processNext(): Unit = {
+      val entry: Cost = q.poll
+      val (col, row, friction1, candidateCost) = entry
+      val currentCost =
+        if (inTile(col, row))
+          costTile.getDouble(col, row)
+        else
+          Double.NaN
+
+      // If the candidate path is an improvement ...
+      if (!isData(currentCost) || candidateCost <= currentCost) {
+        if (inTile(col, row)) costTile.setDouble(col, row, candidateCost)
+        if (onEdge(col, row)) edgeCallback(entry)
+
+        // Compute candidate costs for neighbors and enqueue them
+        if (isPassable(friction1)) {
+          enqueueNeighbor(col-1, row+0, friction1, candidateCost)
+          enqueueNeighbor(col+1, row+0, friction1, candidateCost)
+          enqueueNeighbor(col+0, row+1, friction1, candidateCost)
+          enqueueNeighbor(col+0, row-1, friction1, candidateCost)
+          enqueueNeighbor(col-1, row-1, friction1, candidateCost, math.sqrt(2))
+          enqueueNeighbor(col-1, row+1, friction1, candidateCost, math.sqrt(2))
+          enqueueNeighbor(col+1, row-1, friction1, candidateCost, math.sqrt(2))
+          enqueueNeighbor(col+1, row+1, friction1, candidateCost, math.sqrt(2))
+        }
+      }
+    }
+
+    while (!q.isEmpty) processNext
+
+    costTile
+  }
+}
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 9acb25a99b..4d5dd514b4 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -18,7 +18,7 @@ package geotrellis.spark.costdistance
 
 import geotrellis.proj4.LatLng
 import geotrellis.raster._
-import geotrellis.raster.costdistance.CostDistance
+import geotrellis.raster.costdistance.SimpleCostDistance
 import geotrellis.raster.rasterize.Rasterizer
 import geotrellis.spark._
 import geotrellis.spark.tiling._
@@ -48,7 +48,7 @@ import scala.collection.mutable
   */
 object IterativeCostDistance {
 
-  type KeyCostPair = (SpatialKey, CostDistance.Cost)
+  type KeyCostPair = (SpatialKey, SimpleCostDistance.Cost)
   type Changes = mutable.ArrayBuffer[KeyCostPair]
 
   val logger = Logger.getLogger(IterativeCostDistance.getClass)
@@ -172,7 +172,7 @@ object IterativeCostDistance {
             })
         })
 
-      (k, v, CostDistance.generateEmptyCostTile(cols, rows))
+      (k, v, SimpleCostDistance.generateEmptyCostTile(cols, rows))
     }).persist(StorageLevel.MEMORY_AND_DISK_SER)
 
     costs.count
@@ -180,7 +180,7 @@ object IterativeCostDistance {
     // Repeatedly map over the RDD of cost tiles until no more changes
     // occur on the periphery of any tile.
     do {
-      val _changes: Map[SpatialKey, Seq[CostDistance.Cost]] =
+      val _changes: Map[SpatialKey, Seq[SimpleCostDistance.Cost]] =
         accumulator.value
           .groupBy(_._1)
           .map({ case (k, list) => (k, list.map({ case (_, v) => v })) })
@@ -199,20 +199,20 @@ object IterativeCostDistance {
         val keyRow = key._2
         val frictionTileCols = frictionTile.cols
         val frictionTileRows = frictionTile.rows
-        val localChanges: Option[Seq[CostDistance.Cost]] = changes.value.get(key)
+        val localChanges: Option[Seq[SimpleCostDistance.Cost]] = changes.value.get(key)
 
         localChanges match {
           case Some(localChanges) => {
-            val q: CostDistance.Q = {
-              val q = CostDistance.generateEmptyQueue(frictionTileCols, frictionTileRows)
-              localChanges.foreach({ (entry: CostDistance.Cost) => q.add(entry) })
+            val q: SimpleCostDistance.Q = {
+              val q = SimpleCostDistance.generateEmptyQueue(frictionTileCols, frictionTileRows)
+              localChanges.foreach({ (entry: SimpleCostDistance.Cost) => q.add(entry) })
               q
             }
 
-            val newCostTile = CostDistance.compute(
+            val newCostTile = SimpleCostDistance.compute(
               frictionTile, oldCostTile,
               maxCost, resolution,
-              q, { (entry: CostDistance.Cost) =>
+              q, { (entry: SimpleCostDistance.Cost) =>
                 val (col, row, f, c) = entry
                 if (col == 0 && (minKeyCol <= keyCol-1)) // left
                   accumulator.add((SpatialKey(keyCol-1, keyRow), (frictionTileCols, row, f, c)))

From 30fa62f4f29f26d92b1d256bd182443684e52798 Mon Sep 17 00:00:00 2001
From: James McClain <jmcclain@azavea.com>
Date: Mon, 27 Mar 2017 07:45:26 -0400
Subject: [PATCH 28/28] Other Changes

---
 .../costdistance/SimpleCostDistance.scala     | 23 ++++++-----
 .../costdistance/IterativeCostDistance.scala  | 39 +++++++++----------
 2 files changed, 30 insertions(+), 32 deletions(-)

diff --git a/raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala b/raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala
index 800c70916a..e3295be2c2 100644
--- a/raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala
+++ b/raster/src/main/scala/geotrellis/raster/costdistance/SimpleCostDistance.scala
@@ -64,16 +64,16 @@ object SimpleCostDistance {
   /**
     * Generate a cost-distance raster based on a set of starting
     * points and a friction raster.  This is an implementation of the
-    * standard algorithm mentioned in the "previous work" section of [1].
+    * standard algorithm cited in the "previous work" section of [1].
     *
     * 1. Tomlin, Dana.
     *    "Propagating radial waves of travel cost in a grid."
     *    International Journal of Geographical Information Science 24.9 (2010): 1391-1413.
     *
-    * @param  friction    Friction tile; pixels are interpreted as "second per meter"
-    * @param  points      List of starting points as tuples
-    * @param  maxCost     The maximum cost before pruning a path (in units of "seconds")
-    * @param  resolution  The resolution of the tiles (in units of "meters per pixel")
+    * @param  frictionTile  Friction tile; pixels are interpreted as "second per meter"
+    * @param  points        List of starting points as tuples
+    * @param  maxCost       The maximum cost before pruning a path (in units of "seconds")
+    * @param  resolution    The resolution of the tiles (in units of "meters per pixel")
     */
   def apply(
     frictionTile: Tile,
@@ -86,10 +86,12 @@ object SimpleCostDistance {
     val costTile = generateEmptyCostTile(cols, rows)
     val q: Q = generateEmptyQueue(cols, rows)
 
-    points.foreach({ case (col, row) =>
+    var i = 0; while (i < points.length) {
+      val (col, row) = points(i)
       val entry = (col, row, frictionTile.getDouble(col, row), 0.0)
       q.add(entry)
-    })
+      i += 1
+    }
 
     compute(frictionTile, costTile, maxCost, resolution, q, nop)
   }
@@ -133,6 +135,7 @@ object SimpleCostDistance {
       * @param  col           The column of the given location
       * @param  row           The row of the given location
       * @param  friction1     The instantaneous cost (friction) at the neighboring location
+      * @param  neighborCost  The cost of the neighbor
       * @param  cost          The length of the best-known path from a source to the neighboring location
       * @param  distance      The distance from the neighboring location to this location
       */
@@ -149,12 +152,12 @@ object SimpleCostDistance {
         if (isPassable(friction2)) {
           val step = resolution * distance * (friction1 + friction2) / 2.0
           val candidateCost = neighborCost + step
-          val entry = (col, row, friction2, candidateCost)
 
           // ... and the candidate cost is less than the maximum cost ...
           if (candidateCost <= maxCost) {
             // ... and the candidate is a possible improvement ...
             if (!isData(currentCost) || candidateCost < currentCost) {
+              val entry = (col, row, friction2, candidateCost)
               costTile.setDouble(col, row, candidateCost) // then increase lower bound on pixel,
               q.add(entry) // and enqueue candidate for future processing
             }
@@ -165,10 +168,6 @@ object SimpleCostDistance {
 
     /**
       * Process the candidate path on the top of the queue.
-      *
-      * @param  frictionTile  The friction tile
-      * @param  costTile      The cost tile
-      * @param  q             The priority queue of candidate paths
       */
     def processNext(): Unit = {
       val entry: Cost = q.poll
diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
index 4d5dd514b4..d68fb92334 100644
--- a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
+++ b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -79,9 +79,7 @@ object IterativeCostDistance {
   ) = {
     val md = friction.metadata
     val mt = md.mapTransform
-    val kv = friction.first
-    val key = implicitly[SpatialKey](kv._1)
-    val tile = implicitly[Tile](kv._2)
+    val (key: SpatialKey, tile: Tile) = friction.first
     val extent = mt(key).reproject(md.crs, LatLng)
     val degrees = extent.xmax - extent.xmin
     val meters = degrees * (6378137 * 2.0 * math.Pi) / 360.0
@@ -115,7 +113,6 @@ object IterativeCostDistance {
       .flatMap({ g => geometryToKeys(md, g).map({ k => (k, g) }) })
       .groupBy(_._1)
       .mapValues({ list => list.map({ case (_, v) => v }) })
-      .toMap
   }
 
   /**
@@ -136,13 +133,16 @@ object IterativeCostDistance {
     val resolution = computeResolution(friction)
     logger.debug(s"Computed resolution: $resolution meters/pixel")
 
-    val bounds = friction.metadata.bounds.asInstanceOf[KeyBounds[K]]
-    val minKey = implicitly[SpatialKey](bounds.minKey)
-    val minKeyCol = minKey._1
-    val minKeyRow = minKey._2
-    val maxKey = implicitly[SpatialKey](bounds.maxKey)
-    val maxKeyCol = maxKey._1
-    val maxKeyRow = maxKey._2
+    val bounds = friction.metadata.bounds match {
+      case b: KeyBounds[K] => b
+      case _ => throw new Exception
+    }
+    val minKey: SpatialKey = bounds.minKey
+    val minKeyCol = minKey.col
+    val minKeyRow = minKey.row
+    val maxKey: SpatialKey = bounds.maxKey
+    val maxKeyCol = maxKey.col
+    val maxKeyRow = maxKey.row
 
     val accumulator = new ChangesAccumulator
     sc.register(accumulator)
@@ -153,8 +153,8 @@ object IterativeCostDistance {
     // Create RDD of initial (empty) cost tiles and load the
     // accumulator with the starting values.
     var costs: RDD[(K, V, DoubleArrayTile)] = friction.map({ case (k, v) =>
-      val key = implicitly[SpatialKey](k)
-      val tile = implicitly[Tile](v)
+      val key: SpatialKey = k
+      val tile: Tile = v
       val cols = tile.cols
       val rows = tile.rows
       val extent = mt(key)
@@ -184,7 +184,6 @@ object IterativeCostDistance {
         accumulator.value
           .groupBy(_._1)
           .map({ case (k, list) => (k, list.map({ case (_, v) => v })) })
-          .toMap
       val changes = sc.broadcast(_changes)
       logger.debug(s"At least ${changes.value.size} changed tiles")
 
@@ -193,10 +192,10 @@ object IterativeCostDistance {
       val previous = costs
 
       costs = previous.map({ case (k, v, oldCostTile) =>
-        val key = implicitly[SpatialKey](k)
-        val frictionTile = implicitly[Tile](v)
-        val keyCol = key._1
-        val keyRow = key._2
+        val key: SpatialKey = k
+        val frictionTile: Tile = v
+        val keyCol = key.col
+        val keyRow = key.row
         val frictionTileCols = frictionTile.cols
         val frictionTileRows = frictionTile.rows
         val localChanges: Option[Seq[SimpleCostDistance.Cost]] = changes.value.get(key)
@@ -247,11 +246,11 @@ object IterativeCostDistance {
 
       costs.count
       previous.unpersist()
-    } while (accumulator.value.size > 0)
+    } while (accumulator.value.nonEmpty)
 
     // Construct return value and return it
     val metadata = TileLayerMetadata(DoubleCellType, md.layout, md.extent, md.crs, md.bounds)
-    val rdd = costs.map({ case (k, _, cost) => (k, cost.asInstanceOf[Tile]) })
+    val rdd = costs.map({ case (k, _, cost) => (k, cost: Tile) })
     ContextRDD(rdd, metadata)
   }