Merge pull request #1545 from jpolchlo/feature/Voronoi

Voronoi diagrams and Delaunay triangulations

docs/vector/images/halfedge.fig

@@ -0,0 +1,58 @@
+#FIG 3.2  Produced by xfig version 3.2.5c
+Landscape
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+5 1 1 2 0 7 50 -1 -1 4.000 0 0 1 0 1395.000 5861.250 1845 5805 1620 6255 1170 6255
+	2 1 1.00 60.00 120.00
+5 1 0 2 0 7 50 -1 -1 4.000 0 0 1 0 5265.000 5343.750 4815 5400 5040 4950 5490 4950
+	2 1 1.00 60.00 120.00
+1 3 0 2 0 0 50 -1 20 4.000 1 0.0000 5355 5625 90 90 5355 5625 5355 5535
+1 3 0 1 0 0 50 -1 20 4.000 1 0.0000 1305 5580 90 90 1305 5580 1305 5490
+1 3 0 2 0 0 50 -1 20 4.000 1 0.0000 135 7695 90 90 135 7695 135 7605
+1 3 0 2 0 0 50 -1 20 4.000 1 0.0000 3375 10395 90 90 3375 10395 3375 10305
+1 3 0 2 0 0 50 -1 20 4.000 1 0.0000 6480 7830 90 90 6480 7830 6480 7740
+2 1 1 3 0 0 50 -1 20 6.000 0 0 7 1 0 2
+	2 1 1.00 60.00 120.00
+	 1305 5760 180 7560
+2 1 2 3 0 0 50 -1 20 3.000 0 0 7 1 0 2
+	2 1 1.00 60.00 120.00
+	 180 3600 1305 5400
+2 1 0 3 0 0 50 -1 20 0.000 0 0 7 1 0 2
+	2 1 1.00 60.00 120.00
+	 1530 5490 5130 5490
+2 1 1 3 0 0 50 -1 20 4.000 0 0 7 1 0 2
+	2 1 1.00 60.00 120.00
+	 5130 5715 1530 5715
+2 1 1 3 0 0 50 -1 20 4.000 0 0 7 1 0 2
+	2 1 1.00 60.00 120.00
+	 6480 7605 5355 5805
+2 1 2 3 0 0 50 -1 20 3.000 0 0 7 1 0 2
+	2 1 1.00 60.00 120.00
+	 5355 5400 6480 3600
+2 1 1 3 0 0 50 -1 20 8.000 0 0 -1 1 0 2
+	2 1 1.00 60.00 120.00
+	 3510 10260 6300 7920
+2 1 1 3 0 0 50 -1 20 8.000 0 0 -1 1 0 2
+	2 1 1.00 60.00 120.00
+	 270 7830 3240 10260
+3 0 0 2 0 7 50 -1 -1 4.000 0 1 0 6
+	2 1 1.00 60.00 120.00
+	 2655 5490 2655 5265 2205 5265 2205 5940 2655 5940 2655 5715
+	 0.000 1.000 1.000 1.000 1.000 0.000
+3 0 1 2 0 7 50 -1 -1 4.000 0 1 0 6
+	2 1 1.00 60.00 120.00
+	 3555 5715 3555 5940 4005 5940 4005 5265 3555 5265 3555 5490
+	 0.000 1.000 1.000 1.000 1.000 0.000
+4 0 0 50 -1 4 12 0.0000 0 150 495 1395 6615 f.next\001
+4 0 0 50 -1 4 12 0.0000 0 135 540 4860 4680 e.next\001
+4 0 0 50 -1 4 12 0.0000 0 150 60 3015 5940 f\001
+4 0 0 50 -1 4 12 0.0000 0 105 105 3015 5400 e\001
+4 0 0 50 -1 4 12 0.0000 0 195 720 1935 5220 e.flip = f\001
+4 0 0 50 -1 4 12 0.0000 0 195 720 3555 6165 f.flip = e\001
+4 0 0 50 -1 4 12 0.0000 0 135 840 5535 5580 a = e.vert\001
+4 0 0 50 -1 4 12 0.0000 0 150 795 360 5625 b = f.vert\001

docs/vector/images/stitch-triangles.fig

@@ -0,0 +1,32 @@
+#FIG 3.2  Produced by xfig version 3.2.5c
+Landscape
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+1 3 0 1 0 0 50 -1 20 0.000 1 0.0000 2250 4500 45 45 2250 4500 2295 4500
+1 3 0 1 0 0 50 -1 20 0.000 1 0.0000 3600 4050 45 45 3600 4050 3645 4050
+1 3 0 1 0 0 50 -1 20 0.000 1 0.0000 2700 3150 45 45 2700 3150 2745 3150
+1 3 0 1 0 0 50 -1 20 0.000 1 0.0000 5850 3150 45 45 5850 3150 5895 3150
+1 3 0 1 0 0 50 -1 20 0.000 1 0.0000 5400 2250 45 45 5400 2250 5445 2250
+1 3 0 1 0 0 50 -1 20 0.000 1 0.0000 4050 1350 45 45 4050 1350 4095 1350
+1 4 2 1 1 0 51 -1 -1 3.000 1 0.0000 2812 3937 796 796 2700 3150 2925 4725
+1 4 2 1 4 0 51 -1 -1 3.000 1 0.0000 3375 3825 2565 2565 4050 1350 2700 6300
+1 4 0 2 0 0 50 -1 -1 6.000 1 0.0000 4545 3195 1274 1274 3600 4050 5490 2340
+1 4 2 1 0 0 50 -1 -1 3.000 1 0.0000 5175 4725 1714 1714 3600 4050 6750 5400
+1 4 2 1 0 0 50 -1 -1 3.000 1 0.0000 4815 3825 1236 1236 3600 4050 6030 3600
+1 4 2 1 0 0 50 -1 -1 3.000 1 0.0000 4995 4275 1413 1413 3600 4050 6390 4500
+1 4 2 1 0 0 50 -1 -1 3.000 1 0.0000 4680 3375 1198 1198 3690 4050 5670 2700
+2 3 0 2 4 0 51 -1 -1 0.000 1 0 -1 0 0 4
+	 4050 1350 5850 3150 5400 2250 4050 1350
+2 3 0 2 1 0 51 -1 -1 0.000 1 0 -1 0 0 4
+	 2700 3150 2250 4500 3600 4050 2700 3150
+2 1 1 2 0 0 50 -1 -1 6.000 1 0 -1 0 0 2
+	 3600 4050 5850 3150
+3 0 0 1 0 0 50 -1 -1 4.000 0 1 0 3
+	1 1 1.00 60.00 120.00
+	 4860 7200 5130 3870 4230 2340
+	 0.000 1.000 0.000

docs/vector/voronoi.md

@@ -0,0 +1,26 @@
+# Voronoi Diagrams and Delaunay Triangulations
+
+Voronoi diagrams specify a partitioning of the plane into convex polygonal regions based on an input set of points, with the points being in one-to-one correspondence with the polygons.  Given the set of points `P`, let `p` be a point in that set; then `V(p)` is the Voronoi polygon corresponding to `p`. The interior of `V(p)` contains the part of the plane closer to `p` than any other point in `P`.
+
+To compute the Voronoi diagram, one actually goes about creating a triangulation of the input points called the _Delaunay triangulation_.  In this structure, all the points in `P` are vertices of a set of non-overlapping triangles that comprise the set `T(P)`.  Each triangle `t` in `T(P)` has the property that the unique circle passing through the vertices of `t` has no points of `P` in its interior.
+
+`T(P)` and `V(P)` (with the latter defined as `{V(p) | p in P}`) are _dual_ to each other in the following sense.  Each triangle in `T(P)` corresponds to a vertex in `V(P)` (a corner of some `V(p)`), each vertex in `T(P)` (which is just a point in `P`) corresponds to a polygon in `V(P)`, and each edge in `T(P)` corresponds to an edge in `V(P)`.  The vertices of `V(P)` are defined as the centers of the circumscribing circles of the triangles in `T(P)`.  These vertices are connected by edges such that if `t(p1)` and `t(p2)` share an edge, then the Voronoi vertices corresponding to those two triangles are connected in `V(P)`.  This duality between structures is important because it is much easier to compute the Delaunay triangulation and to take its dual than it is to directly compute the Voronoi diagram.
+
+This PR provides a divide-and-conquer approach to computing the Delaunay triangulation based on Guibas and Stolfi's 1985 ACM Transactions on Graphics paper.  In this case, the oldies are still the goodies, since only minor performance increases have been achieved over this baseline result---hardly worth the increase in complexity.
+
+The triangulation algorithm starts by ordering vertices according to (ascending) x-coordinate, breaking ties with the y-coordinate.  Duplicate vertices are ignored.  Then, the right and left halves of the vertices are recursively triangulated.  To stitch the resulting triangulations, we find a vertex from each of the left and right results so that the connecting edge is guaranteed to be in the convex hull of the merged triangulations; call this edge `base`.  Now, consider a circle that floats upwards and comes into contact with the endpoints of `base`.  This bubble will, by changing its radius, squeeze through the gap between the endpoints of `base`, and rise until it encounters another vertex.  By definition, this ball has no vertices of `P` in its interior, and so the three points on its boundary are the vertices of a Delaunay triangle.  See the following image for clarification:
+<center>
+  ![Merging Delaunay triangulations](images/stitch-triangles.png)
+</center>
+Here, we note that the red triangle's circumscribing ball contains vertices of the blue triangle, and so we will expect that the red triangle will not be part of the final triangulation.  As such the leftmost edge of the red triangle be deleted before the triangulation can be updated to include the triangle circumscribed by the solid black circle.
+
+This process continues, with the newly created edge serving as the new `base`, and the ball rising through until another vertex is encountered and so on, until the ball exits out the top and the triangulation is complete.
+
+## Mesh Representation
+
+The output of Delaunay triangulation and Voronoi diagrams are in the form of meshes represented by the half-edge structure.  These structures can be thought of as directed edges between points in space, where an edge needs two half-edges to complete its representation.  A half-edge, `e`, has three vital pieces of information: a vertex to which it points, `e.vert`; a pointer to its complementary half-edge, `e.flip`; and a pointer to the next half-edge in the polygon, `e.next`.  The following image might be useful:
+<center>
+    ![Half-edge structure](images/halfedge.png)
+</center>
+
+Note that half-edges are only useful for representing orientable manifolds with boundary.  As such, half edge structures couldn't be used to represent a Moebius strip, nor could they be used for meshes where two polygons share a vertex without sharing an edge.  Furthermore, by convention, polygon edges are wound in counter-clockwise order.  We also allow each half-edge to point to an attribute structure for the face that it bounds.  In the case of a Delaunay triangle, that face attribute would be the circumscribing circle's center; edges on the boundary of a mesh have no face attribute (they are stored as `Option[F]` where `F` is the type of the face attribute).

util/src/main/scala/geotrellis/util/Constants.scala

@@ -0,0 +1,17 @@
+package geotrellis.util
+
+object Constants {
+  /**
+   * Defines machine epsilon for double-precision floating point ops.  This is,
+   * roughly speaking, the minimum distance between distinct floating point numbers.
+   * Double values closer than DOUBLE_EPSILON should be considered identical.
+   */
+  val DOUBLE_EPSILON = 1.11e-16
+
+  /**
+   * Defines machine epsilon for single-precision floating point ops.  This is,
+   * roughly speaking, the minimum distance between distinct floating point numbers.
+   * Float values closer than FLOAT_EPSILON should be considered identical.
+   */
+  val FLOAT_EPSILON = 5.96e-8
+}

vector-test/src/test/scala/spec/geotrellis/vector/voronoi/DelaunaySpec.scala

@@ -0,0 +1,114 @@
+package geotrellis.vector.voronoi
+
+import geotrellis.vector._
+import scala.util.Random
+import scala.math.pow
+import org.apache.commons.math3.linear._
+
+import geotrellis.raster._
+import geotrellis.raster.rasterize._
+import geotrellis.raster.render._
+
+import org.scalatest.{FunSpec, Matchers}
+
+class DelaunaySpec extends FunSpec with Matchers {
+
+  val numpts = 2000
+
+  def randomPoint(extent: Extent) : Point = {
+    def randInRange (low : Double, high : Double) : Double = {
+      val x = Random.nextDouble
+      low * (1-x) + high * x
+    }
+    Point(randInRange(extent.xmin, extent.xmax), randInRange(extent.ymin, extent.ymax))
+  }
+
+  def det3 (a11: Double, a12: Double, a13: Double,
+            a21: Double, a22: Double, a23: Double,
+            a31: Double, a32: Double, a33: Double): Double = {
+    val m = MatrixUtils.createRealMatrix(Array(Array(a11, a12, a13),
+                                               Array(a21, a22, a23),
+                                               Array(a31, a32, a33)))
+    (new LUDecomposition(m)).getDeterminant
+  }
+
+  def localInCircle(abc: (Int, Int, Int), di: Int)(implicit trans: Int => Point): Boolean = {
+    val (ai,bi,ci) = abc
+    val a = trans(ai)
+    val b = trans(bi)
+    val c = trans(ci)
+    val d = trans(di)
+    det3(a.x - d.x, a.y - d.y, pow(a.x - d.x, 2) + pow(a.y - d.y, 2),
+         b.x - d.x, b.y - d.y, pow(b.x - d.x, 2) + pow(b.y - d.y, 2),
+         c.x - d.x, c.y - d.y, pow(c.x - d.x, 2) + pow(c.y - d.y, 2)) > 1e-6
+  }
+
+  def rasterizePoly(poly: HalfEdge[Int, Point], tile: MutableArrayTile, re: RasterExtent, erring: Boolean)(implicit trans: Int => Point) = {
+    var pts: List[Point] = Nil
+    var e = poly
+    do {
+      pts = trans(e.vert) :: pts
+      e = e.next
+    } while (e != poly)
+    pts = trans(e.vert) :: pts
+    if (erring) {
+      val p = Polygon(pts)
+      Rasterizer.foreachCellByPolygon(p, re){ (c,r) => tile.set(c, r, 2) }
+    } else {
+      val l = Line(pts)
+      Rasterizer.foreachCellByLineString(l, re){ (c,r) => tile.set(c, r, 1) }
+    }
+  }
+
+  def rasterizeDT(dt: Delaunay)(implicit trans: Int => Point): Unit = {
+    val tile = IntArrayTile.fill(255, 960, 960)
+    val re = RasterExtent(Extent(0,0,1,1),960,960)
+    dt.triangles.foreach{ case ((ai,bi,ci), triEdge) =>
+      val otherPts = (0 until numpts).filter{ i: Int => i != ai && i != bi && i != ci }
+      rasterizePoly(triEdge, tile, re, !otherPts.forall{ i => !localInCircle((ai, bi, ci), i) })
+    }
+    val cm = ColorMap(scala.collection.immutable.Map(1 -> 0x000000ff, 2 -> 0xff0000ff, 255 -> 0xffffffff))
+    tile.renderPng(cm).write("delaunay.png")
+  }
+
+  describe("Delaunay Triangulation") {
+    it("should have a convex boundary") {
+      val range = 0 until numpts
+      val pts = (for (i <- range) yield randomPoint(Extent(0, 0, 1, 1))).toArray
+      implicit val trans = { i: Int => pts(i) }
+      val dt = pts.toList.delaunayTriangulation()
+
+      def boundingEdgeIsConvex(e: HalfEdge[Int, Point]) = {
+        Predicates.isRightOf(e, e.next.vert)
+      }
+      var isConvex = true
+      var e = dt.boundary
+      do {
+        isConvex = isConvex && boundingEdgeIsConvex(e)
+        e = e.next
+      } while (e != dt.boundary)
+
+      isConvex should be (true)
+    }
+
+    it("should preserve Delaunay property") {
+      // Delaunay property: no element of the triangulation should have a 
+      // circumscribing circle that contains another point of the triangulation
+      val range = 0 until numpts
+      val pts = (for (i <- range) yield randomPoint(Extent(0, 0, 1, 1))).toArray
+      val dt = pts.toList.delaunayTriangulation()
+      implicit val trans = { i: Int => pts(i) }
+
+      // NOTE: In the event of failure, the following line will draw the triangulation
+      // to delaunay.png in the working directory, indicating which triangle did not
+      // exhibit the Delaunay property
+      //rasterizeDT(dt)
+
+      (dt.triangles.forall{ case ((ai,bi,ci),_) =>
+        val otherPts = (0 until numpts).filter{ i: Int => i != ai && i != bi && i != ci }
+        otherPts.forall{ i => ! localInCircle((ai,bi,ci),i) }
+      }) should be (true)
+    }
+  }
+
+}

vector-test/src/test/scala/spec/geotrellis/vector/voronoi/VoronoiSpec.scala

@@ -0,0 +1,50 @@
+package geotrellis.vector.voronoi
+
+import geotrellis.vector._
+import scala.util.Random
+import scala.math.pow
+import org.apache.commons.math3.linear._
+
+import geotrellis.raster._
+import geotrellis.raster.rasterize._
+import geotrellis.raster.render._
+
+import org.scalatest.{FunSpec, Matchers}
+
+class VoronoiSpec extends FunSpec with Matchers {
+
+  def rasterizePoly(poly: Polygon, tile: MutableArrayTile, re: RasterExtent, erring: Boolean)(implicit trans: Int => Point) = {
+    if (erring) {
+      Rasterizer.foreachCellByPolygon(poly, re){ (c,r) => tile.set(c, r, 2) }
+    } else {
+      val l = poly.boundary.toGeometry.get
+      Rasterizer.foreachCellByGeometry(l, re){ (c,r) => tile.set(c, r, 1) }
+    }
+  }
+
+  def rasterizeVoronoi(voronoi: Voronoi)(implicit trans: Int => Point): Unit = {
+    val tile = IntArrayTile.fill(255, 325, 600)
+    val re = RasterExtent(voronoi.extent,325,600)
+    voronoi.voronoiCells.foreach{ poly =>
+      rasterizePoly(poly, tile, re, !(poly.isValid && voronoi.extent.covers(poly)))
+    }
+    val cm = ColorMap(scala.collection.immutable.Map(1 -> 0x000000ff, 2 -> 0xff0000ff, 255 -> 0xffffffff))
+    tile.renderPng(cm).write("voronoi.png")
+  }
+
+  describe("Voronoi diagram") {
+    it("should have valid polygons entirely covered by the extent") {
+      val extent = Extent(-2.25, -3, 1, 3)
+      val pts = Array(Point(0,-2), Point(0,0), Point(0,1), Point(-0.5,2), Point(0.5,2))
+      implicit val trans = { i: Int => pts(i) }
+      val voronoi = pts.voronoiDiagram(extent)
+
+      def validCoveredPolygon(poly: Polygon) = {
+        poly.isValid && extent.covers(poly)
+      }
+
+      voronoi.voronoiCells.forall (validCoveredPolygon(_)) should be (true)
+      //rasterizeVoronoi(voronoi)
+    }
+  }
+}

vector/src/main/scala/geotrellis/vector/package.scala

@@ -24,7 +24,8 @@ import scala.collection.mutable
 import scala.collection.JavaConversions._
 
 package object vector extends SeqMethods
-    with reproject.Implicits {
+    with reproject.Implicits 
+    with voronoi.Implicits {
 
   type PointFeature[D] = Feature[Point, D]
   type LineFeature[D] = Feature[Line, D]