Conforming Delaunay

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

@@ -0,0 +1,44 @@
+package geotrellis.vector.voronoi
+
+import geotrellis.vector._
+
+import com.vividsolutions.jts.{ geom => jts }
+import com.vividsolutions.jts.triangulate.ConformingDelaunayTriangulationBuilder
+import org.scalatest.{FunSpec, Matchers}
+
+
+class ConformingDelaunaySpec extends FunSpec with Matchers {
+
+  val factory = GeomFactory.factory
+  val points = Array(Point(-0.001,-0.001), Point(1,0), Point(2,0), Point(2.001,1), Point(2.001,2.001), Point(1,2), Point(1,1), Point(0,1))
+  val emptyCollection = new jts.GeometryCollection(Array[jts.Geometry](), factory)
+  val polygon = Polygon(points ++ points.take(1))
+
+  describe("Conforming Delaunay Triangulation") {
+
+    it("should be the same as standard Delaunay when no constraints are given") {
+      val conformingDelaunay = ConformingDelaunay(points, emptyCollection)
+      val delaunay = Delaunay(points)
+      conformingDelaunay.triangles.toSet should be (delaunay.triangles.toSet)
+    }
+
+    it("should produce same results with redundant constraints as without any") {
+      val conformingDelaunay = ConformingDelaunay(points, List(polygon))
+      val delaunay = Delaunay(points)
+      conformingDelaunay.triangles.toSet should be (delaunay.triangles.toSet)
+    }
+
+    it("should only produce Steiner points on constraints") {
+      val line = Line(points(1), points(4))
+      val conformingDelaunay = ConformingDelaunay(points, List(line))
+      val steiners = conformingDelaunay.steinerPoints
+
+      steiners.foreach({ point =>
+        line.distance(point) should be < (1e-16)
+      })
+
+    }
+
+  }
+
+}

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

@@ -120,26 +120,31 @@ class DelaunaySpec extends FunSpec with Matchers {
   describe("Delaunay Triangulation") {
 
     ignore ("should have a convex boundary") {
-      /* A test to ensure that the boundary of the triangulation was convex
-       * was once included here, but JTS, for reasons relating to numerical 
-       * robustness, does not produce a triangulation with a guaranteed convex
-       * boundary.  This note is here as a suggestion to future developers to 
-       * include such a test.
-       */
+      /**
+        * A test to ensure that the boundary of the triangulation was
+        * convex was once included here, but JTS, for reasons relating
+        * to numerical robustness, does not produce a triangulation
+        * with a guaranteed convex boundary.  This note is here as a
+        * suggestion to future developers to include such a test.
+        */
     }
 
     it("should preserve Delaunay property") {
-      // Delaunay property: no element of the triangulation should have a 
-      // circumscribing circle that contains another point of the triangulation
+      /**
+        * 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.delaunayTriangulation()
 
-      // 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)
-
+      /**
+        * 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)
+        */
       preservesDelaunay(dt) should be (true)
     }
 
@@ -152,6 +157,6 @@ class DelaunaySpec extends FunSpec with Matchers {
 
       preservesDelaunay(dt) should be (true)
     }
-    
+
   }
 }

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

@@ -24,7 +24,7 @@ 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]

vector/src/main/scala/geotrellis/vector/voronoi/ConformingDelaunay.scala

@@ -0,0 +1,80 @@
+/*
+ * 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.vector.voronoi
+
+import geotrellis.vector._
+
+import com.vividsolutions.jts.{ geom => jts }
+import com.vividsolutions.jts.triangulate.ConformingDelaunayTriangulationBuilder
+
+import spire.syntax.cfor._
+
+
+object ConformingDelaunay {
+
+  def apply(xs: Array[Double], ys: Array[Double], constraints: jts.GeometryCollection) = {
+    new ConformingDelaunay(
+      xs.zip(ys).map({ case (x, y) => Point(x,y) }),
+      constraints
+    )
+  }
+
+  def apply(coords: Array[jts.Coordinate], constraints: jts.GeometryCollection) = {
+    new ConformingDelaunay(
+      coords.map({ coord => Point(coord.x, coord.y) }),
+      constraints
+    )
+  }
+
+  def apply(coords: Array[(Double, Double)], constraints: jts.GeometryCollection) = {
+    new ConformingDelaunay(
+      coords.map({ case (x, y) => Point(x, y) }),
+      constraints
+    )
+  }
+
+}
+
+/**
+  * A class for triangulating a set of points to satisfy the Delaunay
+  * property (subject to the given linear constraints).  Each
+  * resulting triangle's circumscribing circle will contain no other
+  * points of the input set.
+  */
+case class ConformingDelaunay(verts: Array[Point], constraints: jts.GeometryCollection) {
+
+  private[voronoi] val gf = new jts.GeometryFactory
+  private val sites = new jts.MultiPoint(verts.map(_.jtsGeom), gf)
+  private val builder = new ConformingDelaunayTriangulationBuilder
+  builder.setSites(sites) ; builder.setConstraints(constraints)
+  private[voronoi] val subd = builder.getSubdivision
+
+  val triangles: Seq[Polygon] = {
+    val tris = subd.getTriangles(gf)
+    val len = tris.getNumGeometries
+    val arr = Array.ofDim[Polygon](len)
+    cfor(0)(_ < len, _ + 1) { i => arr(i) = Polygon(tris.getGeometryN(i).asInstanceOf[jts.Polygon]) }
+    arr
+  }
+
+  val steinerPoints: Seq[Point] = {
+    val allPoints: Set[Point] = triangles.flatMap({ triangle => triangle.vertices }).toSet
+    val givenPoints: Set[Point] = verts.toSet
+    (allPoints &~ givenPoints).toList
+  }
+
+}

vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala

@@ -16,121 +16,32 @@
 
 package geotrellis.vector.voronoi
 
-import geotrellis.util.Constants.{FLOAT_EPSILON => EPSILON}
-import geotrellis.vector.{Point, Polygon}
+import geotrellis.vector._
 
 import com.vividsolutions.jts.geom.{Coordinate, GeometryFactory, MultiPoint, Polygon => JTSPolygon}
 import com.vividsolutions.jts.triangulate.DelaunayTriangulationBuilder
-import com.vividsolutions.jts.triangulate.quadedge.{QuadEdge}
-import org.apache.commons.math3.linear._
 
-import scala.annotation.tailrec
-import scala.collection.JavaConversions._
-import scala.collection.mutable.{Map, Set}
-import scala.math.pow
 import spire.syntax.cfor._
 
-object Predicates {
-  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 isCCW(a: Point, b: Point, c: Point): Boolean = {
-    // det [ a.x-c.x  a.y-c.y ]
-    //     [ b.x-c.x  b.y-c.y ] > 0
-    (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x) > EPSILON
-  }
-
-  def isCCW(ai: Int, bi: Int, ci: Int)(implicit trans: Int => Point): Boolean = {
-    val a = trans(ai)
-    val b = trans(bi)
-    val c = trans(ci)
-    (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x) > EPSILON
-  }
 
-  def isRightOf[T](e: HalfEdge[Int,T], p: Point)(implicit trans: Int => Point) = 
-    isCCW(p, trans(e.vert), trans(e.src))
+object Delaunay {
 
-  def isRightOf[T](e: HalfEdge[Int,T], p: Int)(implicit trans: Int => Point) = 
-    isCCW(trans(p), trans(e.vert), trans(e.src))
+  def apply(xs: Array[Double], ys: Array[Double]) =
+    new Delaunay(xs.zip(ys).map({ case (x, y) => Point(x,y) }))
 
-  def isLeftOf[T](e: HalfEdge[Int,T], p: Point)(implicit trans: Int => Point) = 
-    isCCW(p, trans(e.src), trans(e.vert))
+  def apply(coords: Array[Coordinate]) =
+    new Delaunay(coords.map({ coord => Point(coord.x, coord.y) }))
 
-  def isLeftOf[T](e: HalfEdge[Int,T], p: Int)(implicit trans: Int => Point) = 
-    isCCW(trans(p), trans(e.src), trans(e.vert))
-
-  def inCircle(abc: (Point, Point, Point), d: Point): Boolean = {
-    val (a,b,c) = abc
-    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)) > EPSILON
-  }
-
-  def inCircle(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)) > EPSILON
-  }
+  def apply(coords: Array[(Double, Double)]) =
+    new Delaunay(coords.map({ case (x, y) => Point(x, y) }))
 
-  def circleCenter(a: Point, b: Point, c: Point): Point = {
-    val d = 2.0 * det3(a.x, a.y, 1.0, 
-                       b.x, b.y, 1.0, 
-                       c.x, c.y, 1.0)
-    val h = det3(a.x * a.x + a.y * a.y, a.y, 1.0, 
-                 b.x * b.x + b.y * b.y, b.y, 1.0, 
-                 c.x * c.x + c.y * c.y, c.y, 1.0) / d
-    val k = det3(a.x, a.x * a.x + a.y * a.y, 1.0, 
-                 b.x, b.x * b.x + b.y * b.y, 1.0, 
-                 c.x, c.x * c.x + c.y * c.y, 1.0) / d
-    Point(h,k)
-  }
-
-  def circleCenter(ai: Int, bi: Int, ci: Int)(implicit trans: Int => Point): Point = {
-    val a = trans(ai)
-    val b = trans(bi)
-    val c = trans(ci)
-    val d = 2.0 * det3(a.x, a.y, 1.0, 
-                       b.x, b.y, 1.0, 
-                       c.x, c.y, 1.0)
-    val h = det3(a.x * a.x + a.y * a.y, a.y, 1.0, 
-                 b.x * b.x + b.y * b.y, b.y, 1.0, 
-                 c.x * c.x + c.y * c.y, c.y, 1.0) / d
-    val k = det3(a.x, a.x * a.x + a.y * a.y, 1.0, 
-                 b.x, b.x * b.x + b.y * b.y, 1.0, 
-                 c.x, c.x * c.x + c.y * c.y, 1.0) / d
-    Point(h,k)
-  }
-
-  def isDelaunayEdge[T](e: HalfEdge[Int,T])(implicit trans: Int => Point): Boolean = {
-    // Predicated on the fact that if an edge is Delaunay, then for a
-    // point, A, to the left of edge (X,Y), and a point, B, to the
-    // right of (X,Y), A may not be in the circle defined by points X,
-    // Y, and B.
-    val a = trans(e.next.vert)
-    val b = trans(e.flip.next.vert)
-    val x = trans(e.flip.vert)
-    val y = trans(e.vert)
-    !inCircle((a, x, y), b)
-  }
 }
 
 /**
- * A class for triangulating a set of points to satisfy the delaunay property.
- * Each resulting triangle's circumscribing circle will contain no other points
- * of the input set.
- */
+  * A class for triangulating a set of points to satisfy the delaunay
+  * property.  Each resulting triangle's circumscribing circle will
+  * contain no other points of the input set.
+  */
 case class Delaunay(verts: Array[Point]) {
 
   private[voronoi] val gf = new GeometryFactory
@@ -145,6 +56,6 @@ case class Delaunay(verts: Array[Point]) {
     val arr = Array.ofDim[Polygon](len)
     cfor(0)(_ < len, _ + 1) { i => arr(i) = Polygon(tris.getGeometryN(i).asInstanceOf[JTSPolygon]) }
     arr
-}
+  }
 
 }

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

@@ -0,0 +1,21 @@
+package geotrellis.vector
+
+import geotrellis.vector._
+
+import com.vividsolutions.jts.{ geom => jts }
+
+
+package object voronoi {
+
+  implicit def linesToCollection(lines: Seq[Line]): jts.GeometryCollection  =
+    new jts.GeometryCollection(lines.map({ line => line.jtsGeom }).toArray, GeomFactory.factory)
+
+  implicit def linesToCollection(lines: Array[Line]): jts.GeometryCollection  =
+    new jts.GeometryCollection(lines.map({ line => line.jtsGeom }), GeomFactory.factory)
+
+  implicit def polygonsToCollection(polygons: Seq[Polygon]): jts.GeometryCollection =
+    new jts.GeometryCollection(polygons.map({ poly => poly.jtsGeom }).toArray, GeomFactory.factory)
+
+  implicit def polygonsToCollection(polygons: Array[Polygon]): jts.GeometryCollection =
+    new jts.GeometryCollection(polygons.map({ poly => poly.jtsGeom }), GeomFactory.factory)
+}