diff --git a/geom.go b/geom.go
index 703f532a..eb8e93ab 100644
--- a/geom.go
+++ b/geom.go
@@ -128,3 +128,79 @@ func GetCoordinates(g Geometry) (pts []Point, err error) {
err = getCoordinates(g, &pts)
return pts, err
}
+
+// extractLines is a helper function for ExtractLines to avoid too many
+// array copies and still provide a convenient interface to the user.
+func extractLines(g Geometry, lines *[]Line) error {
+ switch gg := g.(type) {
+
+ default:
+
+ return ErrUnknownGeometry{g}
+
+ case Pointer:
+
+ return nil
+
+ case MultiPointer:
+
+ return nil
+
+ case LineStringer:
+
+ v := gg.Verticies()
+ for i := 0; i < len(v)-1; i++ {
+ *lines = append(*lines, Line{v[i], v[i+1]})
+ }
+ return nil
+
+ case MultiLineStringer:
+
+ for _, ls := range gg.LineStrings() {
+ if err := extractLines(LineString(ls), lines); err != nil {
+ return err
+ }
+ }
+ return nil
+
+ case Polygoner:
+
+ for _, ls := range gg.LinearRings() {
+ if err := extractLines(LineString(ls), lines); err != nil {
+ return err
+ }
+ }
+ return nil
+
+ case MultiPolygoner:
+
+ for _, p := range gg.Polygons() {
+ if err := extractLines(Polygon(p), lines); err != nil {
+ return err
+ }
+ }
+ return nil
+
+ case Collectioner:
+
+ for _, child := range gg.Geometries() {
+ if err := extractLines(child, lines); err != nil {
+ return err
+ }
+ }
+ return nil
+
+ }
+}
+
+/*
+ExtractLines extracts all linear components from a geometry (line segements).
+If the geometry contains no line segements (e.g. empty geometry or
+point), then an empty array will be returned.
+
+Duplicate lines will not be removed.
+*/
+func ExtractLines(g Geometry) (lines []Line, err error) {
+ err = extractLines(g, &lines)
+ return lines, err
+}
diff --git a/geom_test.go b/geom_test.go
index 3efb2d74..90443c5f 100644
--- a/geom_test.go
+++ b/geom_test.go
@@ -24,6 +24,11 @@ func TestGetCoordinates(t *testing.T) {
}
}
testcases := []tcase{
+ {
+ geom: Extent{},
+ expected: nil,
+ err: ErrUnknownGeometry{Extent{}},
+ },
{
geom: Point{10, 20},
expected: []Point{{10, 20}},
@@ -124,3 +129,127 @@ func TestGetCoordinates(t *testing.T) {
t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
}
}
+
+func TestExtractLines(t *testing.T) {
+
+ type tcase struct {
+ geom Geometry
+ expected []Line
+ err error
+ }
+
+ fn := func(t *testing.T, tc tcase) {
+ r, err := ExtractLines(tc.geom)
+ if err != tc.err {
+ t.Errorf("error, expected %v got %v", tc.err, err)
+ }
+ if !(len(r) == 0 && len(tc.expected) == 0) && !reflect.DeepEqual(r, tc.expected) {
+ t.Errorf("error, expected %v got %v", tc.expected, r)
+ }
+ }
+ testcases := []tcase{
+ {
+ geom: Extent{},
+ expected: nil,
+ err: ErrUnknownGeometry{Extent{}},
+ },
+ {
+ geom: Point{10, 20},
+ expected: []Line{},
+ err: nil,
+ },
+ {
+ geom: MultiPoint{
+ {10, 20},
+ {30, 40},
+ {-10, -5},
+ },
+ expected: []Line{},
+ err: nil,
+ },
+ {
+ geom: LineString{
+ {10, 20},
+ {30, 40},
+ {-10, -5},
+ },
+ expected: []Line{{{10, 20}, {30, 40}}, {{30, 40}, {-10, -5}}},
+ err: nil,
+ },
+ {
+ geom: MultiLineString{
+ {
+ {10, 20},
+ {30, 40},
+ },
+ {
+ {-10, -5},
+ {15, 20},
+ },
+ },
+ expected: []Line{{{10, 20}, {30, 40}}, {{-10, -5}, {15, 20}}},
+ err: nil,
+ },
+ {
+ geom: Polygon{
+ {
+ {10, 20},
+ {30, 40},
+ {-10, -5},
+ },
+ {
+ {1, 2},
+ {3, 4},
+ },
+ },
+ expected: []Line{{{10, 20}, {30, 40}}, {{30, 40}, {-10, -5}}, {{1, 2}, {3, 4}}},
+ err: nil,
+ },
+ {
+ geom: MultiPolygon{
+ {
+ {
+ {10, 20},
+ {30, 40},
+ {-10, -5},
+ },
+ {
+ {1, 2},
+ {3, 4},
+ },
+ },
+ {
+ {
+ {5, 6},
+ {7, 8},
+ {9, 10},
+ },
+ },
+ },
+ expected: []Line{{{10, 20}, {30, 40}}, {{30, 40}, {-10, -5}}, {{1, 2}, {3, 4}}, {{5, 6}, {7, 8}}, {{7, 8}, {9, 10}}},
+ err: nil,
+ },
+ {
+ geom: Collection{
+ Point{10, 20},
+ MultiPoint{
+ {10, 20},
+ {30, 40},
+ {-10, -5},
+ },
+ LineString{
+ {1, 2},
+ {3, 4},
+ {5, 6},
+ },
+ },
+ expected: []Line{{{1, 2}, {3, 4}}, {{3, 4}, {5, 6}}},
+ err: nil,
+ },
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}
diff --git a/planar/triangulate/README.md b/planar/triangulate/README.md
index 1035b2b1..4dfa0b48 100644
--- a/planar/triangulate/README.md
+++ b/planar/triangulate/README.md
@@ -7,9 +7,10 @@ porting are:
* Stay true to JTS's implementation, but make the necessary changes to stay
go-ish and fit nicely within the geom package.
* Only implement what is needed as it is needed.
-* When possible, keep modifications to the functionality segregated into
+* When reasonable, keep modifications to the functionality segregated into
different directories. This will help minimize the cost of porting more of
- JTS's functionality or porting JTS's future changes.
+ JTS's functionality or porting JTS's future changes. When it isn't
+ reasonable, make a specific note of the difference in the comments.
To make porting easier, the original Java code will be kept in the source
files until the specific function has been ported at which time the Java
diff --git a/planar/triangulate/constraineddelaunay/LICENSE b/planar/triangulate/constraineddelaunay/LICENSE
new file mode 100644
index 00000000..cbb04dd1
--- /dev/null
+++ b/planar/triangulate/constraineddelaunay/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2017 go-spatial
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/planar/triangulate/constraineddelaunay/README.md b/planar/triangulate/constraineddelaunay/README.md
new file mode 100644
index 00000000..794219e6
--- /dev/null
+++ b/planar/triangulate/constraineddelaunay/README.md
@@ -0,0 +1,12 @@
+
+
+This constrained delaunay triangulation implementation is based on:
+
+Domiter, Vid. "Constrained Delaunay triangulation using plane subdivision."
+Proceedings of the 8th central European seminar on computer graphics.
+Budmerice. 2004.
+http://old.cescg.org/CESCG-2004/web/Domiter-Vid/CDT.pdf
+
+While this code makes heavy use of the JTS port, it is all new code. As such
+it can fall under the MIT license.
+
diff --git a/planar/triangulate/constraineddelaunay/triangle.go b/planar/triangulate/constraineddelaunay/triangle.go
new file mode 100644
index 00000000..820bc1d9
--- /dev/null
+++ b/planar/triangulate/constraineddelaunay/triangle.go
@@ -0,0 +1,183 @@
+package constraineddelaunay
+
+import (
+ "errors"
+ "fmt"
+
+ "github.com/go-spatial/geom/planar/triangulate/quadedge"
+)
+
+var ErrInvalidVertex = errors.New("invalid vertex")
+var ErrNoMatchingEdgeFound = errors.New("no matching edge found")
+
+/*
+Triangle provides operations on a triangle within a
+quadedge.QuadEdgeSubdivision.
+
+This is outside the quadedge package to avoid making changes to the original
+JTS port.
+*/
+type Triangle struct {
+ // the triangle referenced is to the right of this edge
+ qe *quadedge.QuadEdge
+}
+
+/*
+IntersectsPoint returns true if the vertex intersects the given triangle. This
+includes falling on an edge.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangle) IntersectsPoint(v quadedge.Vertex) bool {
+ e := tri.qe
+
+ for i := 0; i < 3; i++ {
+ lc := v.Classify(e.Orig(), e.Dest())
+ switch lc {
+ // return true if v is on the edge
+ case quadedge.ORIGIN:
+ return true
+ case quadedge.DESTINATION:
+ return true
+ case quadedge.BETWEEN:
+ return true
+ // return false if v is well outside the triangle
+ case quadedge.LEFT:
+ return false
+ case quadedge.BEHIND:
+ return false
+ case quadedge.BEYOND:
+ return false
+ }
+ // go to the next edge of the triangle.
+ e = e.RNext()
+ }
+
+ // if v is to the right of all edges, it is inside the triangle.
+ return true
+}
+
+/*
+opposedTriangle returns the triangle opposite to the vertex v.
+ +
+ /|\
+ / | \
+ / | \
+v1 + a | b +
+ \ | /
+ \ | /
+ \|/
+ +
+
+If this method is called on triangle a with v1 as the vertex, the result will be triangle b.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangle) opposedTriangle(v quadedge.Vertex) (*Triangle, error) {
+ qe := tri.qe
+ for qe.Orig().Equals(v) == false {
+
+ qe = qe.RNext()
+
+ if qe == tri.qe {
+ return nil, ErrInvalidVertex
+ }
+ }
+
+ return &Triangle{qe.RNext().RNext().Sym()}, nil
+}
+
+/*
+opposedVertex returns the vertex opposite to this triangle.
+ +
+ /|\
+ / | \
+ / | \
+v1 + a | b + v2
+ \ | /
+ \ | /
+ \|/
+ +
+
+If this method is called as a.opposedVertex(b), the result will be vertex v2.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangle) opposedVertex(other *Triangle) (quadedge.Vertex, error) {
+ ae, err := tri.sharedEdge(other)
+ if err != nil {
+ return quadedge.Vertex{}, err
+ }
+
+ // using the matching edge in triangle a, find the opposed vertex in b.
+ return ae.Sym().ONext().Dest(), nil
+}
+
+/*
+sharedEdge returns the edge that is shared by both a and b. The edge is
+returned with triangle a on the left.
+
+ + l
+ /|\
+ / | \
+ / | \
+ + a | b +
+ \ | /
+ \ | /
+ \|/
+ + r
+
+If this method is called as a.sharedEdge(b), the result will be edge lr.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangle) sharedEdge(other *Triangle) (*quadedge.QuadEdge, error) {
+ ae := tri.qe
+ be := other.qe
+ foundMatch := false
+
+ // search for the matching edge between both triangles
+ for ai := 0; ai < 3; ai++ {
+ for bi := 0; bi < 3; bi++ {
+ if ae.Orig().Equals(be.Dest()) && ae.Dest().Equals(be.Orig()) {
+ foundMatch = true
+ break
+ }
+ be = be.RNext()
+ }
+
+ if foundMatch {
+ break
+ }
+ ae = ae.RNext()
+ }
+
+ if foundMatch == false {
+ // if there wasn't a matching edge
+ return nil, ErrNoMatchingEdgeFound
+ }
+
+ // return the matching edge in triangle a
+ return ae, nil
+}
+
+/*
+String returns a string representation of triangle.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangle) String() string {
+ str := "["
+ e := tri.qe
+ comma := ""
+ for true {
+ str += comma + fmt.Sprintf("%v", e.Orig())
+ comma = ","
+ e = e.RPrev()
+ if e.Orig().Equals(tri.qe.Orig()) {
+ break
+ }
+ }
+ str = str + "]"
+ return str
+}
diff --git a/planar/triangulate/constraineddelaunay/triangulator.go b/planar/triangulate/constraineddelaunay/triangulator.go
new file mode 100644
index 00000000..9eb18b96
--- /dev/null
+++ b/planar/triangulate/constraineddelaunay/triangulator.go
@@ -0,0 +1,801 @@
+package constraineddelaunay
+
+import (
+ "errors"
+ "fmt"
+ "log"
+ "math"
+
+ "github.com/go-spatial/geom"
+ "github.com/go-spatial/geom/cmp"
+ "github.com/go-spatial/geom/planar/triangulate"
+ "github.com/go-spatial/geom/planar/triangulate/quadedge"
+)
+
+var ErrInvalidPointClassification = errors.New("invalid point classification")
+var ErrLinesDoNotIntersect = errors.New("line segments do not intersect")
+
+// these errors indicate a problem with the algorithm.
+var ErrUnableToUpdateVertexIndex = errors.New("unable to update vertex index")
+var ErrUnexpectedDeadNode = errors.New("unexpected dead node")
+var ErrUnsupportedCoincidentEdges = errors.New("unsupported coincident edges")
+
+/*
+Triangulator provides methods for performing a constrainted delaunay
+triangulation.
+
+Domiter, Vid. "Constrained Delaunay triangulation using plane subdivision."
+Proceedings of the 8th central European seminar on computer graphics.
+Budmerice. 2004.
+http://old.cescg.org/CESCG-2004/web/Domiter-Vid/CDT.pdf
+*/
+type Triangulator struct {
+ builder *triangulate.DelaunayTriangulationBuilder
+ // a map of constraints where the segments have the lesser point first.
+ constraints map[triangulate.Segment]bool
+ subdiv *quadedge.QuadEdgeSubdivision
+ tolerance float64
+ // run validation after many modification operations. This is expensive,
+ // but very useful when debugging.
+ validate bool
+ // maintain an index of vertices to quad edges. Each vertex will point to
+ // one quad edge that has the vertex as an origin. The other quad edges
+ // that point to this vertex can be reached from there.
+ vertexIndex map[quadedge.Vertex]*quadedge.QuadEdge
+}
+
+/*
+appendNonRepeat only appends the provided value if it does not repeat the last
+value that was appended onto the array.
+*/
+func appendNonRepeat(arr []quadedge.Vertex, v quadedge.Vertex) []quadedge.Vertex {
+ if len(arr) == 0 || arr[len(arr)-1].Equals(v) == false {
+ arr = append(arr, v)
+ }
+ return arr
+}
+
+/*
+createSegment creates a segment with vertices a & b, if it doesn't already
+exist. All the vertices must already exist in the triangulator.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) createSegment(s triangulate.Segment) error {
+ qe, err := tri.LocateSegment(s.GetStart(), s.GetEnd())
+ if err != nil && err != quadedge.ErrLocateFailure {
+ return err
+ }
+ if qe != nil {
+ // if the segment already exists
+ return nil
+ }
+
+ ct, err := tri.findIntersectingTriangle(s)
+ if err != nil {
+ return err
+ }
+ from := ct.qe.Sym()
+
+ ct, err = tri.findIntersectingTriangle(triangulate.NewSegment(geom.Line{s.GetEnd(), s.GetStart()}))
+ if err != nil {
+ return err
+ }
+ to := ct.qe.OPrev()
+
+ quadedge.Connect(from, to)
+ // since we aren't adding any vertices we don't need to modify the vertex
+ // index.
+ return nil
+}
+
+/*
+createTriangle creates a triangle with vertices a, b and c. All the vertices
+must already exist in the triangulator. Any existing edges that make up the triangle will not be recreated.
+
+This method makes no effort to ensure the resulting changes are a valid
+triangulation.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) createTriangle(a, b, c quadedge.Vertex) error {
+ if err := tri.createSegment(triangulate.NewSegment(geom.Line{a, b})); err != nil {
+ return err
+ }
+
+ if err := tri.createSegment(triangulate.NewSegment(geom.Line{b, c})); err != nil {
+ return err
+ }
+
+ if err := tri.createSegment(triangulate.NewSegment(geom.Line{c, a})); err != nil {
+ return err
+ }
+
+ return nil
+}
+
+/*
+deleteEdge deletes the specified edge and updates all associated neighbors to
+reflect the removal. The local vertex index is also updated to reflect the
+deletion.
+
+It is invalid to call this method on the last edge that links to a vertex.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) deleteEdge(e *quadedge.QuadEdge) error {
+
+ toRemove := make(map[*quadedge.QuadEdge]bool, 4)
+
+ eSym := e.Sym()
+ eRot := e.Rot()
+ eRotSym := e.Rot().Sym()
+
+ // a set of all the edges that will be removed.
+ toRemove[e] = true
+ toRemove[eSym] = true
+ toRemove[eRot] = true
+ toRemove[eRotSym] = true
+
+ // remove this edge from the vertex index.
+ if err := tri.removeEdgesFromVertexIndex(toRemove, e.Orig()); err != nil {
+ return err
+ }
+ if err := tri.removeEdgesFromVertexIndex(toRemove, e.Dest()); err != nil {
+ return err
+ }
+ quadedge.Splice(e.OPrev(), e)
+ quadedge.Splice(eSym.OPrev(), eSym)
+
+ // TODO: this call is horribly inefficient and should be optimized.
+ tri.subdiv.Delete(e)
+
+ return nil
+}
+
+/*
+findIntersectingTriangle finds the triangle that shares the vertex s.GetStart()
+and intersects at least part of the edge that extends from s.GetStart().
+
+Tolerance is not considered when determining if vertices are the same.
+
+Returns a quadedge that has s.GetStart() as the origin and the right face is
+the desired triangle. If the segment falls on an edge, the triangle to the
+right of the segment is returned.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) findIntersectingTriangle(s triangulate.Segment) (*Triangle, error) {
+
+ qe, err := tri.locateEdgeByVertex(s.GetStart())
+ if err != nil {
+ return nil, err
+ }
+
+ left := qe
+
+ // walk around all the triangles that share qe.Orig()
+ for true {
+ if left.IsLive() == false {
+ return nil, ErrUnexpectedDeadNode
+ }
+ // create the two quad edges around s
+ right := left.OPrev()
+
+ lc := s.GetEnd().Classify(left.Orig(), left.Dest())
+ rc := s.GetEnd().Classify(right.Orig(), right.Dest())
+
+ if (lc == quadedge.RIGHT && rc == quadedge.LEFT) || lc == quadedge.BETWEEN || lc == quadedge.DESTINATION || lc == quadedge.BEYOND {
+ // if s is between the two edges, we found our triangle.
+ return &Triangle{left}, nil
+ } else if lc != quadedge.RIGHT && lc != quadedge.LEFT && rc != quadedge.LEFT && rc != quadedge.RIGHT {
+ // if s falls on lc or rc, then throw an error (for now)
+ // TODO: Handle this case
+ return nil, ErrUnsupportedCoincidentEdges
+ }
+ left = right
+
+ if left == qe {
+ // if we've walked all the way around the vertex.
+ return nil, fmt.Errorf("no intersecting triangle: %v", s)
+ }
+ }
+
+ return nil, fmt.Errorf("no intersecting triangle: %v", s)
+}
+
+/*
+GetEdges gets the edges of the computed triangulation as a MultiLineString.
+
+returns the edges of the triangulation
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) GetEdges() geom.MultiLineString {
+ return tri.builder.GetEdges()
+}
+
+/*
+GetTriangles Gets the faces of the computed triangulation as a
+MultiPolygon.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) GetTriangles() (geom.MultiPolygon, error) {
+ return tri.builder.GetTriangles()
+}
+
+/*
+InsertSegments inserts the line segments in the specified geometry and builds
+a triangulation. The line segments are used as constraints in the
+triangulation. If the geometry is made up solely of points, then no
+constraints will be used.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) InsertSegments(g geom.Geometry) error {
+ err := tri.insertSites(g)
+ if err != nil {
+ return err
+ }
+
+ err = tri.insertConstraints(g)
+ if err != nil {
+ return err
+ }
+
+ return nil
+}
+
+/*
+insertSites inserts all of the vertices found in g into a Delaunay
+triangulation. Other steps will modify the Delaunay Triangulation to create
+the constrained Delaunay triangulation.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) insertSites(g geom.Geometry) error {
+ tri.builder = triangulate.NewDelaunayTriangulationBuilder(tri.tolerance)
+ err := tri.builder.SetSites(g)
+ if err != nil {
+ return err
+ }
+ tri.subdiv = tri.builder.GetSubdivision()
+
+ // Add all the edges to a constant time lookup
+ tri.vertexIndex = make(map[quadedge.Vertex]*quadedge.QuadEdge)
+ edges := tri.subdiv.GetEdges()
+ for i := range edges {
+ e := edges[i]
+ if _, ok := tri.vertexIndex[e.Orig()]; ok == false {
+ tri.vertexIndex[e.Orig()] = e
+ }
+ if _, ok := tri.vertexIndex[e.Dest()]; ok == false {
+ tri.vertexIndex[e.Dest()] = e.Sym()
+ }
+ }
+
+ return nil
+}
+
+/*
+insertConstraints modifies the triangulation by incrementally using the
+line segements in g as constraints in the triangulation. After this step
+the triangulation is no longer a proper Delaunay triangulation, but the
+constraints are guaranteed. Some constraints may need to be split (think
+about the case when two constraints intersect).
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) insertConstraints(g geom.Geometry) error {
+ tri.constraints = make(map[triangulate.Segment]bool)
+
+ lines, err := geom.ExtractLines(g)
+ if err != nil {
+ return fmt.Errorf("error adding constraint: %v", err)
+ }
+ constraints := make(map[triangulate.Segment]bool)
+ for _, l := range lines {
+ // make the line ordering consistent
+ if !cmp.PointLess(l[0], l[1]) {
+ l[0], l[1] = l[1], l[0]
+ }
+
+ seg := triangulate.NewSegment(l)
+ // this maintains the constraints and de-dupes
+ constraints[seg] = true
+ tri.constraints[seg] = true
+ }
+
+ for seg := range constraints {
+ if err := tri.insertEdgeCDT(&seg); err != nil {
+ return fmt.Errorf("error adding constraint: %v", err)
+ }
+ if err = tri.Validate(); err != nil {
+ return err
+ }
+ }
+
+ return nil
+}
+
+/*
+intersection calculates the intersection between two line segments. When the
+rest of geom is ported over from spatial, this can be replaced with a more
+generic call.
+
+The tolerance here only acts by extending the lines by tolerance. E.g. if the
+tolerance is 0.1 and you have two lines {{0, 0}, {1, 0}} and
+{{0, 0.01}, {1, 0.01}} then these will not be marked as intersecting lines.
+
+If tolerance is used to mark two lines as intersecting, you are still
+guaranteed that the intersecting point will fall _on_ one of the lines, not in
+the extended region of the line.
+
+Taken from: https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) intersection(l1, l2 triangulate.Segment) (quadedge.Vertex, error) {
+ p := l1.GetStart()
+ r := l1.GetEnd().Sub(p)
+ q := l2.GetStart()
+ s := l2.GetEnd().Sub(q)
+
+ rs := r.CrossProduct(s)
+
+ if rs == 0 {
+ return quadedge.Vertex{}, ErrLinesDoNotIntersect
+ }
+ t := q.Sub(p).CrossProduct(s.Divide(r.CrossProduct(s)))
+ u := p.Sub(q).CrossProduct(r.Divide(s.CrossProduct(r)))
+
+ // calculate the acceptable range of values for t
+ ttolerance := tri.tolerance / r.Magn()
+ tlow := -ttolerance
+ thigh := 1 + ttolerance
+
+ // calculate the acceptable range of values for u
+ utolerance := tri.tolerance / s.Magn()
+ ulow := -utolerance
+ uhigh := 1 + utolerance
+
+ if t < tlow || t > thigh || u < ulow || u > uhigh {
+ return quadedge.Vertex{}, ErrLinesDoNotIntersect
+ }
+ // if t is just out of range, but within the acceptable tolerance, snap
+ // it back to the beginning/end of the line.
+ t = math.Min(1, math.Max(t, 0))
+
+ return p.Sum(r.Times(t)), nil
+}
+
+/*
+IsConstraint returns true if e is a constrained edge.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) IsConstraint(e *quadedge.QuadEdge) bool {
+
+ _, ok := tri.constraints[triangulate.NewSegment(geom.Line{e.Orig(), e.Dest()})]
+ if ok {
+ return true
+ }
+ _, ok = tri.constraints[triangulate.NewSegment(geom.Line{e.Dest(), e.Orig()})]
+ return ok
+}
+
+/*
+insertEdgeCDT attempts to follow the pseudo code in Domiter.
+
+Procedure InsertEdgeCDT(T:CDT, ab:Edge)
+
+There are some deviations that are also mentioned inline in the comments
+
+ - Some aparrent typos that are resolved to give consistent results
+ - Modifications to work with the planar subdivision representation of
+ a triangulation (QuadEdge)
+ - Modification to support the case when two constrained edges intersect
+ at more than the end points.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) insertEdgeCDT(ab *triangulate.Segment) error {
+
+ qe, err := tri.LocateSegment(ab.GetStart(), ab.GetEnd())
+ if qe != nil && err != nil {
+ return fmt.Errorf("error inserting constraint: %v", err)
+ }
+ if qe != nil {
+ // nothing to do, the edge already exists.
+ return nil
+ }
+
+ // Precondition: a,b in T and ab not in T
+ // Find the triangle t ∈ T that contains a and is cut by ab
+ t, err := tri.findIntersectingTriangle(*ab)
+ if err != nil {
+ return err
+ }
+
+ removalList := make([]*quadedge.QuadEdge, 0)
+
+ // PU:=EmptyList
+ pu := make([]quadedge.Vertex, 0)
+ // PL:=EmptyList
+ pl := make([]quadedge.Vertex, 0)
+ // v:=a
+ v := ab.GetStart()
+ b := ab.GetEnd()
+
+ // While v not in t do -- should this be 'b not in t'!? -JRS
+ for t.IntersectsPoint(b) == false {
+ // tseq:=OpposedTriangle(t,v)
+ tseq, err := t.opposedTriangle(v)
+ if err != nil {
+ return err
+ }
+ // vseq:=OpposesdVertex(tseq,t)
+ vseq, err := tseq.opposedVertex(t)
+ if err != nil {
+ return err
+ }
+
+ shared, err := t.sharedEdge(tseq)
+ if err != nil {
+ return err
+ }
+
+ c := vseq.Classify(ab.GetStart(), ab.GetEnd())
+
+ // should we remove the edge shared between t & tseq?
+ flagEdgeForRemoval := false
+
+ switch {
+
+ case tri.subdiv.IsOnLine(ab.GetLineSegment(), shared.Orig()):
+ // InsertEdgeCDT(T, vseqb)
+ vb := triangulate.NewSegment(geom.Line{shared.Orig(), ab.GetEnd()})
+ tri.insertEdgeCDT(&vb)
+ // a:=vseq -- Should this be b:=vseq!? -JRS
+ b = shared.Orig()
+ *ab = triangulate.NewSegment(geom.Line{ab.GetStart(), b})
+
+ case tri.subdiv.IsOnLine(ab.GetLineSegment(), shared.Dest()):
+ // InsertEdgeCDT(T, vseqb)
+ vb := triangulate.NewSegment(geom.Line{shared.Dest(), ab.GetEnd()})
+ tri.insertEdgeCDT(&vb)
+ // a:=vseq -- Should this be b:=vseq!? -JRS
+ b = shared.Dest()
+ *ab = triangulate.NewSegment(geom.Line{ab.GetStart(), b})
+
+ // if the constrained edge is passing through another constrained edge
+ case tri.IsConstraint(shared):
+ // find the point of intersection
+ iv, err := tri.intersection(*ab, triangulate.NewSegment(geom.Line{shared.Orig(), shared.Dest()}))
+ if err != nil {
+ return err
+ }
+
+ // split the constrained edge we interesect
+ if err := tri.splitEdge(shared, iv); err != nil {
+ return err
+ }
+ tri.deleteEdge(shared)
+ tseq, err = t.opposedTriangle(v)
+ if err != nil {
+ return err
+ }
+
+ // create a new edge for the rest of this segment and recursively
+ // insert the new edge.
+ vb := triangulate.NewSegment(geom.Line{iv, ab.GetEnd()})
+ tri.insertEdgeCDT(&vb)
+
+ // the current insertion will stop at the interesction point
+ b = iv
+ *ab = triangulate.NewSegment(geom.Line{ab.GetStart(), iv})
+
+ // If vseq above the edge ab then
+ case c == quadedge.LEFT:
+ // v:=Vertex shared by t and tseq above ab
+ v = shared.Orig()
+ pu = appendNonRepeat(pu, v)
+ // AddList(PU ,vseq)
+ pu = appendNonRepeat(pu, vseq)
+ flagEdgeForRemoval = true
+
+ // Else If vseq below the edge ab
+ case c == quadedge.RIGHT:
+ // v:=Vertex shared by t and tseq below ab
+ v = shared.Dest()
+ pl = appendNonRepeat(pl, v)
+ // AddList(PL, vseq)
+ pl = appendNonRepeat(pl, vseq)
+ flagEdgeForRemoval = true
+
+ case c == quadedge.DESTINATION:
+ flagEdgeForRemoval = true
+
+ default:
+ return ErrInvalidPointClassification
+ }
+
+ if flagEdgeForRemoval {
+ // "Remove t from T" -- We are just removing the edge intersected
+ // by ab, which in effect removes the triangle.
+ removalList = append(removalList, shared)
+ }
+
+ t = tseq
+ }
+ // EndWhile
+
+ // remove the previously marked edges
+ for i := range removalList {
+ tri.deleteEdge(removalList[i])
+ }
+
+ // TriangulatePseudoPolygon(PU,ab,T)
+ if err := tri.triangulatePseudoPolygon(pu, *ab); err != nil {
+ return err
+ }
+ // TriangulatePseudoPolygon(PL,ab,T)
+ if err := tri.triangulatePseudoPolygon(pl, *ab); err != nil {
+ return err
+ }
+
+ if err := tri.Validate(); err != nil {
+ return err
+ }
+
+ // Add edge ab to T
+ if err := tri.createSegment(*ab); err != nil {
+ return err
+ }
+
+ return nil
+}
+
+/*
+locateEdgeByVertex finds a quad edge that has this vertex as Orig(). This will
+not be a unique edge.
+
+This is looking for an exact match and tolerance will not be considered.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) locateEdgeByVertex(v quadedge.Vertex) (*quadedge.QuadEdge, error) {
+ qe := tri.vertexIndex[v]
+
+ if qe == nil {
+ return nil, quadedge.ErrLocateFailure
+ }
+ return qe, nil
+}
+
+/*
+locateEdgeByVertex finds a quad edge that has this vertex as Orig(). This will
+not be a unique edge.
+
+This is looking for an exact match and tolerance will not be considered.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) LocateSegment(v1 quadedge.Vertex, v2 quadedge.Vertex) (*quadedge.QuadEdge, error) {
+ qe := tri.vertexIndex[v1]
+
+ if qe == nil {
+ return nil, quadedge.ErrLocateFailure
+ }
+
+ start := qe
+ for true {
+ if qe == nil || qe.IsLive() == false {
+ log.Printf("unexpected dead node: %v", qe)
+ return nil, fmt.Errorf("nil or dead qe when locating segment %v %v", v1, v2)
+ }
+ if v2.Equals(qe.Dest()) {
+ return qe, nil
+ }
+
+ qe = qe.ONext()
+ if qe == start {
+ return nil, quadedge.ErrLocateFailure
+ }
+ }
+
+ return qe, nil
+}
+
+/*
+removeConstraintEdge removes any constraints that share the same Orig() and
+Dest() as the edge provided. If there are none, no changes are made.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) removeConstraintEdge(e *quadedge.QuadEdge) {
+ delete(tri.constraints, triangulate.NewSegment(geom.Line{e.Orig(), e.Dest()}))
+ delete(tri.constraints, triangulate.NewSegment(geom.Line{e.Dest(), e.Orig()}))
+}
+
+/*
+removeEdgesFromVertexIndex will remove a set of QuadEdges from the vertex index
+for the specified vertex. If the operation cannot be completed an error will be
+returned and the index will not be modified.
+
+The vertex index maps from a vertex to an arbitrary QuadEdges. This method is
+helpful in modifying the index after an edge has been deleted.
+
+toRemove - a set of QuadEdges that should be removed from the index. These
+QuadEdges don't necessarily have to link to the provided vertex.
+v - The vertex to modify in the index.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) removeEdgesFromVertexIndex(toRemove map[*quadedge.QuadEdge]bool, v quadedge.Vertex) error {
+ ve := tri.vertexIndex[v]
+ if toRemove[ve] {
+ for testEdge := ve.ONext(); ; testEdge = testEdge.ONext() {
+ if testEdge == ve {
+ // if we made it all the way around the vertex without finding
+ // a valid edge to reference from this vertex
+ return ErrUnableToUpdateVertexIndex
+ }
+ if toRemove[testEdge] == false {
+ tri.vertexIndex[v] = testEdge
+ return nil
+ }
+ }
+ }
+ // this should happen if the vertex doesn't need to be updated.
+ return nil
+}
+
+/*
+splitEdge splits the given edge at the vertex v.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) splitEdge(e *quadedge.QuadEdge, v quadedge.Vertex) error {
+ constraint := tri.IsConstraint(e)
+
+ ePrev := e.OPrev()
+ eSym := e.Sym()
+ eSymPrev := eSym.OPrev()
+
+ tri.removeConstraintEdge(e)
+
+ e1 := tri.subdiv.MakeEdge(e.Orig(), v)
+ e2 := tri.subdiv.MakeEdge(e.Dest(), v)
+
+ if _, ok := tri.vertexIndex[v]; ok == false {
+ tri.vertexIndex[v] = e1.Sym()
+ }
+
+ // splice e1 on
+ quadedge.Splice(ePrev, e1)
+ // splice e2 on
+ quadedge.Splice(eSymPrev, e2)
+
+ // splice e1 and e2 together
+ quadedge.Splice(e1.Sym(), e2.Sym())
+
+ if constraint {
+ tri.constraints[triangulate.NewSegment(geom.Line{e1.Orig(), e1.Dest()})] = true
+ tri.constraints[triangulate.NewSegment(geom.Line{e2.Dest(), e2.Orig()})] = true
+ }
+
+ // since we aren't adding any vertices we don't need to modify the vertex
+ // index.
+ return nil
+}
+
+/*
+triangulatePseudoPolygon is taken from the pseudocode TriangulatePseudoPolygon
+from Figure 10 in Domiter.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) triangulatePseudoPolygon(p []quadedge.Vertex, ab triangulate.Segment) error {
+ a := ab.GetStart()
+ b := ab.GetEnd()
+ var c quadedge.Vertex
+ // If P has more than one element then
+ if len(p) > 1 {
+ // c:=First vertex of P
+ c = p[0]
+ ci := 0
+ // For each vertex v in P do
+ for i, v := range p {
+ // If v ∈ CircumCircle (a, b, c) then
+ if quadedge.TrianglePredicate.IsInCircleRobust(a, b, c, v) {
+ c = v
+ ci = i
+ }
+ }
+ // Divide P into PE and PD giving P=PE+c+PD
+ pe := p[0:ci]
+ pd := p[ci+1:]
+ // TriangulatePseudoPolygon(PE, ac, T)
+ if err := tri.triangulatePseudoPolygon(pe, triangulate.NewSegment(geom.Line{a, c})); err != nil {
+ return err
+ }
+ // TriangulatePseudoPolygon(PD, cd, T) (cb instead of cd? -JRS)
+ if err := tri.triangulatePseudoPolygon(pd, triangulate.NewSegment(geom.Line{c, b})); err != nil {
+ return err
+ }
+ } else if len(p) == 1 {
+ c = p[0]
+ }
+
+ // If P is not empty then
+ if len(p) > 0 {
+ // Add triangle with vertices a, b, c into T
+ if err := tri.createTriangle(a, c, b); err != nil {
+ return err
+ }
+ }
+
+ return nil
+}
+
+/*
+validate runs a number of self consistency checks against a triangulation and
+reports the first error.
+
+This is most useful when testing/debugging.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) Validate() error {
+ if tri.validate == false {
+ return nil
+ }
+ err := tri.subdiv.Validate()
+ if err != nil {
+ return err
+ }
+ return tri.validateVertexIndex()
+}
+
+/*
+validateVertexIndex self consistency checks against a triangulation and the
+subdiv and reports the first error.
+
+If tri is nil a panic will occur.
+*/
+func (tri *Triangulator) validateVertexIndex() error {
+ // collect a set of all edges
+ edgeSet := make(map[*quadedge.QuadEdge]bool)
+ vertexSet := make(map[quadedge.Vertex]bool)
+ edges := tri.subdiv.GetEdges()
+ for i := range edges {
+ edgeSet[edges[i]] = true
+ edgeSet[edges[i].Sym()] = true
+ vertexSet[edges[i].Orig()] = true
+ vertexSet[edges[i].Dest()] = true
+ }
+
+ // verify the vertex index points to appropriate edges and vertices
+ for v, e := range tri.vertexIndex {
+ if _, ok := vertexSet[v]; ok == false {
+ return fmt.Errorf("vertex index contains an unexpected vertex: %v", v)
+ }
+ if _, ok := edgeSet[e]; ok == false {
+ return fmt.Errorf("vertex index contains an unexpected edge: %v", e)
+ }
+ if v.Equals(e.Orig()) == false {
+ return fmt.Errorf("vertex index points to an incorrect edge, expected %v got %v", e.Orig(), v)
+ }
+ }
+
+ // verify all vertices are in the vertex index
+ for v, _ := range vertexSet {
+ if _, ok := tri.vertexIndex[v]; ok == false {
+ return fmt.Errorf("vertex index is missing a vertex: %v", v)
+ }
+ }
+
+ return nil
+}
diff --git a/planar/triangulate/constraineddelaunay/triangulator_test.go b/planar/triangulate/constraineddelaunay/triangulator_test.go
new file mode 100644
index 00000000..e7baf6c9
--- /dev/null
+++ b/planar/triangulate/constraineddelaunay/triangulator_test.go
@@ -0,0 +1,358 @@
+package constraineddelaunay
+
+import (
+ "encoding/hex"
+ "log"
+ "strconv"
+ "testing"
+
+ "github.com/go-spatial/geom"
+ "github.com/go-spatial/geom/encoding/wkb"
+ "github.com/go-spatial/geom/encoding/wkt"
+ "github.com/go-spatial/geom/planar/triangulate"
+ "github.com/go-spatial/geom/planar/triangulate/quadedge"
+)
+
+func TestFindIntersectingTriangle(t *testing.T) {
+ type tcase struct {
+ // provided for readability
+ inputWKT string
+ // this can be removed if/when geom has a WKT decoder.
+ // A simple website for performing conversions:
+ // https://rodic.fr/blog/online-conversion-between-geometric-formats/
+ inputWKB string
+ searchFrom geom.Line
+ expectedTriangle string
+ err error
+ }
+
+ fn := func(t *testing.T, tc tcase) {
+ bytes, err := hex.DecodeString(tc.inputWKB)
+ if err != nil {
+ t.Fatalf("error decoding hex string: %v", err)
+ return
+ }
+ g, err := wkb.DecodeBytes(bytes)
+ if err != nil {
+ t.Fatalf("error decoding WKB: %v", err)
+ return
+ }
+
+ uut := new(Triangulator)
+ uut.tolerance = 1e-6
+ // perform self consistency validation while building the
+ // triangulation.
+ uut.validate = true
+ uut.insertSites(g)
+
+ // find the triangle
+ tri, err := uut.findIntersectingTriangle(triangulate.NewSegment(tc.searchFrom))
+ if err != tc.err {
+ t.Fatalf("error, expected %v got %v", tc.err, err)
+ return
+ }
+
+ if tc.err == nil {
+ qeStr := tri.String()
+ if qeStr != tc.expectedTriangle {
+ t.Fatalf("error, expected %v got %v", tc.expectedTriangle, qeStr)
+ }
+ }
+ }
+ testcases := []tcase{
+ {
+ inputWKT: `MULTIPOINT (10 10, 10 20, 20 20, 20 10, 20 0, 10 0, 0 0, 0 10, 0 20)`,
+ inputWKB: `010400000009000000010100000000000000000024400000000000002440010100000000000000000024400000000000003440010100000000000000000034400000000000003440010100000000000000000034400000000000002440010100000000000000000034400000000000000000010100000000000000000024400000000000000000010100000000000000000000000000000000000000010100000000000000000000000000000000002440010100000000000000000000000000000000003440`,
+ searchFrom: geom.Line{{0, 0}, {10, 10}},
+ expectedTriangle: `[[0 0],[0 10],[10 0]]`,
+ },
+ {
+ inputWKT: `MULTIPOINT (10 10, 10 20, 20 20, 20 10, 20 0, 10 0, 0 0, 0 10, 0 20)`,
+ inputWKB: `010400000009000000010100000000000000000024400000000000002440010100000000000000000024400000000000003440010100000000000000000034400000000000003440010100000000000000000034400000000000002440010100000000000000000034400000000000000000010100000000000000000024400000000000000000010100000000000000000000000000000000000000010100000000000000000000000000000000002440010100000000000000000000000000000000003440`,
+ searchFrom: geom.Line{{10, 0}, {0, 20}},
+ expectedTriangle: `[[10 0],[0 10],[10 10]]`,
+ },
+ {
+ inputWKT: `MULTIPOINT (10 10, 10 20, 20 20, 20 10, 20 0, 10 0, 0 0, 0 10, 0 20)`,
+ inputWKB: `010400000009000000010100000000000000000024400000000000002440010100000000000000000024400000000000003440010100000000000000000034400000000000003440010100000000000000000034400000000000002440010100000000000000000034400000000000000000010100000000000000000024400000000000000000010100000000000000000000000000000000000000010100000000000000000000000000000000002440010100000000000000000000000000000000003440`,
+ searchFrom: geom.Line{{10, 10}, {0, 0}},
+ expectedTriangle: `[[10 10],[10 0],[0 10]]`,
+ },
+ {
+ inputWKT: `MULTIPOINT (10 10, 10 20, 20 20, 20 10, 20 0, 10 0, 0 0, 0 10, 0 20)`,
+ inputWKB: `010400000009000000010100000000000000000024400000000000002440010100000000000000000024400000000000003440010100000000000000000034400000000000003440010100000000000000000034400000000000002440010100000000000000000034400000000000000000010100000000000000000024400000000000000000010100000000000000000000000000000000000000010100000000000000000000000000000000002440010100000000000000000000000000000000003440`,
+ searchFrom: geom.Line{{10, 10}, {10, 20}},
+ expectedTriangle: `[[10 10],[10 20],[20 10]]`,
+ },
+ {
+ inputWKT: `MULTIPOINT (10 10, 10 20, 20 20, 20 10, 20 0, 10 0, 0 0, 0 10, 0 20)`,
+ inputWKB: `010400000009000000010100000000000000000024400000000000002440010100000000000000000024400000000000003440010100000000000000000034400000000000003440010100000000000000000034400000000000002440010100000000000000000034400000000000000000010100000000000000000024400000000000000000010100000000000000000000000000000000000000010100000000000000000000000000000000002440010100000000000000000000000000000000003440`,
+ searchFrom: geom.Line{{1000, 1000}, {10000, 20000}},
+ expectedTriangle: ``,
+ err: quadedge.ErrLocateFailure,
+ },
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}
+
+func TestDeleteEdge(t *testing.T) {
+ type tcase struct {
+ // provided for readability
+ inputWKT string
+ // this can be removed if/when geom has a WKT decoder.
+ // A simple website for performing conversions:
+ // https://rodic.fr/blog/online-conversion-between-geometric-formats/
+ inputWKB string
+ deleteMe geom.Line
+ }
+
+ fn := func(t *testing.T, tc tcase) {
+ bytes, err := hex.DecodeString(tc.inputWKB)
+ if err != nil {
+ t.Fatalf("error decoding hex string, expected nil got %v", err)
+ return
+ }
+ g, err := wkb.DecodeBytes(bytes)
+ if err != nil {
+ t.Fatalf("error decoding WKB, expected nil got %v", err)
+ return
+ }
+
+ uut := new(Triangulator)
+ uut.tolerance = 1e-6
+ // perform self consistency validation while building the
+ // triangulation.
+ uut.validate = true
+ uut.InsertSegments(g)
+ e, err := uut.LocateSegment(quadedge.Vertex(tc.deleteMe[0]), quadedge.Vertex(tc.deleteMe[1]))
+ if err != nil {
+ t.Fatalf("error locating segment, expected nil got %v", err)
+ return
+ }
+
+ err = uut.Validate()
+ if err != nil {
+ t.Errorf("error validating triangulation, expected nil got %v", err)
+ return
+ }
+
+ if err = uut.deleteEdge(e); err != nil {
+ t.Errorf("error deleting edge, expected nil got %v", err)
+ }
+ err = uut.Validate()
+ if err != nil {
+ t.Errorf("error validating triangulation after delete, expected nil got %v", err)
+ return
+ }
+
+ // this edge shouldn't exist anymore.
+ _, err = uut.LocateSegment(quadedge.Vertex(tc.deleteMe[0]), quadedge.Vertex(tc.deleteMe[1]))
+ if err == nil {
+ t.Fatalf("error locating segment, expected %v got nil", quadedge.ErrLocateFailure)
+ return
+ }
+
+ }
+ testcases := []tcase{
+ {
+ inputWKT: `MULTIPOINT (10 10, 10 20, 20 20, 20 10, 20 0, 10 0, 0 0, 0 10, 0 20)`,
+ inputWKB: `010400000009000000010100000000000000000024400000000000002440010100000000000000000024400000000000003440010100000000000000000034400000000000003440010100000000000000000034400000000000002440010100000000000000000034400000000000000000010100000000000000000024400000000000000000010100000000000000000000000000000000000000010100000000000000000000000000000000002440010100000000000000000000000000000000003440`,
+ deleteMe: geom.Line{{0, 10}, {10, 0}},
+ },
+ {
+ inputWKT: `MULTIPOINT (10 10, 10 20, 20 20, 20 10, 20 0, 10 0, 0 0, 0 10, 0 20)`,
+ inputWKB: `010400000009000000010100000000000000000024400000000000002440010100000000000000000024400000000000003440010100000000000000000034400000000000003440010100000000000000000034400000000000002440010100000000000000000034400000000000000000010100000000000000000024400000000000000000010100000000000000000000000000000000000000010100000000000000000000000000000000002440010100000000000000000000000000000000003440`,
+ deleteMe: geom.Line{{0, 10}, {10, 0}},
+ },
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}
+
+func TestIntersection(t *testing.T) {
+ type tcase struct {
+ l1 triangulate.Segment
+ l2 triangulate.Segment
+ intersection quadedge.Vertex
+ expectedError error
+ }
+
+ fn := func(t *testing.T, tc tcase) {
+ uut := new(Triangulator)
+ uut.tolerance = 1e-2
+ v, err := uut.intersection(tc.l1, tc.l2)
+ if err != tc.expectedError {
+ t.Errorf("error intersecting line segments, expected %v got %v", tc.expectedError, err)
+ return
+ }
+
+ if v.Equals(tc.intersection) == false {
+ t.Errorf("error validating intersection, expected %v got %v", tc.intersection, v)
+ }
+ }
+ testcases := []tcase{
+ {
+ l1: triangulate.NewSegment(geom.Line{{0, 1}, {2, 3}}),
+ l2: triangulate.NewSegment(geom.Line{{1, 1}, {0, 2}}),
+ intersection: quadedge.Vertex{0.5, 1.5},
+ expectedError: nil,
+ },
+ {
+ l1: triangulate.NewSegment(geom.Line{{0, 1}, {2, 4}}),
+ l2: triangulate.NewSegment(geom.Line{{1, 1}, {0, 2}}),
+ intersection: quadedge.Vertex{0.4, 1.6},
+ expectedError: nil,
+ },
+ {
+ l1: triangulate.NewSegment(geom.Line{{0, 1}, {2, 3}}),
+ l2: triangulate.NewSegment(geom.Line{{1, 1}, {2, 2}}),
+ intersection: quadedge.Vertex{0, 0},
+ expectedError: ErrLinesDoNotIntersect,
+ },
+ {
+ l1: triangulate.NewSegment(geom.Line{{3, 5}, {3, 6}}),
+ l2: triangulate.NewSegment(geom.Line{{1, 4.995}, {4, 4.995}}),
+ intersection: quadedge.Vertex{3, 5},
+ expectedError: nil,
+ },
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}
+
+/*
+TestTriangulation test cases test for small constrained triangulations and
+edge cases
+*/
+func TestTriangulation(t *testing.T) {
+ type tcase struct {
+ // provided for readability
+ inputWKT string
+ // this can be removed if/when geom has a WKT decoder.
+ // A simple website for performing conversions:
+ // https://rodic.fr/blog/online-conversion-between-geometric-formats/
+ inputWKB string
+ expectedEdges string
+ expectedTris string
+ }
+
+ // to change the flags on the default logger
+ log.SetFlags(log.LstdFlags | log.Lshortfile)
+
+ fn := func(t *testing.T, tc tcase) {
+ bytes, err := hex.DecodeString(tc.inputWKB)
+ if err != nil {
+ t.Fatalf("error decoding hex string: %v", err)
+ return
+ }
+ g, err := wkb.DecodeBytes(bytes)
+ if err != nil {
+ t.Fatalf("error decoding WKB: %v", err)
+ return
+ }
+
+ uut := new(Triangulator)
+ uut.tolerance = 1e-6
+ err = uut.InsertSegments(g)
+ if err != nil {
+ t.Fatalf("error inserting segments, expected nil got %v", err)
+ }
+
+ edges := uut.GetEdges()
+ edgesWKT, err := wkt.Encode(edges)
+ if err != nil {
+ t.Errorf("error, expected nil got %v", err)
+ return
+ }
+ if edgesWKT != tc.expectedEdges {
+ t.Errorf("error, expected %v got %v", tc.expectedEdges, edgesWKT)
+ return
+ }
+
+ tris, err := uut.GetTriangles()
+ if err != nil {
+ t.Errorf("error, expected nil got %v", err)
+ return
+ }
+ trisWKT, err := wkt.Encode(tris)
+ if err != nil {
+ t.Errorf("error, expected nil got %v", err)
+ return
+ }
+ if trisWKT != tc.expectedTris {
+ t.Errorf("error, expected %v got %v", tc.expectedTris, trisWKT)
+ return
+ }
+ }
+ testcases := []tcase{
+ {
+ // should create a triangulation w/ a vertical line (2 5, 2 -5).
+ // The unconstrained version has a horizontal line
+ inputWKT: `LINESTRING(0 0, 2 5, 2 -5, 5 0)`,
+ inputWKB: `0102000000040000000000000000000000000000000000000000000000000000400000000000001440000000000000004000000000000014c000000000000014400000000000000000`,
+ expectedEdges: `MULTILINESTRING ((2 5,5 0),(0 0,2 5),(0 0,2 -5),(2 -5,5 0),(2 -5,2 5))`,
+ expectedTris: `MULTIPOLYGON (((0 0,2 -5,2 5,0 0)),((2 5,2 -5,5 0,2 5)))`,
+ },
+ {
+ // a horizontal rectangle w/ one diagonal line. The diagonal line
+ // should be maintained and the top/bottom re-triangulated.
+ inputWKT: `MULTILINESTRING((0 0,0 1,1 1.1,2 1,2 0,1 0.1,0 0),(0 1,2 0))`,
+ inputWKB: `010500000002000000010200000007000000000000000000000000000000000000000000000000000000000000000000f03f000000000000f03f9a9999999999f13f0000000000000040000000000000f03f00000000000000400000000000000000000000000000f03f9a9999999999b93f000000000000000000000000000000000102000000020000000000000000000000000000000000f03f00000000000000400000000000000000`,
+ expectedEdges: `MULTILINESTRING ((1 1.1,2 1),(0 1,1 1.1),(0 0,0 1),(0 0,2 0),(2 0,2 1),(1 1.1,2 0),(0 1,2 0),(1 0.1,2 0),(0 1,1 0.1),(0 0,1 0.1))`,
+ expectedTris: `MULTIPOLYGON (((0 1,0 0,1 0.1,0 1)),((0 1,1 0.1,2 0,0 1)),((0 1,2 0,1 1.1,0 1)),((1 1.1,2 0,2 1,1 1.1)),((0 0,2 0,1 0.1,0 0)))`,
+ },
+ {
+ // an egg shape with one horizontal line. The horizontal line
+ // should be maintained and the top/bottom re-triangulated.
+ inputWKT: `MULTILINESTRING((0 0,-0.1 0.5,0 1,0.5 1.2,1 1.3,1.5 1.2,2 1,2.1 0.5,2 0,1.5 -0.2,1 -0.3,0.5 -0.2,0 0),(-0.1 0.5,2.1 0.5))`,
+ inputWKB: `01050000000200000001020000000d000000000000000000000000000000000000009a9999999999b9bf000000000000e03f0000000000000000000000000000f03f000000000000e03f333333333333f33f000000000000f03fcdccccccccccf43f000000000000f83f333333333333f33f0000000000000040000000000000f03fcdcccccccccc0040000000000000e03f00000000000000400000000000000000000000000000f83f9a9999999999c9bf000000000000f03f333333333333d3bf000000000000e03f9a9999999999c9bf000000000000000000000000000000000102000000020000009a9999999999b9bf000000000000e03fcdcccccccccc0040000000000000e03f`,
+ expectedEdges: `MULTILINESTRING ((1.5 1.2,2 1),(1 1.3,1.5 1.2),(0.5 1.2,1 1.3),(0 1,0.5 1.2),(-0.1 0.5,0 1),(-0.1 0.5,0 0),(0 0,0.5 -0.2),(0.5 -0.2,1 -0.3),(1 -0.3,1.5 -0.2),(1.5 -0.2,2 0),(2 0,2.1 0.5),(2 1,2.1 0.5),(1.5 1.2,2.1 0.5),(1 1.3,2.1 0.5),(-0.1 0.5,2.1 0.5),(-0.1 0.5,1 1.3),(-0.1 0.5,0.5 1.2),(1.5 -0.2,2.1 0.5),(-0.1 0.5,1.5 -0.2),(0.5 -0.2,1.5 -0.2),(-0.1 0.5,0.5 -0.2))`,
+ expectedTris: `MULTIPOLYGON (((0 1,-0.1 0.5,0.5 1.2,0 1)),((0.5 1.2,-0.1 0.5,1 1.3,0.5 1.2)),((1 1.3,-0.1 0.5,2.1 0.5,1 1.3)),((1 1.3,2.1 0.5,1.5 1.2,1 1.3)),((1.5 1.2,2.1 0.5,2 1,1.5 1.2)),((1 -0.3,1.5 -0.2,0.5 -0.2,1 -0.3)),((0.5 -0.2,1.5 -0.2,-0.1 0.5,0.5 -0.2)),((0.5 -0.2,-0.1 0.5,0 0,0.5 -0.2)),((-0.1 0.5,1.5 -0.2,2.1 0.5,-0.1 0.5)),((2.1 0.5,1.5 -0.2,2 0,2.1 0.5)))`,
+ },
+ {
+ // a triangle with a line intersecting the top vertex. Where the
+ // line intersects the vertex, the line should be broken into two
+ // pieces and triangulated properly.
+ inputWKT: `MULTILINESTRING((0 0,-0.1 0.5,0 1,0.5 1.2,1 1.3,1.5 1.2,2 1,2.1 0.5,2 0,1.5 -0.2,1 -0.3,0.5 -0.2,0 0),(-0.1 0.5,2.1 0.5))`,
+ inputWKB: `01050000000200000001020000000400000000000000000000000000000000000000000000000000f03f000000000000f03f00000000000000400000000000000000000000000000000000000000000000000102000000020000000000000000000000000000000000f03f0000000000000040000000000000f03f`,
+ expectedEdges: `MULTILINESTRING ((1 1,2 1),(0 1,1 1),(0 0,0 1),(0 0,2 0),(2 0,2 1),(1 1,2 0),(0 0,1 1))`,
+ expectedTris: `MULTIPOLYGON (((0 1,0 0,1 1,0 1)),((1 1,0 0,2 0,1 1)),((1 1,2 0,2 1,1 1)))`,
+ },
+ {
+ // a figure eight with a duplicate constrained line.
+ inputWKT: `MULTIPOLYGON (((0 0,0 1,1 1,1 0,0 0,0 -1,1 -1,1 0,0 0)))`,
+ inputWKB: `01060000000100000001030000000100000009000000000000000000000000000000000000000000000000000000000000000000f03f000000000000f03f000000000000f03f000000000000f03f0000000000000000000000000000000000000000000000000000000000000000000000000000f0bf000000000000f03f000000000000f0bf000000000000f03f000000000000000000000000000000000000000000000000`,
+ expectedEdges: `MULTILINESTRING ((0 1,1 1),(0 0,0 1),(0 -1,0 0),(0 -1,1 -1),(1 -1,1 0),(1 0,1 1),(0 1,1 0),(0 0,1 0),(0 0,1 -1))`,
+ expectedTris: `MULTIPOLYGON (((0 1,0 0,1 0,0 1)),((0 1,1 0,1 1,0 1)),((0 -1,1 -1,0 0,0 -1)),((0 0,1 -1,1 0,0 0)))`,
+ },
+ {
+ // A constraint line that overlaps with another edge
+ inputWKT: `MULTIPOLYGON (((0 0,1 1,2 1,3 0,3 1,0 1,0 0)))`,
+ inputWKB: `0106000000010000000103000000010000000700000000000000000000000000000000000000000000000000f03f000000000000f03f0000000000000040000000000000f03f000000000000084000000000000000000000000000000840000000000000f03f0000000000000000000000000000f03f00000000000000000000000000000000`,
+ expectedEdges: `MULTILINESTRING ((2 1,3 1),(1 1,2 1),(0 1,1 1),(0 0,0 1),(0 0,3 0),(3 0,3 1),(2 1,3 0),(0 0,2 1),(0 0,1 1))`,
+ expectedTris: `MULTIPOLYGON (((0 1,0 0,1 1,0 1)),((1 1,0 0,2 1,1 1)),((2 1,0 0,3 0,2 1)),((2 1,3 0,3 1,2 1)))`,
+ },
+ {
+ // bow-tie
+ inputWKT: `MULTIPOLYGON (((0 0,1 1,1 0,0 1,0 0)))`,
+ inputWKB: `0106000000010000000103000000010000000500000000000000000000000000000000000000000000000000f03f000000000000f03f000000000000f03f00000000000000000000000000000000000000000000f03f00000000000000000000000000000000`,
+ expectedEdges: `MULTILINESTRING ((0 1,1 1),(0 0,0 1),(0 0,1 0),(1 0,1 1),(0.5 0.5,1 0),(0.5 0.5,1 1),(0 1,0.5 0.5),(0 0,0.5 0.5))`,
+ expectedTris: `MULTIPOLYGON (((0 1,0 0,0.5 0.5,0 1)),((0 1,0.5 0.5,1 1,0 1)),((1 1,0.5 0.5,1 0,1 1)),((0 0,1 0,0.5 0.5,0 0)))`,
+ },
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}
diff --git a/planar/triangulate/delaunay_test.go b/planar/triangulate/delaunay_test.go
index a96c8061..af8d69fb 100644
--- a/planar/triangulate/delaunay_test.go
+++ b/planar/triangulate/delaunay_test.go
@@ -30,6 +30,8 @@ func TestDelaunayTriangulation(t *testing.T) {
// provided for readability
inputWKT string
// this can be removed if/when geom has a WKT decoder.
+ // A simple website for performing conversions:
+ // https://rodic.fr/blog/online-conversion-between-geometric-formats/
inputWKB string
expectedEdges string
expectedTris string
@@ -47,11 +49,17 @@ func TestDelaunayTriangulation(t *testing.T) {
return
}
- builder := new(DelaunayTriangulationBuilder)
- builder.tolerance = 1e-6
+ builder := NewDelaunayTriangulationBuilder(1e-6)
builder.SetSites(sites)
+ if builder.create() == false {
+ t.Errorf("error building triangulation, expected true got false")
+ }
+ err = builder.subdiv.Validate()
+ if err != nil {
+ t.Errorf("error, expected nil got %v", err)
+ }
- edges := builder.getEdges()
+ edges := builder.GetEdges()
edgesWKT, err := wkt.Encode(edges)
if err != nil {
t.Errorf("error, expected nil got %v", err)
diff --git a/planar/triangulate/delaunaytriangulationbuilder.go b/planar/triangulate/delaunaytriangulationbuilder.go
index 441e9dc2..d14f8003 100644
--- a/planar/triangulate/delaunaytriangulationbuilder.go
+++ b/planar/triangulate/delaunaytriangulationbuilder.go
@@ -40,6 +40,10 @@ func (xy PointByXY) Less(i, j int) bool { return cmp.XYLessPoint(xy[i], xy[j]) }
func (xy PointByXY) Swap(i, j int) { xy[i], xy[j] = xy[j], xy[i] }
func (xy PointByXY) Len() int { return len(xy) }
+func NewDelaunayTriangulationBuilder(tolerance float64) *DelaunayTriangulationBuilder {
+ return &DelaunayTriangulationBuilder{tolerance: tolerance}
+}
+
/*
extractUniqueCoordinates extracts the unique points from the given Geometry.
@@ -208,7 +212,7 @@ returns the edges of the triangulation
If dtb is nil a panic will occur.
*/
-func (dtb *DelaunayTriangulationBuilder) getEdges() geom.MultiLineString {
+func (dtb *DelaunayTriangulationBuilder) GetEdges() geom.MultiLineString {
if !dtb.create() {
return geom.MultiLineString{}
}
diff --git a/planar/triangulate/delaunaytriangulationbuilder_test.go b/planar/triangulate/delaunaytriangulationbuilder_test.go
index 8a70465e..dfcc4cc7 100644
--- a/planar/triangulate/delaunaytriangulationbuilder_test.go
+++ b/planar/triangulate/delaunaytriangulationbuilder_test.go
@@ -32,6 +32,10 @@ func TestUnique(t *testing.T) {
if reflect.DeepEqual(result, tc.expected) == false {
t.Errorf("error, expected %v got %v", tc.expected, result)
}
+ // This shouldn't exist with no data
+ if uut.GetSubdivision() != nil {
+ t.Errorf("error, expected nil got not nil")
+ }
}
testcases := []tcase{
{
diff --git a/planar/triangulate/incrementaldelaunaytriangulator.go b/planar/triangulate/incrementaldelaunaytriangulator.go
index f4441817..1768a877 100644
--- a/planar/triangulate/incrementaldelaunaytriangulator.go
+++ b/planar/triangulate/incrementaldelaunaytriangulator.go
@@ -100,12 +100,12 @@ func (idt *IncrementalDelaunayTriangulator) InsertSite(v quadedge.Vertex) (*quad
// log.Printf("Made Edge: %v -> %v", base.Orig(), base.Dest());
quadedge.Splice(base, e)
startEdge := base
- done := false
- for !done {
+
+ for true {
base = idt.subdiv.Connect(e, base.Sym())
e = base.OPrev()
if e.LNext() == startEdge {
- done = true
+ break
}
}
diff --git a/planar/triangulate/quadedge/quadedge.go b/planar/triangulate/quadedge/quadedge.go
index 8a90b6ea..239b4207 100644
--- a/planar/triangulate/quadedge/quadedge.go
+++ b/planar/triangulate/quadedge/quadedge.go
@@ -71,6 +71,9 @@ func MakeEdge(o Vertex, d Vertex) *QuadEdge {
base := q0
base.setOrig(o)
base.setDest(d)
+ base.rot.setOrig(o)
+ base.rot.setDest(d)
+
return base
}
@@ -432,7 +435,7 @@ public LineSegment toLineSegment()
*/
/*
-Converts this edge to a WKT two-point LINESTRING indicating
+String Converts this edge to a WKT two-point LINESTRING indicating
the geometry of this edge.
Unlike JTS, if IsLive() is false, a deleted string is returned.
diff --git a/planar/triangulate/quadedge/quadedgesubdivision.go b/planar/triangulate/quadedge/quadedgesubdivision.go
index c2c43e5d..69921429 100644
--- a/planar/triangulate/quadedge/quadedgesubdivision.go
+++ b/planar/triangulate/quadedge/quadedgesubdivision.go
@@ -15,6 +15,7 @@ package quadedge
import (
"errors"
"fmt"
+ "log"
"github.com/go-spatial/geom"
"github.com/go-spatial/geom/planar"
@@ -224,19 +225,32 @@ func (qes *QuadEdgeSubdivision) Delete(e *QuadEdge) {
eRot := e.Rot()
eRotSym := e.Rot().Sym()
+ e.Delete()
+ eSym.Delete()
+ eRot.Delete()
+ eRotSym.Delete()
+
// this is inefficient on an array, but this method should be called
// infrequently
newArray := make([]*QuadEdge, 0, len(qes.quadEdges))
for _, ele := range qes.quadEdges {
- if ele != e && ele != eSym && ele != eRot && ele != eRotSym {
+ if ele.IsLive() {
newArray = append(newArray, ele)
+
+ if ele.next.IsLive() == false {
+ log.Fatalf("a dead edge is still linked: %v", ele)
+ }
}
}
+ qes.quadEdges = newArray
- e.Delete()
- eSym.Delete()
- eRot.Delete()
- eRotSym.Delete()
+ if qes.startingEdge.IsLive() == false {
+ if len(qes.quadEdges) > 0 {
+ qes.startingEdge = qes.quadEdges[0]
+ } else {
+ qes.startingEdge = nil
+ }
+ }
}
/*
@@ -268,21 +282,16 @@ func (qes *QuadEdgeSubdivision) LocateFromEdge(v Vertex, startEdge *QuadEdge) (*
iter++
/*
- So far it has always been the case that failure to locate indicates an
- invalid subdivision. So just fail completely. (An alternative would be
- to perform an exhaustive search for the containing triangle, but this
- would mask errors in the subdivision topology)
+ So far it has always been the case that failure to locate indicates an
+ invalid subdivision. So just fail completely. (An alternative would be
+ to perform an exhaustive search for the containing triangle, but this
+ would mask errors in the subdivision topology)
- This can also happen if two vertices are located very close together,
- since the orientation predicates may experience precision failures.
+ This can also happen if two vertices are located very close together,
+ since the orientation predicates may experience precision failures.
*/
if iter > maxIter {
return nil, ErrLocateFailure
- // String msg = "Locate failed to converge (at edge: " + e + ").
- // Possible causes include invalid Subdivision topology or very close
- // sites";
- // System.err.println(msg);
- // dumpTriangles();
}
if v.Equals(e.Orig()) || v.Equals(e.Dest()) {
@@ -476,7 +485,22 @@ func (qes *QuadEdgeSubdivision) IsOnEdge(e *QuadEdge, p geom.Pointer) bool {
}
/*
-IsVertexOfEdge tests whether a {@link Vertex} is the start or end vertex of a
+IsOnLine Tests whether a point lies on a segment, up to a tolerance
+determined by the subdivision tolerance.
+
+Returns true if the vertex lies on the edge
+
+If qes is nil a panic will occur.
+*/
+func (qes *QuadEdgeSubdivision) IsOnLine(l geom.Line, p geom.Pointer) bool {
+ dist := planar.DistanceToLineSegment(p, geom.Point(l[0]), geom.Point(l[1]))
+
+ // heuristic (hack?)
+ return dist < qes.edgeCoincidenceTolerance
+}
+
+/*
+IsVertexOfEdge tests whether a Vertex is the start or end vertex of a
QuadEdge, up to the subdivision tolerance distance.
Returns true if the vertex is a endpoint of the edge
@@ -626,6 +650,9 @@ func (qes *QuadEdgeSubdivision) GetPrimaryEdges(includeFrame bool) []*QuadEdge {
qes.visitedKey++
var edges []*QuadEdge
+ if qes.startingEdge == nil {
+ return edges
+ }
var stack edgeStack
stack.push(qes.startingEdge)
@@ -690,7 +717,9 @@ func (qes *QuadEdgeSubdivision) visitTriangles(triVisitor func(triEdges []*QuadE
// visited flag is used to record visited edges of triangles
// setVisitedAll(false);
var stack *edgeStack = new(edgeStack)
- stack.push(qes.startingEdge)
+ if qes.startingEdge != nil {
+ stack.push(qes.startingEdge)
+ }
visitedEdges := make(edgeSet)
@@ -718,10 +747,14 @@ func (qes *QuadEdgeSubdivision) fetchTriangleToVisit(edge *QuadEdge, stack *edge
triEdges := make([]*QuadEdge, 0, 3)
curr := edge
var isFrame bool
- var done bool
- for !done {
+
+ for true {
triEdges = append(triEdges, curr)
+ if curr.IsLive() == false {
+ log.Fatal("traversing dead edge")
+ }
+
if qes.isFrameEdge(curr) {
isFrame = true
}
@@ -738,7 +771,7 @@ func (qes *QuadEdgeSubdivision) fetchTriangleToVisit(edge *QuadEdge, stack *edge
curr = curr.LNext()
if curr == edge {
- done = true
+ break
}
}
diff --git a/planar/triangulate/quadedge/quadedgesubdivision_test.go b/planar/triangulate/quadedge/quadedgesubdivision_test.go
index e3ca3b49..de489eb3 100644
--- a/planar/triangulate/quadedge/quadedgesubdivision_test.go
+++ b/planar/triangulate/quadedge/quadedgesubdivision_test.go
@@ -77,6 +77,12 @@ func TestQuadEdgeSubdivisionDelete(t *testing.T) {
if err != nil {
t.Errorf("expected nil got %v", err)
}
+
+ _, err = uut.LocateSegment(Vertex{100,1000}, tc.b)
+ if err == nil {
+ t.Errorf("expected %v got %v", ErrLocateFailure, err)
+ }
+
if qe.Orig().Equals(tc.a) == false || qe.Dest().Equals(tc.b) == false {
t.Errorf("expected true got false")
}
@@ -166,6 +172,12 @@ func TestQuadEdgeSubdivisionIsOnEdge(t *testing.T) {
if onEdge != tc.expected {
t.Fatalf("expected %v got %v", tc.expected, onEdge)
}
+
+ onLine := uut.IsOnLine(geom.Line{tc.e1, tc.e2}, tc.p)
+
+ if onLine != tc.expected {
+ t.Fatalf("expected %v got %v", tc.expected, onEdge)
+ }
}
testcases := []tcase{
{Vertex{0, 0}, Vertex{5, 5}, Vertex{3, 3}, true},
diff --git a/planar/triangulate/quadedge/vertex.go b/planar/triangulate/quadedge/vertex.go
index c7acda40..7640ac0f 100644
--- a/planar/triangulate/quadedge/vertex.go
+++ b/planar/triangulate/quadedge/vertex.go
@@ -120,13 +120,23 @@ func (u Vertex) Dot(v Vertex) float64 {
/*
Times computes the scalar product c(v)
-@param v a vertex
-@return returns the scaled vector
+Return the scaled vector
*/
func (u Vertex) Times(c float64) Vertex {
return Vertex{u.X() * c, u.Y() * c}
}
+/*
+Divide computes the scalar division v / c
+
+Returns the scaled vector
+
+This is not part of the original JTS code.
+*/
+func (u Vertex) Divide(c float64) Vertex {
+ return Vertex{u.X() / c, u.Y() / c}
+}
+
// Sum u + v and return the new Vertex
func (u Vertex) Sum(v Vertex) Vertex {
return Vertex{u.X() + v.X(), u.Y() + v.Y()}
@@ -142,6 +152,15 @@ func (u Vertex) Magn() float64 {
return math.Sqrt(u.X()*u.X() + u.Y()*u.Y())
}
+/*
+Normalize scales the vector so the length is one.
+
+This is not part of the original JTS code.
+*/
+func (u Vertex) Normalize() Vertex {
+ return u.Divide(u.Magn())
+}
+
/*
Cross returns k X v (cross product). this is a vector perpendicular to v
*/
diff --git a/planar/triangulate/quadedge/vertex_test.go b/planar/triangulate/quadedge/vertex_test.go
index 466b0ce6..60e3bf5d 100644
--- a/planar/triangulate/quadedge/vertex_test.go
+++ b/planar/triangulate/quadedge/vertex_test.go
@@ -157,6 +157,7 @@ func TestVertexScalar(t *testing.T) {
v Vertex
scalar float64
times Vertex
+ divide Vertex
}
fn := func(t *testing.T, tc tcase) {
@@ -166,16 +167,82 @@ func TestVertexScalar(t *testing.T) {
return
}
- r = tc.v.Cross()
- c := Vertex{tc.v.Y(), -tc.v.X()}
- if c.Equals(r) == false {
- t.Errorf("error, expected %v got %v", c, r)
+ r = tc.v.Divide(tc.scalar)
+ if r.Equals(tc.divide) == false {
+ t.Errorf("error, expected %v got %v", tc.divide, r)
+ return
+ }
+ }
+ testcases := []tcase{
+ {Vertex{1, 2}, 3, Vertex{3, 6}, Vertex{0.3333333333333333, 0.6666666666666666}},
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}
+
+func TestVertexUnary(t *testing.T) {
+ type tcase struct {
+ v Vertex
+ cross Vertex
+ magn float64
+ normalize Vertex
+ }
+
+ fn := func(t *testing.T, tc tcase) {
+ r := tc.v.Cross()
+ if tc.cross.Equals(r) == false {
+ t.Errorf("error, expected %v got %v", tc.cross, r)
+ return
+ }
+
+ m := tc.v.Magn()
+ if tc.magn != m {
+ t.Errorf("error, expected %v got %v", tc.magn, m)
+ return
+ }
+
+ r = tc.v.Normalize()
+ if tc.normalize.Equals(r) == false {
+ t.Errorf("error, expected %v got %v", tc.normalize, r)
+ return
+ }
+ }
+ testcases := []tcase{
+ {Vertex{1, 2}, Vertex{2, -1}, 2.23606797749979, Vertex{0.4472135954999579, 0.8944271909999159}},
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}
+
+func TestVertexVertex(t *testing.T) {
+ type tcase struct {
+ v1 Vertex
+ v2 Vertex
+ dot float64
+ sum Vertex
+ }
+
+ fn := func(t *testing.T, tc tcase) {
+ s := tc.v1.Dot(tc.v2)
+ if s != tc.dot {
+ t.Errorf("error, expected %v got %v", tc.dot, s)
return
}
+ r := tc.v1.Sum(tc.v2)
+ if r.Equals(tc.sum) == false {
+ t.Errorf("error, expected %v got %v", tc.sum, r)
+ return
+ }
}
testcases := []tcase{
- {Vertex{1, 2}, 3, Vertex{3, 6}},
+ {Vertex{1, 2}, Vertex{3, 4}, 11, Vertex{4, 6}},
}
for i, tc := range testcases {
diff --git a/planar/triangulate/segment.go b/planar/triangulate/segment.go
new file mode 100644
index 00000000..2c77e3df
--- /dev/null
+++ b/planar/triangulate/segment.go
@@ -0,0 +1,192 @@
+/*
+Copyright (c) 2016 Vivid Solutions.
+
+All rights reserved. This program and the accompanying materials
+are made available under the terms of the Eclipse Public License v1.0
+and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v10.html
+and the Eclipse Distribution License is available at
+
+http://www.eclipse.org/org/documents/edl-v10.php.
+*/
+
+package triangulate
+
+import (
+ "github.com/go-spatial/geom"
+ "github.com/go-spatial/geom/planar/triangulate/quadedge"
+)
+
+/*
+Segment models a constraint segment in a triangulation.
+
+A constraint segment is an oriented straight line segment between a start point
+and an end point.
+
+Author David Skea
+Author Martin Davis
+Ported to Go by Jason R. Surratt
+*/
+type Segment struct {
+ ls geom.Line
+ data interface{}
+}
+
+func NewSegment(l geom.Line) Segment {
+ return Segment{ls: l}
+}
+
+/**
+ * Creates a new instance for the given ordinates.
+ public Segment(double x1, double y1, double z1, double x2, double y2, double z2) {
+ this(new Coordinate(x1, y1, z1), new Coordinate(x2, y2, z2));
+ }
+*/
+
+/**
+ * Creates a new instance for the given ordinates, with associated external data.
+ public Segment(double x1, double y1, double z1, double x2, double y2, double z2, Object data) {
+ this(new Coordinate(x1, y1, z1), new Coordinate(x2, y2, z2), data);
+ }
+*/
+
+/**
+ * Creates a new instance for the given points.
+ *
+ * @param p0 the start point
+ * @param p1 the end point
+ public Segment(Coordinate p0, Coordinate p1) {
+ ls = new LineSegment(p0, p1);
+ }
+*/
+
+/*
+Gets the start coordinate of the segment
+
+Returns the starting vertex
+*/
+func (seg *Segment) GetStart() quadedge.Vertex {
+ return quadedge.Vertex(seg.ls[0])
+}
+
+/*
+Gets the end coordinate of the segment
+
+Return a Coordinate
+*/
+func (seg *Segment) GetEnd() quadedge.Vertex {
+ return quadedge.Vertex(seg.ls[1])
+}
+
+/**
+ * Gets the start X ordinate of the segment
+ *
+ * @return the X ordinate value
+ public double getStartX() {
+ Coordinate p = ls.getCoordinate(0);
+ return p.x;
+ }
+*/
+
+/**
+ * Gets the start Y ordinate of the segment
+ *
+ * @return the Y ordinate value
+ public double getStartY() {
+ Coordinate p = ls.getCoordinate(0);
+ return p.y;
+ }
+*/
+
+/**
+ * Gets the start Z ordinate of the segment
+ *
+ * @return the Z ordinate value
+ public double getStartZ() {
+ Coordinate p = ls.getCoordinate(0);
+ return p.z;
+ }
+*/
+
+/**
+ * Gets the end X ordinate of the segment
+ *
+ * @return the X ordinate value
+ public double getEndX() {
+ Coordinate p = ls.getCoordinate(1);
+ return p.x;
+ }
+*/
+
+/**
+ * Gets the end Y ordinate of the segment
+ *
+ * @return the Y ordinate value
+ public double getEndY() {
+ Coordinate p = ls.getCoordinate(1);
+ return p.y;
+ }
+*/
+
+/**
+ * Gets the end Z ordinate of the segment
+ *
+ * @return the Z ordinate value
+ public double getEndZ() {
+ Coordinate p = ls.getCoordinate(1);
+ return p.z;
+ }
+*/
+
+// GetLineSegment gets a Line modelling this segment.
+func (seg *Segment) GetLineSegment() geom.Line {
+ return seg.ls
+}
+
+/**
+ * Gets the external data associated with this segment
+ *
+ * @return a data object
+ public Object getData() {
+ return data;
+ }
+*/
+
+/**
+ * Sets the external data to be associated with this segment
+ *
+ * @param data a data object
+ public void setData(Object data) {
+ this.data = data;
+ }
+*/
+
+/**
+ * Determines whether two segments are topologically equal.
+ * I.e. equal up to orientation.
+ *
+ * @param s a segment
+ * @return true if the segments are topologically equal
+ public boolean equalsTopo(Segment s) {
+ return ls.equalsTopo(s.getLineSegment());
+ }
+*/
+
+/**
+ * Computes the intersection point between this segment and another one.
+ *
+ * @param s a segment
+ * @return the intersection point, or null if there is none
+ public Coordinate intersection(Segment s) {
+ return ls.intersection(s.getLineSegment());
+ }
+*/
+
+/**
+ * Computes a string representation of this segment.
+ *
+ * @return a string
+ public String toString() {
+ return ls.toString();
+ }
+*/
diff --git a/planar/triangulate/segment_test.go b/planar/triangulate/segment_test.go
new file mode 100644
index 00000000..763b2cc3
--- /dev/null
+++ b/planar/triangulate/segment_test.go
@@ -0,0 +1,52 @@
+/*
+Copyright (c) 2016 Vivid Solutions.
+
+All rights reserved. This program and the accompanying materials
+are made available under the terms of the Eclipse Public License v1.0
+and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v10.html
+and the Eclipse Distribution License is available at
+
+http://www.eclipse.org/org/documents/edl-v10.php.
+*/
+
+package triangulate
+
+import (
+ "strconv"
+ "testing"
+
+ "github.com/go-spatial/geom"
+)
+
+/*
+TestSegmentDummy keeps coveralls from complaining. Not really necessary tests.
+*/
+func TestSegmentDummy(t *testing.T) {
+ type tcase struct {
+ line geom.Line
+ }
+
+ fn := func(t *testing.T, tc tcase) {
+ s := NewSegment(tc.line)
+ if s.GetStart().Equals(tc.line[0]) == false {
+ t.Errorf("error, expected %v got %v", tc.line[0], s.GetStart())
+ }
+ if s.GetEnd().Equals(tc.line[1]) == false {
+ t.Errorf("error, expected %v got %v", tc.line[1], s.GetEnd())
+ }
+ if s.GetLineSegment() != tc.line {
+ t.Errorf("error, expected %v got %v", tc.line, s.GetLineSegment())
+ }
+ }
+ testcases := []tcase{
+ {
+ line: geom.Line{{1, 2}, {3, 4}},
+ },
+ }
+
+ for i, tc := range testcases {
+ tc := tc
+ t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) })
+ }
+}