diff --git a/planar/makevalid/hitmap/hitmap.go b/planar/makevalid/hitmap/hitmap.go new file mode 100644 index 00000000..d8480d47 --- /dev/null +++ b/planar/makevalid/hitmap/hitmap.go @@ -0,0 +1,142 @@ +package hitmap + +import ( + "errors" + "sort" + + "github.com/go-spatial/geom" +) + +var ErrInvalidLineString = errors.New("invalid linestring") + +// label is the he label for the triangle. Is in "inside" or "outside". +// TODO: gdey — would be make more sense to just have a bool here? IsInside or somthing like that? +type Label uint8 + +func (l Label) String() string { + switch l { + case Outside: + return "outside" + case Inside: + return "inside" + default: + return "unknown" + } +} + +const ( + Unknown Label = iota + Outside + Inside +) + +type Interface interface { + LabelFor(pt [2]float64) Label + Extent() [4]float64 + Area() float64 +} + +func asGeomExtent(e [4]float64) *geom.Extent { + ee := geom.Extent(e) + return &ee + +} + +// PolygonHMSliceByAreaDec will allow you to sort a slice of PolygonHM in decending order +type ByAreaDec []Interface + +func (hm ByAreaDec) Len() int { return len(hm) } +func (hm ByAreaDec) Swap(i, j int) { hm[i], hm[j] = hm[j], hm[i] } +func (hm ByAreaDec) Less(i, j int) bool { + ia, ja := hm[i].Area(), hm[j].Area() + return ia < ja +} + +// OrderedHM will iterate through a set of HitMaps looking for the first one to return +// inside, if none of the hitmaps return inside it will return outside. +type OrderedHM []Interface + +func (hms OrderedHM) LabelFor(pt [2]float64) Label { + for i := range hms { + if hms[i].LabelFor(pt) == Inside { + return Inside + } + } + return Outside +} + +// Extent is the accumlative extent of all the extens in the slice. +func (hms OrderedHM) Extent() [4]float64 { + e := new(geom.Extent) + for i := range hms { + e.Add(asGeomExtent(hms[i].Extent())) + } + return e.Extent() +} + +// Area returns the area of the total extent of the hitmaps that are contain in the slice. +func (hms OrderedHM) Area() float64 { + return asGeomExtent(hms.Extent()).Area() +} + +// NewOrderdHM will add the provided hitmaps in reverse order so that the last hit map is always tried first. +func NewOrderedHM(hms ...Interface) OrderedHM { + ohm := make(OrderedHM, len(hms)) + size := len(hms) - 1 + for i := size; i >= 0; i-- { + ohm[size-i] = hms[i] + } + return ohm +} + +// NewHitMap will return a Polygon Hit map, a Ordered Hit Map, or a nil Hit map based on the geomtry type. +func New(clipbox *geom.Extent, geo geom.Geometry) (Interface, error) { + var err error + switch g := geo.(type) { + case geom.Polygoner: + + ghm, err := NewFromPolygons(clipbox, g.LinearRings()) + if err != nil { + return nil, err + } + + return ghm, nil + + case geom.MultiPolygoner: + + polygons := g.Polygons() + ghms := make([]Interface, len(polygons)) + for i := range polygons { + ghms[i], err = NewFromPolygons(clipbox, polygons[i]) + if err != nil { + return nil, err + } + } + sort.Sort(ByAreaDec(ghms)) + return NewOrderedHM(ghms...), nil + + case geom.Collectioner: + + geometries := g.Geometries() + ghms := make([]Interface, 0, len(geometries)) + for i := range geometries { + g, err := New(clipbox, geometries[i]) + if err != nil { + return nil, err + } + // skip empty hitmaps. + if g == nil { + continue + } + ghms = append(ghms, g) + } + sort.Sort(ByAreaDec(ghms)) + return NewOrderedHM(ghms...), nil + + case geom.Pointer, geom.MultiPointer, geom.LineStringer, geom.MultiLineStringer: + return nil, nil + + default: + return nil, geom.ErrUnknownGeometry{geo} + } +} diff --git a/planar/makevalid/hitmap/hitmap_test.go b/planar/makevalid/hitmap/hitmap_test.go new file mode 100644 index 00000000..66f26529 --- /dev/null +++ b/planar/makevalid/hitmap/hitmap_test.go @@ -0,0 +1,164 @@ +package hitmap + +import ( + "strconv" + "testing" + + "github.com/go-spatial/geom" +) + +func TestRingContains(t *testing.T) { + type pt struct { + pt [2]float64 + contained bool + } + cpt := func(x, y float64) pt { return pt{pt: [2]float64{x, y}, contained: true} } + opt := func(x, y float64) pt { return pt{pt: [2]float64{x, y}, contained: false} } + type tcase struct { + linestring [][2]float64 + pts []pt + } + fn := func(t *testing.T, tc tcase) { + ring, err := NewRing(tc.linestring, Outside) + if err != nil { + // panic for now, should have test for error generated by create segments. + panic(err) + } + for i := range tc.pts { + i := i + t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { + c := ring.Contains(tc.pts[i].pt) + if c != tc.pts[i].contained { + t.Errorf("containes %v, expected %v got %v", tc.pts[i].pt, tc.pts[i].contained, c) + } + }) + } + } + tests := map[string]tcase{ + "simple 6 seg": { + linestring: [][2]float64{{1, 1}, {4, -4}, {8, -4}, {8, 5}, {3, 5}, {1, 3}}, + pts: []pt{ + cpt(2, 2), cpt(2, 3), opt(2, 5), opt(2, -4), opt(6, -4), opt(8, 0), + opt(20, -4), opt(2, -2), cpt(4, -1), cpt(3, 3)}, + }, + "complicated shape 20x20": { + linestring: [][2]float64{ + {2, 3}, {4, 3}, {4, 4}, {6, 6}, {9, 6}, {8, 4}, {6, 4}, + {4, 2}, {10, 2}, {10, 4}, {12, 6}, {16, 3}, {16, 4}, + {18, 6}, {18, 8}, {16, 12}, {14, 10}, {16, 8}, {16, 6}, + {12, 11}, {10, 8}, {10, 7}, {8, 7}, {8, 10}, {6, 10}, + {6, 8}, {4, 8}, {4, 12}, {18, 18}, {8, 18}, {2, 12}, + {2, 8}, {4, 6}, {2, 4}, + }, + pts: []pt{ + opt(1, 1), opt(1, 2), opt(1, 3), opt(1, 4), opt(1, 5), opt(1, 6), opt(1, 7), opt(1, 8), opt(1, 9), opt(1, 10), + opt(1, 11), opt(1, 12), opt(1, 13), opt(1, 14), opt(1, 15), opt(1, 16), opt(1, 17), opt(1, 18), opt(1, 19), opt(1, 20), + + opt(2, 1), opt(2, 2), opt(2, 3), opt(2, 4), opt(2, 5), opt(2, 6), opt(2, 7), opt(2, 8), opt(2, 9), opt(2, 10), + opt(2, 11), opt(2, 12), opt(2, 13), opt(2, 14), opt(2, 15), opt(2, 16), opt(2, 17), opt(2, 18), opt(2, 19), opt(2, 20), + + opt(3, 1), opt(3, 2), opt(3, 3), cpt(3, 4), opt(3, 5), opt(3, 6), opt(3, 7), cpt(3, 8), cpt(3, 9), cpt(3, 10), + cpt(3, 11), cpt(3, 12), opt(3, 13), opt(3, 14), opt(3, 15), opt(3, 16), opt(3, 17), opt(3, 18), opt(3, 19), opt(3, 20), + + opt(4, 1), opt(4, 2), opt(4, 3), opt(4, 4), cpt(4, 5), opt(4, 6), cpt(4, 7), opt(4, 8), opt(4, 9), opt(4, 10), + opt(4, 11), opt(4, 12), cpt(4, 13), opt(4, 14), opt(4, 15), opt(4, 16), opt(4, 17), opt(4, 18), opt(4, 19), opt(4, 20), + + opt(5, 1), opt(5, 2), opt(5, 3), opt(5, 4), opt(5, 5), cpt(5, 6), cpt(5, 7), opt(5, 8), opt(5, 9), opt(5, 10), + opt(5, 11), opt(5, 12), cpt(5, 13), cpt(5, 14), opt(5, 15), opt(5, 16), opt(5, 17), opt(5, 18), opt(5, 19), opt(5, 20), + + opt(6, 1), opt(6, 2), cpt(6, 3), opt(6, 4), opt(6, 5), opt(6, 6), cpt(6, 7), opt(6, 8), opt(6, 9), opt(6, 10), + opt(6, 11), opt(6, 12), cpt(6, 13), cpt(6, 14), cpt(6, 15), opt(6, 16), opt(6, 17), opt(6, 18), opt(6, 19), opt(6, 20), + + opt(7, 1), opt(7, 2), cpt(7, 3), opt(7, 4), opt(7, 5), opt(7, 6), cpt(7, 7), cpt(7, 8), cpt(7, 9), opt(7, 10), + opt(7, 11), opt(7, 12), opt(7, 13), cpt(7, 14), cpt(7, 15), cpt(7, 16), opt(7, 17), opt(7, 18), opt(7, 19), opt(7, 20), + + opt(8, 1), opt(8, 2), cpt(8, 3), opt(8, 4), opt(8, 5), opt(8, 6), opt(8, 7), opt(8, 8), opt(8, 9), opt(8, 10), + opt(8, 11), opt(8, 12), opt(8, 13), cpt(8, 14), cpt(8, 15), cpt(8, 16), cpt(8, 17), opt(8, 18), opt(8, 19), opt(8, 20), + + opt(9, 1), opt(9, 2), cpt(9, 3), cpt(9, 4), cpt(9, 5), opt(9, 6), opt(9, 7), opt(9, 8), opt(9, 9), opt(9, 10), + opt(9, 11), opt(9, 12), opt(9, 13), opt(9, 14), cpt(9, 15), cpt(9, 16), cpt(9, 17), opt(9, 18), opt(9, 19), opt(9, 20), + + opt(10, 1), opt(10, 2), opt(10, 3), opt(10, 4), cpt(10, 5), cpt(10, 6), opt(10, 7), opt(10, 8), opt(10, 9), opt(10, 10), + opt(10, 11), opt(10, 12), opt(10, 13), opt(10, 14), cpt(10, 15), cpt(10, 16), cpt(10, 17), opt(10, 18), opt(10, 19), opt(10, 20), + + opt(11, 1), opt(11, 2), opt(11, 3), opt(11, 4), opt(11, 5), cpt(11, 6), cpt(11, 7), cpt(11, 8), cpt(11, 9), opt(11, 10), + opt(11, 11), opt(11, 12), opt(11, 13), opt(11, 14), opt(11, 15), cpt(11, 16), cpt(11, 17), opt(11, 18), opt(11, 19), opt(11, 20), + + opt(12, 1), opt(12, 2), opt(12, 3), opt(12, 4), opt(12, 5), opt(12, 6), cpt(12, 7), cpt(12, 8), cpt(12, 9), cpt(12, 10), + opt(12, 11), opt(12, 12), opt(12, 13), opt(12, 14), opt(12, 15), cpt(12, 16), cpt(12, 17), opt(12, 18), opt(12, 19), opt(12, 20), + + opt(13, 1), opt(13, 2), opt(13, 3), opt(13, 4), opt(13, 5), cpt(13, 6), cpt(13, 7), cpt(13, 8), cpt(13, 9), opt(13, 10), + opt(13, 11), opt(13, 12), opt(13, 13), opt(13, 14), opt(13, 15), cpt(13, 16), cpt(13, 17), opt(13, 18), opt(13, 19), opt(13, 20), + + opt(14, 1), opt(14, 2), opt(14, 3), opt(14, 4), cpt(14, 5), cpt(14, 6), cpt(14, 7), cpt(14, 8), opt(14, 9), opt(14, 10), + opt(14, 11), opt(14, 12), opt(14, 13), opt(14, 14), opt(14, 15), opt(14, 16), cpt(14, 17), opt(14, 18), opt(14, 19), opt(14, 20), + + opt(15, 1), opt(15, 2), opt(15, 3), cpt(15, 4), cpt(15, 5), cpt(15, 6), cpt(15, 7), opt(15, 8), opt(15, 9), cpt(15, 10), + opt(15, 11), opt(15, 12), opt(15, 13), opt(15, 14), opt(15, 15), opt(15, 16), cpt(15, 17), opt(15, 18), opt(15, 19), opt(15, 20), + + opt(16, 1), opt(16, 2), opt(16, 3), opt(16, 4), cpt(16, 5), opt(16, 6), opt(16, 7), opt(16, 8), cpt(16, 9), cpt(16, 10), + cpt(16, 11), opt(16, 12), opt(16, 13), opt(16, 14), opt(16, 15), opt(16, 16), opt(16, 17), opt(16, 18), opt(16, 19), opt(16, 20), + + opt(17, 1), opt(17, 2), opt(17, 3), opt(17, 4), opt(17, 5), cpt(17, 6), cpt(17, 7), cpt(17, 8), cpt(17, 9), opt(17, 10), + opt(17, 11), opt(17, 12), opt(17, 13), opt(17, 14), opt(17, 15), opt(17, 16), opt(17, 17), opt(17, 18), opt(17, 19), opt(17, 20), + + opt(18, 1), opt(18, 2), opt(18, 3), opt(18, 4), opt(18, 5), opt(18, 6), opt(18, 7), opt(18, 8), opt(18, 9), opt(18, 10), + opt(18, 11), opt(18, 12), opt(18, 13), opt(18, 14), opt(18, 15), opt(18, 16), opt(18, 17), opt(18, 18), opt(18, 19), opt(18, 20), + + opt(19, 1), opt(19, 2), opt(19, 3), opt(19, 4), opt(19, 5), opt(19, 6), opt(19, 7), opt(19, 8), opt(19, 9), opt(19, 10), + opt(19, 11), opt(19, 12), opt(19, 13), opt(19, 14), opt(19, 15), opt(19, 16), opt(19, 17), opt(19, 18), opt(19, 19), opt(19, 20), + + opt(20, 1), opt(20, 2), opt(20, 3), opt(20, 4), opt(20, 5), opt(20, 6), opt(20, 7), opt(20, 8), opt(20, 9), opt(20, 10), + opt(20, 11), opt(20, 12), opt(20, 13), opt(20, 14), opt(20, 15), opt(20, 16), opt(20, 17), opt(20, 18), opt(20, 19), opt(20, 20), + }, + }, + } + for name, tc := range tests { + tc := tc + t.Run(name, func(t *testing.T) { fn(t, tc) }) + } +} + +func TestNewFromPolygon(t *testing.T) { + + type tcase struct { + polygon geom.Polygon + err error + } + fn := func(t *testing.T, tc tcase) { + t.Parallel() + defer func() { + if r := recover(); r != nil { + t.Errorf("panic, expected nil got %v", r) + } + }() + // Don't want hm to be optimized away. + hm, err := NewFromPolygons(nil, tc.polygon.LinearRings()) + if tc.err != nil { + if err == nil { + t.Errorf("error, expected %v got nil", tc.err) + + } + if err.Error() != tc.err.Error() { + t.Errorf("error, expected %v got %v", tc.err.Error(), err.Error()) + } + return + } + if err != nil { + t.Errorf("error, expected nil got %v", err) + return + } + _ = hm + } + tests := map[string]tcase{ + "Nil Polygon": {polygon: nil}, + "Basic Polygon": {polygon: geom.Polygon{}}, + "With one nil": {polygon: geom.Polygon{nil}, err: ErrInvalidLineString}, + "With one empty line": {polygon: geom.Polygon{[][2]float64{}}, err: ErrInvalidLineString}, + "With one one non-empty line": {polygon: geom.Polygon{[][2]float64{{10, 10}, {20, 10}, {20, 20}, {10, 20}}}}, + } + for name, tc := range tests { + tc := tc + t.Run(name, func(t *testing.T) { fn(t, tc) }) + } +} diff --git a/planar/makevalid/hitmap/polygon_hitmap.go b/planar/makevalid/hitmap/polygon_hitmap.go new file mode 100644 index 00000000..9ae8d000 --- /dev/null +++ b/planar/makevalid/hitmap/polygon_hitmap.go @@ -0,0 +1,84 @@ +package hitmap + +import ( + "sort" + + "github.com/go-spatial/geom" +) + +// PolygonHM implements a basic hit map that gives the label for a point based on the order of the rings. +type PolygonHM struct { + // clipBox this is going to be either the clipping area or the bouding box of all the rings. + // This allows us to quickly determine if a point is outside the set of rings. + clipBox *geom.Extent + // These are the rings + rings []*Ring +} + +// NewFromPolygons assumes that the outer ring of each polygon is inside, and each inner ring is inside. +func NewFromPolygons(clipbox *geom.Extent, plys ...[][][2]float64) (*PolygonHM, error) { + + hm := &PolygonHM{ + clipBox: new(geom.Extent), + } + for i := range plys { + if len(plys[i]) == 0 { + continue + } + { + ring, err := NewRing(plys[i][0], Inside) + if err != nil { + return nil, err + } + if clipbox == nil { + // add to the bb of ring to the hm clipbox + hm.clipBox.Add(ring.bbox) + } + hm.rings = append(hm.rings, ring) + } + if len(plys[i]) <= 1 { + continue + } + for j := range plys[i][1:] { + // plys we assume the first ring is inside, and all other rings are outside. + ring, err := NewRing(plys[i][j+1], Outside) + if err != nil { + return nil, err + } + if clipbox == nil { + // add to the bb of ring to the hm clipbox + hm.clipBox.Add(ring.bbox) + } + hm.rings = append(hm.rings, ring) + } + } + sort.Sort(bySmallestBBArea(hm.rings)) + return hm, nil +} + +// LabelFor returns the label for the given point. +func (hm *PolygonHM) LabelFor(pt [2]float64) Label { + // nil clipBox contains all points. + if hm == nil || !hm.clipBox.ContainsPoint(pt) { + return Outside + } + + // TODO(gdey): See if it make sense to change data structures here. + // For now we iterate through all the rings, but maybe an r-tree would make + // sense here, or some "smart" container that would use an r-tree or iterate + // through all the points depending on the number of things. + + // We assume the []*Rings are sorted in from smallest area to largest area. + for i := range hm.rings { + if hm.rings[i].Contains(pt) { + return hm.rings[i].Label + } + } + return Outside +} + +// Extent returns the extent of the hitmap. +func (hm *PolygonHM) Extent() [4]float64 { return hm.clipBox.Extent() } + +// Area returns the area covered by the hitmap. +func (hm *PolygonHM) Area() float64 { return hm.clipBox.Area() } diff --git a/planar/makevalid/hitmap/ring.go b/planar/makevalid/hitmap/ring.go new file mode 100644 index 00000000..b37dd41b --- /dev/null +++ b/planar/makevalid/hitmap/ring.go @@ -0,0 +1,186 @@ +package hitmap + +import ( + "github.com/go-spatial/geom" + "github.com/go-spatial/geom/internal/rtreego" + "github.com/go-spatial/geom/planar" +) + +type segRect struct { + cx float64 + deltax float64 + seg [2][2]float64 + isvert bool + m, b float64 + rect *rtreego.Rect +} + +func (sr *segRect) Bounds() *rtreego.Rect { + if sr == nil { + return nil + } + return sr.rect +} + +func (sr *segRect) y4x(x float64) float64 { return sr.m*x + sr.b } + +const smallep = 0.00001 + +func minmaxbyx(seg [2][2]float64) [2][2]float64 { + if seg[0][0] == seg[1][0] { + if seg[0][1] > seg[1][1] { + return [2][2]float64{seg[1], seg[0]} + } + return seg + } + if seg[0][0] > seg[1][0] { + return [2][2]float64{seg[1], seg[0]} + } + return seg +} + +func newSegRects(segs ...[2][2]float64) (segRects []*segRect) { + for i := range segs { + minx, maxx := segs[i][0][0], segs[i][1][0] + if minx > maxx { + minx, maxx = maxx, minx + } + deltax := maxx - minx + cx := minx + (deltax / 2) + if deltax <= 0 { + // https://github.com/dhconnelly/rtreego/issues/18 + deltax = smallep + } + + m, b, defined := planar.Slope(segs[i]) + + rect, _ := rtreego.NewRect(rtreego.Point{cx}, []float64{deltax}) + segRects = append(segRects, &segRect{ + cx: cx, + deltax: deltax, + // We store the segment with the smaller x point first. + seg: minmaxbyx(segs[i]), + m: m, + b: b, + isvert: !defined, + rect: rect, + }) + } + return segRects + +} + +func createSegments(ls [][2]float64, isClosed bool) (segs [][2][2]float64, err error) { + if len(ls) <= 1 { + return nil, ErrInvalidLineString + + } + i := 0 + for j := 1; j < len(ls); j++ { + segs = append(segs, [2][2]float64{ls[i], ls[j]}) + i = j + } + if isClosed { + segs = append(segs, [2][2]float64{ls[len(ls)-1], ls[0]}) + } + return segs, nil +} + +type Ring struct { + tree *rtreego.Rtree + bbox *geom.Extent + verts map[[2]float64]struct{} + Label Label +} + +// Contains returns weather the point is contained by the ring, if the point is on the border it is considered not contained. +func (r Ring) Contains(pt [2]float64) bool { + if r.tree == nil { + return false + } + if !r.bbox.ContainsPoint(pt) { + return false + } + bb, _ := rtreego.NewRect(rtreego.Point{pt[0]}, []float64{r.bbox.XSpan()}) + results := r.tree.SearchIntersect(bb) + if len(results) == 0 { + return false + } + + var segs []*segRect + for i := range results { + seg, ok := results[i].(*segRect) + // don't know how to deal with this rect, ignore. + if !ok { + continue + } + // Is the x of the point even on the line? + if seg.seg[0][0] > pt[0] || pt[0] > seg.seg[1][0] { + continue + } + if seg.isvert { + if seg.cx == pt[0] && seg.seg[0][1] <= pt[1] && pt[1] <= seg.seg[1][1] { + // We ended up on a boundry, assume it's not contained. + return false + } + // ignore this vertical line. + continue + } + ly := seg.y4x(pt[0]) + + // This may need a different comparison method due to floating point precision issues. + if pt[1] == ly { + // We ended up on a boundry, assume it's not contained. + return false + } + if pt[1] < ly { + continue + } + // Save line to process later. + segs = append(segs, seg) + } + + found := 0 + // seen := map[float64]struct{}{} + for _, seg := range segs { + if _, vok := r.verts[[2]float64{pt[0], seg.y4x(pt[0])}]; vok { + return r.Contains([2]float64{pt[0] + 0.0001, pt[1]}) + } + found++ + } + // if it's odd the point is contained. + return found%2 != 0 +} + +func NewRing(ring [][2]float64, label Label) (*Ring, error) { + bb := geom.NewExtent(ring...) + verts := make(map[[2]float64]struct{}) + for i := range ring { + verts[ring[i]] = struct{}{} + } + segs, err := createSegments(ring, true) + if err != nil { + return nil, err + } + rects := newSegRects(segs...) + var rs = make([]rtreego.Spatial, len(rects)) + for i := range rects { + rs[i] = rects[i] + } + rtree := rtreego.NewTree(1, 2, 5, rs...) + return &Ring{tree: rtree, bbox: bb, verts: verts, Label: label}, nil +} + +type bySmallestBBArea []*Ring + +// Sort Interface + +func (rs bySmallestBBArea) Len() int { return len(rs) } +func (rs bySmallestBBArea) Less(i, j int) bool { + ia, ja := rs[i].bbox.Area(), rs[j].bbox.Area() + if ia == ja { + return rs[i].Label == Outside + } + return ia < ja +} +func (rs bySmallestBBArea) Swap(i, j int) { rs[i], rs[j] = rs[j], rs[i] }