From c215dafd780e03d3370803c0bcbc055ab3792097 Mon Sep 17 00:00:00 2001 From: James McClain Date: Thu, 17 Nov 2016 12:04:44 -0500 Subject: [PATCH 1/5] Format, Remove Excess Whitespace --- .../vector/voronoi/DelaunaySpec.scala | 33 +++++++++-------- .../scala/geotrellis/vector/package.scala | 2 +- .../geotrellis/vector/voronoi/Delaunay.scala | 35 +++++++++---------- 3 files changed, 37 insertions(+), 33 deletions(-) diff --git a/vector-test/src/test/scala/spec/geotrellis/vector/voronoi/DelaunaySpec.scala b/vector-test/src/test/scala/spec/geotrellis/vector/voronoi/DelaunaySpec.scala index 99c0fe6596..87a1225fce 100644 --- a/vector-test/src/test/scala/spec/geotrellis/vector/voronoi/DelaunaySpec.scala +++ b/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) } - + } } diff --git a/vector/src/main/scala/geotrellis/vector/package.scala b/vector/src/main/scala/geotrellis/vector/package.scala index d8bc94a77b..aad91d8514 100644 --- a/vector/src/main/scala/geotrellis/vector/package.scala +++ b/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] diff --git a/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala b/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala index 414f9d78eb..115981ca39 100644 --- a/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala +++ b/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala @@ -39,7 +39,6 @@ object Predicates { 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 ] @@ -54,16 +53,16 @@ object Predicates { (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) = + def isRightOf[T](e: HalfEdge[Int,T], p: Point)(implicit trans: Int => Point) = isCCW(p, trans(e.vert), trans(e.src)) - def isRightOf[T](e: HalfEdge[Int,T], p: Int)(implicit trans: Int => Point) = + def isRightOf[T](e: HalfEdge[Int,T], p: Int)(implicit trans: Int => Point) = isCCW(trans(p), trans(e.vert), trans(e.src)) - def isLeftOf[T](e: HalfEdge[Int,T], p: Point)(implicit trans: Int => Point) = + def isLeftOf[T](e: HalfEdge[Int,T], p: Point)(implicit trans: Int => Point) = isCCW(p, trans(e.src), trans(e.vert)) - def isLeftOf[T](e: HalfEdge[Int,T], p: Int)(implicit trans: Int => Point) = + 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 = { @@ -85,14 +84,14 @@ object Predicates { } 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, + 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, + 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, + 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) } @@ -101,14 +100,14 @@ object Predicates { 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, + 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, + 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, + 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) } @@ -145,6 +144,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 -} + } } From 5796f8dfcb55e198446f701243426bf800a18b70 Mon Sep 17 00:00:00 2001 From: James McClain Date: Thu, 17 Nov 2016 12:42:30 -0500 Subject: [PATCH 2/5] Remove Delaunay Predicates --- .../geotrellis/vector/voronoi/Delaunay.scala | 111 +----------------- 1 file changed, 5 insertions(+), 106 deletions(-) diff --git a/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala b/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala index 115981ca39..dce9eb4091 100644 --- a/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala +++ b/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala @@ -16,120 +16,19 @@ 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)) - - def isRightOf[T](e: HalfEdge[Int,T], p: Int)(implicit trans: Int => Point) = - isCCW(trans(p), trans(e.vert), trans(e.src)) - - def isLeftOf[T](e: HalfEdge[Int,T], p: Point)(implicit trans: Int => Point) = - isCCW(p, trans(e.src), trans(e.vert)) - - 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 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 From 584b30948840777a397ad9b3e56bfb2a6e47aa0f Mon Sep 17 00:00:00 2001 From: James McClain Date: Thu, 17 Nov 2016 12:41:35 -0500 Subject: [PATCH 3/5] Add Additional Delaunay Constructors --- .../scala/geotrellis/vector/voronoi/Delaunay.scala | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala b/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala index dce9eb4091..b9d5fcdd9c 100644 --- a/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala +++ b/vector/src/main/scala/geotrellis/vector/voronoi/Delaunay.scala @@ -24,6 +24,19 @@ import com.vividsolutions.jts.triangulate.DelaunayTriangulationBuilder import spire.syntax.cfor._ +object Delaunay { + + def apply(xs: Array[Double], ys: Array[Double]) = + new Delaunay(xs.zip(ys).map({ case (x, y) => Point(x,y) })) + + def apply(coords: Array[Coordinate]) = + new Delaunay(coords.map({ coord => Point(coord.x, coord.y) })) + + def apply(coords: Array[(Double, Double)]) = + new Delaunay(coords.map({ case (x, y) => Point(x, y) })) + +} + /** * A class for triangulating a set of points to satisfy the delaunay * property. Each resulting triangle's circumscribing circle will From 32d8343236f29f2dd93f93b17b9a7a34aa1a2348 Mon Sep 17 00:00:00 2001 From: James McClain Date: Tue, 22 Nov 2016 14:11:49 -0500 Subject: [PATCH 4/5] Conforming Delaunay Triangulation Code --- .../vector/voronoi/ConformingDelaunay.scala | 80 +++++++++++++++++++ .../geotrellis/vector/voronoi/package.scala | 21 +++++ 2 files changed, 101 insertions(+) create mode 100644 vector/src/main/scala/geotrellis/vector/voronoi/ConformingDelaunay.scala create mode 100644 vector/src/main/scala/geotrellis/vector/voronoi/package.scala diff --git a/vector/src/main/scala/geotrellis/vector/voronoi/ConformingDelaunay.scala b/vector/src/main/scala/geotrellis/vector/voronoi/ConformingDelaunay.scala new file mode 100644 index 0000000000..8a0cce9068 --- /dev/null +++ b/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 + } + +} diff --git a/vector/src/main/scala/geotrellis/vector/voronoi/package.scala b/vector/src/main/scala/geotrellis/vector/voronoi/package.scala new file mode 100644 index 0000000000..c5a795e32a --- /dev/null +++ b/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) +} From f147614524f72239a561a958dfe8f79a40c48de5 Mon Sep 17 00:00:00 2001 From: James McClain Date: Tue, 22 Nov 2016 16:29:47 -0500 Subject: [PATCH 5/5] Conforming Delaunay Unit Tests --- .../voronoi/ConformingDelaunaySpec.scala | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 vector-test/src/test/scala/spec/geotrellis/vector/voronoi/ConformingDelaunaySpec.scala diff --git a/vector-test/src/test/scala/spec/geotrellis/vector/voronoi/ConformingDelaunaySpec.scala b/vector-test/src/test/scala/spec/geotrellis/vector/voronoi/ConformingDelaunaySpec.scala new file mode 100644 index 0000000000..2da2e35311 --- /dev/null +++ b/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) + }) + + } + + } + +}