Skip to content

Commit

Permalink
Merge pull request #339 from peterstace/use_spanning_tree_for_ghosts
Browse files Browse the repository at this point in the history
Use Spanning Tree for Ghost Lines
  • Loading branch information
peterstace authored Jan 29, 2021
2 parents c13a053 + 08b6aec commit bc4715a
Show file tree
Hide file tree
Showing 6 changed files with 259 additions and 261 deletions.
9 changes: 5 additions & 4 deletions geom/alg_binary_op.go
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,7 @@ func binaryOp(a, b Geometry, include func(uint8) bool) (Geometry, error) {
if err != nil {
return Geometry{}, fmt.Errorf("internal error creating overlay: %v", err)
}

g, err := overlay.extractGeometry(include)
if err != nil {
return Geometry{}, fmt.Errorf("internal error extracting geometry: %v", err)
Expand All @@ -78,10 +79,10 @@ func createOverlay(a, b Geometry) (*doublyConnectedEdgeList, error) {
return nil, errors.New("GeometryCollection argument not supported")
}

aGhost := connectGeometry(a)
bGhost := connectGeometry(b)
joinGhost := connectGeometries(a, b)
ghosts := mergeMultiLineStrings([]MultiLineString{aGhost, bGhost, joinGhost.AsMultiLineString()})
var points []XY
points = appendComponentPoints(points, a)
points = appendComponentPoints(points, b)
ghosts := spanningTree(points)

a, b, ghosts, err := reNodeGeometries(a, b, ghosts)
if err != nil {
Expand Down
19 changes: 12 additions & 7 deletions geom/alg_binary_op_test.go
Original file line number Diff line number Diff line change
Expand Up @@ -137,11 +137,11 @@ func TestBinaryOp(t *testing.T) {
*/
input1: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,2 3,2 1,0 1)))",
input2: "MULTIPOLYGON(((4 0,4 1,5 1,5 0,4 0)),((1 0,1 2,3 2,3 0,1 0)))",
union: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,1 3,2 3,2 2,3 2,3 1,3 0,1 0,1 1,0 1)),((4 0,4 1,5 1,5 0,4 0)))",
union: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,2 3,2 2,3 2,3 0,1 0,1 1,0 1)),((4 0,4 1,5 1,5 0,4 0)))",
inter: "POLYGON((2 2,2 1,1 1,1 2,2 2))",
fwdDiff: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,1 3,2 3,2 2,1 2,1 1,0 1)))",
revDiff: "MULTIPOLYGON(((4 0,4 1,5 1,5 0,4 0)),((1 0,1 1,2 1,2 2,3 2,3 1,3 0,1 0)))",
symDiff: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,1 3,2 3,2 2,1 2,1 1,0 1)),((1 1,2 1,2 2,3 2,3 1,3 0,1 0,1 1)),((4 0,4 1,5 1,5 0,4 0)))",
fwdDiff: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,2 3,2 2,1 2,1 1,0 1)))",
revDiff: "MULTIPOLYGON(((4 0,4 1,5 1,5 0,4 0)),((1 0,1 1,2 1,2 2,3 2,3 0,1 0)))",
symDiff: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,2 3,2 2,1 2,1 1,0 1)),((1 1,2 1,2 2,3 2,3 0,1 0,1 1)),((4 0,4 1,5 1,5 0,4 0)))",
},
{
/*
Expand Down Expand Up @@ -722,11 +722,11 @@ func TestBinaryOp(t *testing.T) {
{
input1: "LINESTRING(0 1,1 0)",
input2: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))",
union: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))",
union: "MULTIPOLYGON(((0 0,0 1,1 1,1 0.5,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))",
inter: "LINESTRING(0 1,1 0)",
fwdDiff: "GEOMETRYCOLLECTION EMPTY",
revDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))",
symDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))",
revDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0.5,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))",
symDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0.5,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))",
},
{
input1: "POLYGON((1 0,0 1,1 1,1 0))",
Expand Down Expand Up @@ -805,6 +805,11 @@ func TestBinaryOp(t *testing.T) {
input2: "MULTIPOLYGON(((0 0,2 0,2 2,0 2,0 0),(0.5 0.5,1 0.5,1 1.5,0.5 1.5,0.5 0.5)))",
union: "POLYGON((0 0,1 0,2 0,2 2,0 2,0 1,0 0),(0.5000000000000001 0.5,1 0.5,1 1.5,0.5 1.5,0.5000000000000001 0.5))",
},
{
input1: "MULTILINESTRING((2 0,2 1),(2 2,2 3))",
input2: "POLYGON((0 0,0 10,10 10,10 0,0 0),(1.5 1.5,8.5 1.5,8.5 8.5,1.5 8.5,1.5 1.5))",
union: "GEOMETRYCOLLECTION(POLYGON((2 0,10 0,10 10,0 10,0 0,2 0),(1.5 1.5,1.5 8.5,8.5 8.5,8.5 1.5,1.5 1.5)),LINESTRING(2 2,2 3))",
},
} {
t.Run(strconv.Itoa(i), func(t *testing.T) {
g1 := geomFromWKT(t, geomCase.input1)
Expand Down
194 changes: 82 additions & 112 deletions geom/dcel_ghosts.go
Original file line number Diff line number Diff line change
@@ -1,152 +1,122 @@
package geom

import "fmt"
import (
"fmt"

func connectGeometry(g Geometry) MultiLineString {
var ghostLSs []LineString
var seenFirst bool
var first XY
addComponent := func(pt Point) {
xy, ok := pt.XY()
if !ok {
return
}
if seenFirst {
if first != xy {
seq := NewSequence([]float64{first.X, first.Y, xy.X, xy.Y}, DimXY)
ghostLS, err := NewLineString(seq)
if err != nil {
// Can't happen, since first and pt are not the same.
panic(fmt.Sprintf("could not construct LineString: %v", err))
}
ghostLSs = append(ghostLSs, ghostLS)
}
} else {
seenFirst = true
first = xy
}
"github.com/peterstace/simplefeatures/rtree"
)

// spanningTree creates a near-minimum spanning tree (using the euclidean
// distance metric) over the supplied points. The tree will consist of N-1
// lines, where N is the number of _distinct_ xys supplied.
//
// It's a 'near' minimum spanning tree rather than a spanning tree, because we
// use a simple greedy algorithm rather than a proper minimum spanning tree
// algorithm.
func spanningTree(xys []XY) MultiLineString {
if len(xys) <= 1 {
return MultiLineString{}
}

switch g.Type() {
case TypePoint:
// A single Point is already trivially connected.
case TypeMultiPoint:
mp := g.AsMultiPoint()
n := mp.NumPoints()
for i := 0; i < n; i++ {
addComponent(mp.PointN(i))
}
case TypeLineString:
// LineStrings are already connected.
case TypeMultiLineString:
mls := g.AsMultiLineString()
n := mls.NumLineStrings()
for i := 0; i < n; i++ {
ls := mls.LineStringN(i)
addComponent(ls.StartPoint())
}
case TypePolygon:
poly := g.AsPolygon()
addComponent(poly.ExteriorRing().StartPoint())
n := poly.NumInteriorRings()
for i := 0; i < n; i++ {
addComponent(poly.InteriorRingN(i).StartPoint())
// Load points into r-tree.
xys = sortAndUniquifyXYs(xys)
items := make([]rtree.BulkItem, len(xys))
for i, xy := range xys {
items[i] = rtree.BulkItem{Box: xy.box(), RecordID: i}
}
tree := rtree.BulkLoad(items)

// The disjoint set keeps track of which points have been joined together
// so far. Two entries in dset are in the same set iff they are connected
// in the incrementally-built spanning tree.
dset := newDisjointSet(len(xys))
lss := make([]LineString, 0, len(xys)-1)

for i, xyi := range xys {
if i == len(xys)-1 {
// Skip the last point, since a tree is formed from N-1 edges
// rather than N edges. The last point will be included by virtue
// of being the closest to another point.
continue
}
case TypeMultiPolygon:
mp := g.AsMultiPolygon()
n := mp.NumPolygons()
for i := 0; i < n; i++ {
poly := mp.PolygonN(i)
addComponent(poly.ExteriorRing().StartPoint())
m := poly.NumInteriorRings()
for j := 0; j < m; j++ {
addComponent(poly.InteriorRingN(j).StartPoint())
tree.PrioritySearch(xyi.box(), func(j int) error {
// We don't want to include a new edge in the spanning tree if it
// would cause a cycle (i.e. the two endpoints are already in the
// same tree). This is checked via dset.
if i == j || dset.find(i) == dset.find(j) {
return nil
}
}
case TypeGeometryCollection:
gc := g.AsGeometryCollection()
n := gc.NumGeometries()
for i := 0; i < n; i++ {
addComponent(pointOnGeometry(gc.GeometryN(i)))
}
default:
panic(fmt.Sprintf("unknown geometry type: %v", g.Type()))
dset.union(i, j)
xyj := xys[j]
lss = append(lss, line{xyi, xyj}.asLineString())
return rtree.Stop
})
}

return NewMultiLineStringFromLineStrings(ghostLSs)
return NewMultiLineStringFromLineStrings(lss)
}

func connectGeometries(g1, g2 Geometry) LineString {
pt1 := pointOnGeometry(g1)
pt2 := pointOnGeometry(g2)

xy1, ok1 := pt1.XY()
xy2, ok2 := pt2.XY()
if !ok1 || !ok2 || xy1 == xy2 {
return LineString{}
func appendXYForPoint(xys []XY, pt Point) []XY {
if xy, ok := pt.XY(); ok {
xys = append(xys, xy)
}
return xys
}

func appendXYForLineString(xys []XY, ls LineString) []XY {
return appendXYForPoint(xys, ls.StartPoint())
}

coords := []float64{xy1.X, xy1.Y, xy2.X, xy2.Y}
ls, err := NewLineString(NewSequence(coords, DimXY))
if err != nil {
// Can't happen, since we have already checked that xy1 != xy2.
panic(fmt.Sprintf("could not create lines string: %v", err))
func appendXYsForPolygon(xys []XY, poly Polygon) []XY {
xys = appendXYForLineString(xys, poly.ExteriorRing())
n := poly.NumInteriorRings()
for i := 0; i < n; i++ {
xys = appendXYForLineString(xys, poly.InteriorRingN(i))
}
return ls
return xys
}

func pointOnGeometry(g Geometry) Point {
func appendComponentPoints(xys []XY, g Geometry) []XY {
switch g.Type() {
case TypePoint:
return g.AsPoint()
return appendXYForPoint(xys, g.AsPoint())
case TypeMultiPoint:
mp := g.AsMultiPoint()
n := mp.NumPoints()
for i := 0; i < n; i++ {
pt := mp.PointN(i)
if !pt.IsEmpty() {
return pt
}
xys = appendXYForPoint(xys, mp.PointN(i))
}
return Point{}
return xys
case TypeLineString:
return g.AsLineString().StartPoint()
ls := g.AsLineString()
return appendXYForLineString(xys, ls)
case TypeMultiLineString:
mls := g.AsMultiLineString()
n := mls.NumLineStrings()
for i := 0; i < n; i++ {
pt := mls.LineStringN(i).StartPoint()
if !pt.IsEmpty() {
return pt
}
ls := mls.LineStringN(i)
xys = appendXYForLineString(xys, ls)
}
return Point{}
return xys
case TypePolygon:
return pointOnGeometry(g.Boundary())
poly := g.AsPolygon()
return appendXYsForPolygon(xys, poly)
case TypeMultiPolygon:
return pointOnGeometry(g.Boundary())
mp := g.AsMultiPolygon()
n := mp.NumPolygons()
for i := 0; i < n; i++ {
poly := mp.PolygonN(i)
xys = appendXYsForPolygon(xys, poly)
}
return xys
case TypeGeometryCollection:
gc := g.AsGeometryCollection()
n := gc.NumGeometries()
for i := 0; i < n; i++ {
pt := pointOnGeometry(gc.GeometryN(i))
if !pt.IsEmpty() {
return pt
}
xys = appendComponentPoints(xys, gc.GeometryN(i))
}
return Point{}
return xys
default:
panic(fmt.Sprintf("unknown geometry type: %v", g.Type()))
}
}

func mergeMultiLineStrings(mlss []MultiLineString) MultiLineString {
var lss []LineString
for _, mls := range mlss {
n := mls.NumLineStrings()
for i := 0; i < n; i++ {
lss = append(lss, mls.LineStringN(i))
}
}
return NewMultiLineStringFromLineStrings(lss)
}
47 changes: 47 additions & 0 deletions geom/dcel_ghosts_test.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
package geom

import (
"strconv"
"testing"
)

func TestSpanningTree(t *testing.T) {
for i, tc := range []struct {
xys []XY
wantWKT string
}{
{
xys: nil,
wantWKT: "MULTILINESTRING EMPTY",
},
{
xys: []XY{{1, 1}},
wantWKT: "MULTILINESTRING EMPTY",
},
{
xys: []XY{{2, 1}, {1, 2}},
wantWKT: "MULTILINESTRING((2 1,1 2))",
},
{
xys: []XY{{2, 0}, {2, 2}, {0, 0}, {1.5, 1.5}},
wantWKT: "MULTILINESTRING((0 0,2 0),(1.5 1.5,2 2),(2 0,1.5 1.5))",
},
{
xys: []XY{{-0.5, 0.5}, {0, 0}, {0, 1}, {1, 0}},
wantWKT: "MULTILINESTRING((-0.5 0.5,0 0),(0 0,0 1),(0 1,1 0))",
},
} {
t.Run(strconv.Itoa(i), func(t *testing.T) {
want, err := UnmarshalWKT(tc.wantWKT)
if err != nil {
t.Fatal(err)
}
got := spanningTree(tc.xys)
if !ExactEquals(want, got.AsGeometry(), IgnoreOrder) {
t.Logf("got: %v", got.AsText())
t.Logf("want: %v", want.AsText())
t.Fatal("mismatch")
}
})
}
}
Loading

0 comments on commit bc4715a

Please sign in to comment.