Use Bresenham-Inspired Approach For Line Drawing

locationtech · Dec 6, 2016 · 94aed1d · 94aed1d
1 parent 8beca2c
commit 94aed1d
Show file tree

Hide file tree

Showing 2 changed files with 145 additions and 22 deletions.
diff --git a/...er-test/src/test/scala/geotrellis/raster/rasterize/polygon/FractionalRasterizerSpec.scala b/...er-test/src/test/scala/geotrellis/raster/rasterize/polygon/FractionalRasterizerSpec.scala
@@ -91,5 +91,57 @@ class FractionalRasterizerSpec extends FunSpec with Matchers {
       round(actual * 1000000) should be (round(expected * 1000000))
     }
 
+    it("should efficiently handle a long diagonal line (1/4)") {
+      val re = RasterExtent(e, 30, 30)
+      val poly = Polygon(Point(0,0), Point(3,0), Point(3,3), Point(0,0))
+      var actual = 0.0
+      val expected = (30*30)/2.0
+
+      FractionalRasterizer.foreachCellByPolygon(poly, re) { (col: Int, row: Int, p: Double) =>
+        actual += p
+      }
+
+      round(actual * 1000000) should be (round(expected * 1000000))
+    }
+
+    it("should efficiently handle a long diagonal line (2/4)") {
+      val re = RasterExtent(e, 30, 30)
+      val poly = Polygon(Point(0,0), Point(3,0), Point(0,3), Point(0,0))
+      var actual = 0.0
+      val expected = (30*30)/2.0
+
+      FractionalRasterizer.foreachCellByPolygon(poly, re) { (col: Int, row: Int, p: Double) =>
+        actual += p
+      }
+
+      round(actual * 1000000) should be (round(expected * 1000000))
+    }
+
+    it("should efficiently handle a long diagonal line (3/4)") {
+      val re = RasterExtent(Extent(0, 0, 3, 3.1), 30, 30)
+      val poly = Polygon(Point(0,0), Point(3,0), Point(0,3), Point(0,0))
+      var actual = 0.0
+      val expected = 435.483871
+
+      FractionalRasterizer.foreachCellByPolygon(poly, re) { (col: Int, row: Int, p: Double) =>
+        actual += p
+      }
+
+      round(actual * 1000000) should be (round(expected * 1000000))
+    }
+
+    it("should efficiently handle a long diagonal line (4/4)") {
+      val re = RasterExtent(Extent(0, 0, 3.1, 3), 30, 30)
+      val poly = Polygon(Point(0,0), Point(3,0), Point(0,3), Point(0,0))
+      var actual = 0.0
+      val expected = 435.483871
+
+      FractionalRasterizer.foreachCellByPolygon(poly, re) { (col: Int, row: Int, p: Double) =>
+        actual += p
+      }
+
+      round(actual * 1000000) should be (round(expected * 1000000))
+    }
+
   }
 }
diff --git a/raster/src/main/scala/geotrellis/raster/rasterize/polygon/FractionalRasterizer.scala b/raster/src/main/scala/geotrellis/raster/rasterize/polygon/FractionalRasterizer.scala
@@ -8,7 +8,7 @@ import com.vividsolutions.jts.geom.Envelope
 import spire.syntax.cfor._
 
 import scala.collection.mutable
-import scala.math.{min, max, ceil, floor}
+import scala.math.{min, max, ceil, floor, abs}
 
 
 object FractionalRasterizer {
@@ -31,7 +31,7 @@ object FractionalRasterizer {
       val row2 = re.mapYToGridDouble(coord2.y)
 
       val segment =
-        if (row1 < row2) (col1, row1, col2, row2)
+        if (col1 < col2) (col1, row1, col2, row2)
         else (col2, row2, col1, row1)
 
       arrayBuffer += segment
@@ -50,7 +50,7 @@ object FractionalRasterizer {
         val row2 = re.mapYToGridDouble(coord2.y)
 
         val segment =
-          if (row1 < row2) (col1, row1, col2, row2)
+          if (col1 < col2) (col1, row1, col2, row2)
           else (col2, row2, col1, row1)
 
         arrayBuffer += segment
@@ -62,35 +62,106 @@ object FractionalRasterizer {
 
   private def renderEdge(
     edge: Segment,
+    re: RasterExtent,
     poly: Polygon,
     set: mutable.Set[(Int, Int)],
     fn: FractionCallback
   ): Unit = {
-    val xmin = floor(min(edge._1, edge._3)).toInt
-    val ymin = floor(min(edge._2, edge._4)).toInt
-    val xmax = ceil(max(edge._1, edge._3)).toInt
-    val ymax = ceil(max(edge._2, edge._4)).toInt
-
-    val envelope = Polygon(Point(xmin, ymin), Point(xmin, ymax), Point(xmax, ymax), Point(xmax, ymin), Point(xmin, ymin)).jtsGeom
+    // Screen coordinates
+    val (x0, y0, x1, y1) = edge
+    val xmin = min(x0, x1)
+    val ymin = min(y0, y1)
+    val xmax = max(x0, x1)
+    val ymax = max(y0, y1)
+    val m = (y1 - y0) / (x1 - x0)
+
+    // Integral screen coordinates
+    val xminint = floor(xmin).toInt
+    val yminint = floor(ymin).toInt
+    val xmaxint =  ceil(xmax).toInt
+    val ymaxint =  ceil(ymax).toInt
+
+    // Map coordinates
+    val xminmap = re.gridColToMap(xminint) - re.cellwidth/2
+    val yminmap = re.gridRowToMap(ymaxint) + re.cellheight/2
+    val xmaxmap = re.gridColToMap(xmaxint) - re.cellwidth/2
+    val ymaxmap = re.gridRowToMap(yminint) + re.cellheight/2
+
+    // Envelope around the edge (in map space)
+    val envelope = Polygon(
+      Point(xminmap, yminmap),
+      Point(xminmap, ymaxmap),
+      Point(xmaxmap, ymaxmap),
+      Point(xmaxmap, yminmap),
+      Point(xminmap, yminmap)
+    ).jtsGeom
+
+    // Intersection of envelope and polygon (in map space)
     val localPoly = poly.jtsGeom.intersection(envelope)
 
-    var x = xmin; while (x <= xmax) {
-      var y = ymin; while (y <= ymax) {
-        val pair = (x, y)
-        val pixel = Polygon(Point(x, y), Point(x, y+1), Point(x+1, y+1), Point(x+1, y), Point(x, y)).jtsGeom
-        val fraction = (pixel.intersection(localPoly)).getArea
-
-        if (fraction > 0.0) {
-          synchronized {
-            if (!set.contains(pair)) {
-              fn(x, y, fraction)
-              set += ((x, y))
+    if (abs(m) <= 1) { // Mostly horizontal
+      var x = xminint; while (x <= xmaxint) {
+        val _y = floor(m * (x + 0.5 - x0) + y0).toInt
+        var i = -1; while (i <= 1) {
+          val y = _y + i
+          val pair = (x, y)
+          val xmap0 = re.gridColToMap(x+0) - re.cellwidth/2
+          val xmap1 = re.gridColToMap(x+1) - re.cellwidth/2
+          val ymap0 = re.gridRowToMap(y+0) + re.cellheight/2
+          val ymap1 = re.gridRowToMap(y+1) + re.cellheight/2
+          val pixel = Polygon(
+            Point(xmap0, ymap0),
+            Point(xmap0, ymap1),
+            Point(xmap1, ymap1),
+            Point(xmap1, ymap0),
+            Point(xmap0, ymap0)
+          ).jtsGeom
+          val fraction = (pixel.intersection(localPoly)).getArea / pixel.getArea
+
+          if (fraction > 0.0) {
+            synchronized {
+              if (!set.contains(pair)) {
+                fn(x, y, fraction)
+                set += ((x, y))
+              }
+            }
+          }
+          i += 1
+        }
+        x += 1
+      }
+    } else { // Mostly vertical
+      val m = (x1 - x0) / (y1 - y0)
+      var y = yminint; while (y <= ymaxint) {
+        val _x = floor(m * (y + 0.5 - y0) + x0).toInt
+        var i = -1; while (i <= 1) {
+          val x = _x + i
+          val pair = (x, y)
+          val xmap0 = re.gridColToMap(x+0) - re.cellwidth/2
+          val xmap1 = re.gridColToMap(x+1) - re.cellwidth/2
+          val ymap0 = re.gridRowToMap(y+0) + re.cellheight/2
+          val ymap1 = re.gridRowToMap(y+1) + re.cellheight/2
+          val pixel = Polygon(
+            Point(xmap0, ymap0),
+            Point(xmap0, ymap1),
+            Point(xmap1, ymap1),
+            Point(xmap1, ymap0),
+            Point(xmap0, ymap0)
+          ).jtsGeom
+          val fraction = (pixel.intersection(localPoly)).getArea / pixel.getArea
+
+          if (fraction > 0.0) {
+            synchronized {
+              if (!set.contains(pair)) {
+                fn(x, y, fraction)
+                set += ((x, y))
+              }
             }
           }
+          i += 1
         }
         y += 1
       }
-      x += 1
     }
   }
 
@@ -103,7 +174,7 @@ object FractionalRasterizer {
 
     polygonToEdges(poly, re)
       .par
-      .foreach({ edge => renderEdge(edge, poly, seen, fn) })
+      .foreach({ edge => renderEdge(edge, re, poly, seen, fn) })
 
     PolygonRasterizer.foreachCellByPolygon(poly, re) {(col: Int, row: Int) =>
       val pair = (col, row)