From 6449bd7e0f2bcc0eb1cf7d24f16cd02d1ea87e98 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Wed, 2 May 2018 10:32:44 -0700 Subject: [PATCH 1/6] Basic Planar fns, Simplifier interface Added some handy planar functions, and defined an Interface for Simplifers. This way we can add and change out simplifiers easily. --- encoding/error.go | 14 ----- encoding/geojson/geojson.go | 3 +- encoding/geojson/geojson_test.go | 3 +- encoding/wkb/wkb.go | 10 +--- encoding/wkt/wkt.go | 10 +--- encoding/wkt/wkt_test.go | 2 +- errors.go | 12 +++++ geom.go | 5 -- planar/planar.go | 86 +++++++++++++++++++++++++++++++ planar/planar_test.go | 46 +++++++++++++++++ planar/simplify.go | 87 ++++++++++++++++++++++++++++++++ 11 files changed, 236 insertions(+), 42 deletions(-) create mode 100644 errors.go create mode 100644 planar/planar.go create mode 100644 planar/planar_test.go create mode 100644 planar/simplify.go diff --git a/encoding/error.go b/encoding/error.go index ac2dd82a..523b90ce 100644 --- a/encoding/error.go +++ b/encoding/error.go @@ -1,15 +1 @@ package encoding - -import ( - "fmt" - - "github.com/go-spatial/geom" -) - -type ErrUnknownGeometry struct { - Geom geom.Geometry -} - -func (e ErrUnknownGeometry) Error() string { - return fmt.Sprintf("unknown geometry: %T", e.Geom) -} diff --git a/encoding/geojson/geojson.go b/encoding/geojson/geojson.go index 19b48fd7..7ad80823 100644 --- a/encoding/geojson/geojson.go +++ b/encoding/geojson/geojson.go @@ -4,7 +4,6 @@ import ( "encoding/json" "github.com/go-spatial/geom" - "github.com/go-spatial/geom/encoding" ) type GeoJSONType string @@ -95,7 +94,7 @@ func (geo Geometry) MarshalJSON() ([]byte, error) { }) default: - return nil, encoding.ErrUnknownGeometry{g} + return nil, geom.ErrUnknownGeometry{g} } } diff --git a/encoding/geojson/geojson_test.go b/encoding/geojson/geojson_test.go index 41314381..9da9cfe4 100644 --- a/encoding/geojson/geojson_test.go +++ b/encoding/geojson/geojson_test.go @@ -6,7 +6,6 @@ import ( "testing" "github.com/go-spatial/geom" - "github.com/go-spatial/geom/encoding" "github.com/go-spatial/geom/encoding/geojson" ) @@ -109,7 +108,7 @@ func TestFeatureMarshalJSON(t *testing.T) { geom: nil, expectedErr: json.MarshalerError{ Type: reflect.TypeOf(geojson.Geometry{}), - Err: encoding.ErrUnknownGeometry{nil}, + Err: geom.ErrUnknownGeometry{nil}, }, }, } diff --git a/encoding/wkb/wkb.go b/encoding/wkb/wkb.go index d6bf9c59..4298113d 100644 --- a/encoding/wkb/wkb.go +++ b/encoding/wkb/wkb.go @@ -17,14 +17,6 @@ import ( "github.com/go-spatial/geom/encoding/wkb/internal/encode" ) -type ErrUnknownGeometry struct { - Geom geom.Geometry -} - -func (e ErrUnknownGeometry) Error() string { - return fmt.Sprintf("Unknown Geometry! %v", e.Geom) -} - type ErrUnknownGeometryType struct { Typ uint32 } @@ -108,7 +100,7 @@ func _encode(en *encode.Encoder, g geom.Geometry) error { } } default: - return ErrUnknownGeometry{g} + return geom.ErrUnknownGeometry{g} } return en.Err() } diff --git a/encoding/wkt/wkt.go b/encoding/wkt/wkt.go index 3bb00b73..aec22d4b 100644 --- a/encoding/wkt/wkt.go +++ b/encoding/wkt/wkt.go @@ -104,14 +104,6 @@ use to take a tagola.Geometry and convert it to a wkt string. It will, also, contain functions to parse a wkt string into a wkb.Geometry. */ -type ErrUnknownGeometry struct { - Geom geom.Geometry -} - -func (e ErrUnknownGeometry) Error() string { - return fmt.Sprintf("Unknown Geometry! %v", e.Geom) -} - func _encode(geo geom.Geometry) string { switch g := geo.(type) { @@ -173,7 +165,7 @@ func _encode(geo geom.Geometry) string { func Encode(geo geom.Geometry) (string, error) { switch g := geo.(type) { default: - return "", ErrUnknownGeometry{geo} + return "", geom.ErrUnknownGeometry{geo} case geom.Pointer: // POINT( 10 10) if isNil(g) { diff --git a/encoding/wkt/wkt_test.go b/encoding/wkt/wkt_test.go index d5ab3d59..f3df2d8b 100644 --- a/encoding/wkt/wkt_test.go +++ b/encoding/wkt/wkt_test.go @@ -33,7 +33,7 @@ func TestEncode(t *testing.T) { tests := map[string]map[string]tcase{ "Point": { "empty nil": { - Err: ErrUnknownGeometry{nil}, + Err: geom.ErrUnknownGeometry{nil}, }, "empty": { Geom: (*geom.Point)(nil), diff --git a/errors.go b/errors.go new file mode 100644 index 00000000..7fead4be --- /dev/null +++ b/errors.go @@ -0,0 +1,12 @@ +package geom + +import "fmt" + +// ErrUnknownGeometry represents an objects that is not a known geom geometry. +type ErrUnknownGeometry struct { + Geom Geometry +} + +func (e ErrUnknownGeometry) Error() string { + return fmt.Sprintf("unknown geometry: %T", e.Geom) +} diff --git a/geom.go b/geom.go index 17a7c6e1..0045119a 100644 --- a/geom.go +++ b/geom.go @@ -1,11 +1,6 @@ // Package geom describes geometry interfaces. package geom -import "errors" - -// ErrUnknownGeometry is returned when the geometry type is unknown or unsupported. -var ErrUnknownGeometry = errors.New("unknown geometry") - // Geometry is an object with a spatial reference. // if a method accepts a Geometry type it's only expected to support the geom types in this package type Geometry interface{} diff --git a/planar/planar.go b/planar/planar.go new file mode 100644 index 00000000..c6b20fbe --- /dev/null +++ b/planar/planar.go @@ -0,0 +1,86 @@ +package planar + +import ( + "math" + "math/big" +) + +const Rad = math.Pi / 180 + +type PointLineDistanceFunc func(line [2][2]float64, point [2]float64) float64 + +// PerpendicularDistance provides the distance between a line and a point in Euclidean space. +func PerpendicularDistance(line [2][2]float64, point [2]float64) float64 { + + deltaX := line[1][0] - line[0][0] + deltaY := line[1][1] - line[0][1] + denom := math.Abs((deltaY * point[0]) - (deltaX * point[1]) + (line[1][0] * line[0][1]) - (line[1][1] * line[0][0])) + num := math.Sqrt(math.Pow(deltaY, 2) + math.Pow(deltaX, 2)) + if num == 0 { + return 0 + } + return denom / num +} + +type Haversine struct { + Semimajor *float64 +} + +// PerpendicularDistance returns the distance between a point and a line in meters, using the Harversine distance. +func (hs Haversine) PerpendicularDistance(line [2][2]float64, point [2]float64) float64 { + dist := hs.Distance(line[0], point) + angle := Angle3Pts(line[0], line[1], point) + return math.Abs(dist * math.Sin(angle)) +} + +// Distance returns the Haversine distance between two points in meters +func (hs Haversine) Distance(pt1, pt2 [2]float64) float64 { + rpt1x, rpt2x := pt1[0]*Rad, pt2[0]*Rad + distx := rpt1x - rpt2x + disty := (pt1[1] * Rad) - (pt2[1] * Rad) + dist := 2 * math.Asin( + math.Sqrt( + + (math.Pow(math.Sin(distx/2), 2)+math.Cos(rpt1x))* + math.Cos(rpt2x)* + (math.Pow(math.Sin(disty/2), 2)), + ), + ) + if math.IsNaN(dist) { + dist = 0 + } + + // https://resources.arcgis.com/en/help/main/10.1/index.html#//003r00000003000000 + // we will default to the semimajor of WGS84 + semimajor := big.NewFloat(6378137) + if hs.Semimajor != nil { + semimajor = big.NewFloat(*hs.Semimajor) + } + + // 4 digits of precision + d, _ := new(big.Float).SetPrec(16).Mul(big.NewFloat(dist), semimajor).Float64() + return d +} + +// Slope — finds the Slope of a line +func Slope(line [2][2]float64) (m, b float64, defined bool) { + dx := line[1][0] - line[0][0] + dy := line[1][1] - line[0][1] + if dx == 0 || dy == 0 { + // if dx == 0 then m == 0; and the intercept is y. + // However if the lines are verticle then the slope is not defined. + return 0, line[0][1], dx != 0 + } + m = dy / dx + b = line[0][1] - (m * line[0][0]) + return m, b, true +} + +// Angle3Pts returns the angle between three points… +func Angle3Pts(pt1, pt2, pt3 [2]float64) float64 { + m1, _, _ := Slope([2][2]float64{pt2, pt1}) + m2, _, _ := Slope([2][2]float64{pt3, pt1}) + // 6 digits of prec + d, _ := new(big.Float).SetPrec(6*4).Sub(big.NewFloat(m2), big.NewFloat(m1)).Float64() + return d +} diff --git a/planar/planar_test.go b/planar/planar_test.go new file mode 100644 index 00000000..296aa204 --- /dev/null +++ b/planar/planar_test.go @@ -0,0 +1,46 @@ +package planar + +import ( + "strconv" + "testing" +) + +func TestSlope(t *testing.T) { + type tcase struct { + line [2][2]float64 + m, b float64 + defined bool + } + + fn := func(t *testing.T, tc tcase) { + t.Parallel() + gm, gb, gd := Slope(tc.line) + if tc.defined != gd { + t.Errorf("sloped defined, expected %v got %v", tc.defined, gd) + return + } + // if the slope is not defined, line is verticle and m,b don't have good values. + if !tc.defined { + return + } + if tc.m != gm { + t.Errorf("sloped, expected %v got %v", tc.m, gm) + + } + if tc.b != gb { + t.Errorf("sloped intercept, expected %v got %v", tc.b, gb) + } + } + tests := []tcase{ + { + line: [2][2]float64{{0, 0}, {10, 10}}, + m: 1, + b: 0, + defined: true, + }, + } + for i := range tests { + tc := tests[i] + t.Run(strconv.FormatInt(int64(i), 10), func(t *testing.T) { fn(t, tc) }) + } +} diff --git a/planar/simplify.go b/planar/simplify.go new file mode 100644 index 00000000..f72f57d8 --- /dev/null +++ b/planar/simplify.go @@ -0,0 +1,87 @@ +package planar + +import "github.com/go-spatial/geom" + +// Simplifer is an interface for Simplifying geometries. +type Simplifer interface { + Simplify(linestring [][2]float64, isClosed bool) ([][2]float64, error) +} + +func simplifyPolygon(simplifer Simplifer, plg [][][2]float64, isClosed bool) (ret [][][2]float64, err error) { + ret = make([][][2]float64, len(plg)) + for i := range plg { + ls, err := simplifer.Simplify(plg[i], isClosed) + if err != nil { + return nil, err + } + ret[i] = ls + } + return ret, nil + +} + +// Simplify will simplify the provided geometry using the provided simplifer. +// If the simplifer is nil, no simplification will be attempted. +func Simplify(simplifer Simplifer, geometry geom.Geometry) (geom.Geometry, error) { + + if simplifer == nil { + return geometry, nil + } + + switch gg := geometry.(type) { + + case geom.Collectioner: + + geos := gg.Geometries() + coll := make([]geom.Geometry, len(geos)) + for i := range geos { + geo, err := Simplify(simplifer, geos[i]) + if err != nil { + return nil, err + } + coll[i] = geo + } + return geom.Collection(coll), nil + + case geom.MultiPolygoner: + + plys := gg.Polygons() + mply := make([][][][2]float64, len(plys)) + for i := range plys { + ply, err := simplifyPolygon(simplifer, plys[i], true) + if err != nil { + return nil, err + } + mply[i] = ply + } + return geom.MultiPolygon(mply), nil + + case geom.Polygoner: + + ply, err := simplifyPolygon(simplifer, gg.LinearRings(), true) + if err != nil { + return nil, err + } + return geom.Polygon(ply), nil + + case geom.MultiLineStringer: + + mls, err := simplifyPolygon(simplifer, gg.LineStrings(), false) + if err != nil { + return nil, err + } + return geom.MultiLineString(mls), nil + + case geom.LineStringer: + + ls, err := simplifer.Simplify(gg.Verticies(), false) + if err != nil { + return nil, err + } + return geom.LineString(ls), nil + + default: // Points, MutliPoints or anything else. + return geometry, nil + + } +} From 5ae79529e89f4d4ae2949d8f6ffaee0671b54cf1 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Wed, 2 May 2018 11:50:40 -0700 Subject: [PATCH 2/6] Added doublaspeucker simplification algo. --- planar/simplify/douglaspeucker.go | 49 +++++++++++++++++++++ planar/simplify/douglaspeucker_test.go | 61 ++++++++++++++++++++++++++ 2 files changed, 110 insertions(+) create mode 100644 planar/simplify/douglaspeucker.go create mode 100644 planar/simplify/douglaspeucker_test.go diff --git a/planar/simplify/douglaspeucker.go b/planar/simplify/douglaspeucker.go new file mode 100644 index 00000000..ca071040 --- /dev/null +++ b/planar/simplify/douglaspeucker.go @@ -0,0 +1,49 @@ +package simplify + +import "github.com/go-spatial/geom/planar" + +type DouglasPeucker struct { + + // Tolerance is the tolerance used to eliminate points, a tolerance of zero is not eliminate any points. + Tolerance float64 + + // Dist is the distance function to use, defaults to planar.PerpendicularDistance + Dist planar.PointLineDistanceFunc +} + +func (dp DouglasPeucker) Simplify(linestring [][2]float64, isClosed bool) ([][2]float64, error) { + + if dp.Tolerance <= 0 || len(linestring) <= 2 { + return linestring, nil + } + + dmax, idx := 0.0, 0 + dist := planar.PerpendicularDistance + if dp.Dist != nil { + dist = dp.Dist + } + + line := [2][2]float64{linestring[0], linestring[len(linestring)-1]} + + // Find the point that is the furthest away. + for i := 1; i <= len(linestring)-2; i++ { + d := dist(line, linestring[i]) + if d > dmax { + dmax, idx = d, i + } + } + + // If the furtherest point is greater then tolerance, we split at that point, and look again at each + // subsections. + if dmax > dp.Tolerance { + if len(linestring) <= 3 { + return linestring, nil + } + rec1, _ := dp.Simplify(linestring[0:idx], isClosed) + rec2, _ := dp.Simplify(linestring[idx:], isClosed) + return append(rec1, rec2...), nil + } + + // Drop all points between the end points. + return line[:], nil +} diff --git a/planar/simplify/douglaspeucker_test.go b/planar/simplify/douglaspeucker_test.go new file mode 100644 index 00000000..7748e876 --- /dev/null +++ b/planar/simplify/douglaspeucker_test.go @@ -0,0 +1,61 @@ +package simplify + +import ( + "reflect" + "testing" + + "github.com/go-spatial/geom/planar" +) + +func TestDouglasPeucker(t *testing.T) { + type tcase struct { + l [][2]float64 + dp DouglasPeucker + el [][2]float64 + } + + fn := func(t *testing.T, tc tcase) { + gl, err := tc.dp.Simplify(tc.l, false) + // Douglas Peucker should never return an error. + // This is more of a sanity check. + if err != nil { + t.Errorf("Douglas Peucker error, expected nil got %v", err) + return + } + if !reflect.DeepEqual(tc.el, gl) { + t.Errorf("simplified points, expected %v got %v", tc.el, gl) + return + } + // Let's try it with true, it should not matter, as DP does not care. + // More sanity checking. + gl, _ = tc.dp.Simplify(tc.l, true) + + if !reflect.DeepEqual(tc.el, gl) { + t.Errorf("simplified points (true), expected %v got %v", tc.el, gl) + return + } + } + + tests := map[string]tcase{ + "simple box": { + l: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + dp: DouglasPeucker{ + Tolerance: 0.001, + }, + el: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + }, + "haversine simple box": { + l: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + dp: DouglasPeucker{ + Tolerance: 0.001, + Dist: planar.Haversine{}.PerpendicularDistance, + }, + el: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + }, + } + + for name, tc := range tests { + tc := tc + t.Run(name, func(t *testing.T) { fn(t, tc) }) + } +} From 14233ae61f9be9fdc1b53b754b009c431d791d29 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Wed, 2 May 2018 13:22:44 -0700 Subject: [PATCH 3/6] Added reference to the formula for PerpendicularDisance. --- planar/planar.go | 2 ++ 1 file changed, 2 insertions(+) diff --git a/planar/planar.go b/planar/planar.go index c6b20fbe..1893ba7c 100644 --- a/planar/planar.go +++ b/planar/planar.go @@ -10,10 +10,12 @@ const Rad = math.Pi / 180 type PointLineDistanceFunc func(line [2][2]float64, point [2]float64) float64 // PerpendicularDistance provides the distance between a line and a point in Euclidean space. +// ref: https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points func PerpendicularDistance(line [2][2]float64, point [2]float64) float64 { deltaX := line[1][0] - line[0][0] deltaY := line[1][1] - line[0][1] + denom := math.Abs((deltaY * point[0]) - (deltaX * point[1]) + (line[1][0] * line[0][1]) - (line[1][1] * line[0][0])) num := math.Sqrt(math.Pow(deltaY, 2) + math.Pow(deltaX, 2)) if num == 0 { From 0da3a097dd1c0d79cedc964f8101fe0982b817b0 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Thu, 3 May 2018 13:38:50 -0700 Subject: [PATCH 4/6] Updated according to reviews. Moved spherical functions to the spherical namespace. Added references to docs used for algos. --- planar/planar.go | 60 ++++--------------------------------- planar/planar_test.go | 4 +++ spherical/spherical.go | 58 +++++++++++++++++++++++++++++++++++ spherical/spherical_test.go | 19 ++++++++++++ 4 files changed, 87 insertions(+), 54 deletions(-) create mode 100644 spherical/spherical.go create mode 100644 spherical/spherical_test.go diff --git a/planar/planar.go b/planar/planar.go index 1893ba7c..b48b346e 100644 --- a/planar/planar.go +++ b/planar/planar.go @@ -2,7 +2,6 @@ package planar import ( "math" - "math/big" ) const Rad = math.Pi / 180 @@ -15,53 +14,15 @@ func PerpendicularDistance(line [2][2]float64, point [2]float64) float64 { deltaX := line[1][0] - line[0][0] deltaY := line[1][1] - line[0][1] + deltaXSq := deltaX * deltaX + deltaYSq := deltaY * deltaY - denom := math.Abs((deltaY * point[0]) - (deltaX * point[1]) + (line[1][0] * line[0][1]) - (line[1][1] * line[0][0])) - num := math.Sqrt(math.Pow(deltaY, 2) + math.Pow(deltaX, 2)) - if num == 0 { + num := math.Abs((deltaY * point[0]) - (deltaX * point[1]) + (line[1][0] * line[0][1]) - (line[1][1] * line[0][0])) + denom := math.Sqrt(deltaYSq + deltaXSq) + if denom == 0 { return 0 } - return denom / num -} - -type Haversine struct { - Semimajor *float64 -} - -// PerpendicularDistance returns the distance between a point and a line in meters, using the Harversine distance. -func (hs Haversine) PerpendicularDistance(line [2][2]float64, point [2]float64) float64 { - dist := hs.Distance(line[0], point) - angle := Angle3Pts(line[0], line[1], point) - return math.Abs(dist * math.Sin(angle)) -} - -// Distance returns the Haversine distance between two points in meters -func (hs Haversine) Distance(pt1, pt2 [2]float64) float64 { - rpt1x, rpt2x := pt1[0]*Rad, pt2[0]*Rad - distx := rpt1x - rpt2x - disty := (pt1[1] * Rad) - (pt2[1] * Rad) - dist := 2 * math.Asin( - math.Sqrt( - - (math.Pow(math.Sin(distx/2), 2)+math.Cos(rpt1x))* - math.Cos(rpt2x)* - (math.Pow(math.Sin(disty/2), 2)), - ), - ) - if math.IsNaN(dist) { - dist = 0 - } - - // https://resources.arcgis.com/en/help/main/10.1/index.html#//003r00000003000000 - // we will default to the semimajor of WGS84 - semimajor := big.NewFloat(6378137) - if hs.Semimajor != nil { - semimajor = big.NewFloat(*hs.Semimajor) - } - - // 4 digits of precision - d, _ := new(big.Float).SetPrec(16).Mul(big.NewFloat(dist), semimajor).Float64() - return d + return num / denom } // Slope — finds the Slope of a line @@ -77,12 +38,3 @@ func Slope(line [2][2]float64) (m, b float64, defined bool) { b = line[0][1] - (m * line[0][0]) return m, b, true } - -// Angle3Pts returns the angle between three points… -func Angle3Pts(pt1, pt2, pt3 [2]float64) float64 { - m1, _, _ := Slope([2][2]float64{pt2, pt1}) - m2, _, _ := Slope([2][2]float64{pt3, pt1}) - // 6 digits of prec - d, _ := new(big.Float).SetPrec(6*4).Sub(big.NewFloat(m2), big.NewFloat(m1)).Float64() - return d -} diff --git a/planar/planar_test.go b/planar/planar_test.go index 296aa204..f050d0ac 100644 --- a/planar/planar_test.go +++ b/planar/planar_test.go @@ -38,6 +38,10 @@ func TestSlope(t *testing.T) { b: 0, defined: true, }, + { + line: [2][2]float64{{1, 7}, {1, 17}}, + defined: false, + }, } for i := range tests { tc := tests[i] diff --git a/spherical/spherical.go b/spherical/spherical.go new file mode 100644 index 00000000..f4902db5 --- /dev/null +++ b/spherical/spherical.go @@ -0,0 +1,58 @@ +package spherical + +import ( + "math" +) + +const Rad = math.Pi / 180 + +// https://resources.arcgis.com/en/help/main/10.1/index.html#//003r00000003000000 +// we will default to the semimajor of WGS84 +const wgs84Semimajor = 6378137.0 + +// Angle3Pts returns the angle between three points… +func Angle3Pts(pt1, pt2, pt3 [2]float64) float64 { + m1 := math.Atan2((pt1[1] - pt2[1]), (pt1[0] - pt2[0])) + m2 := math.Atan2((pt1[1] - pt3[1]), (pt1[0] - pt3[0])) + return m2 - m1 +} + +type Haversine struct { + Semimajor *float64 +} + +// semimajor returns a default value if the main one is not set. +func (hs Haversine) semimajor() float64 { + if hs.Semimajor == nil { + return wgs84Semimajor + } + return *hs.Semimajor +} + +// Distance returns the Haversine distance between two points in meters +// ref: https://en.wikipedia.org/wiki/Haversine_formula +func (hs Haversine) Distance(pt1, pt2 [2]float64) float64 { + rpt1x, rpt2x := pt1[0]*Rad, pt2[0]*Rad + distx := rpt1x - rpt2x + disty := (pt1[1] * Rad) - (pt2[1] * Rad) + dist := 2 * math.Asin( + math.Sqrt( + (math.Pow(math.Sin(distx/2), 2)+math.Cos(rpt1x))* + math.Cos(rpt2x)* + (math.Pow(math.Sin(disty/2), 2)), + ), + ) + if math.IsNaN(dist) { + dist = 0 + } + + return dist * hs.semimajor() +} + +// PerpendicularDistance returns the distance between a point and a line in meters, using the Harversine distance. +// cross track distance: https://www.movable-type.co.uk/scripts/latlong.html +func (hs Haversine) PerpendicularDistance(line [2][2]float64, point [2]float64) float64 { + dist := hs.Distance(line[0], point) + angle := Angle3Pts(line[0], line[1], point) + return math.Abs(dist * math.Sin(angle)) +} diff --git a/spherical/spherical_test.go b/spherical/spherical_test.go new file mode 100644 index 00000000..80166a3e --- /dev/null +++ b/spherical/spherical_test.go @@ -0,0 +1,19 @@ +package spherical + +import "testing" + +// TODO(gdey): Need to add real tests. + +func TestSemimajor(t *testing.T) { + hs := Haversine{} + sm := hs.semimajor() + if sm != wgs84Semimajor { + t.Errorf("seimimajor, expected %v got %v", wgs84Semimajor, sm) + } + f := 1.0 + hs.Semimajor = &f + sm = hs.semimajor() + if sm != 1.0 { + t.Errorf("seimimajor, expected %v got %v", 1.0, sm) + } +} From 0fdf75481abc68c0a652b94d145bec4fae658124 Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Thu, 3 May 2018 15:44:44 -0700 Subject: [PATCH 5/6] Removed haversine test for simplify. Because not sure if it is working. --- planar/simplify/douglaspeucker_test.go | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/planar/simplify/douglaspeucker_test.go b/planar/simplify/douglaspeucker_test.go index 7748e876..3a1842bb 100644 --- a/planar/simplify/douglaspeucker_test.go +++ b/planar/simplify/douglaspeucker_test.go @@ -3,8 +3,6 @@ package simplify import ( "reflect" "testing" - - "github.com/go-spatial/geom/planar" ) func TestDouglasPeucker(t *testing.T) { @@ -44,14 +42,19 @@ func TestDouglasPeucker(t *testing.T) { }, el: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, }, - "haversine simple box": { - l: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, - dp: DouglasPeucker{ - Tolerance: 0.001, - Dist: planar.Haversine{}.PerpendicularDistance, + /* + This isn't working: I don't if this is a problem or not; we will not be using Haversine + for things + + "haversine simple box": { + l: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, + dp: DouglasPeucker{ + Tolerance: 0.001, + Dist: spherical.Haversine{}.PerpendicularDistance, + }, + el: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, }, - el: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, - }, + */ } for name, tc := range tests { From abb6d65a4109e5b87406cef9e232f7004ba2916d Mon Sep 17 00:00:00 2001 From: Gautam Dey Date: Tue, 29 May 2018 10:16:02 -0700 Subject: [PATCH 6/6] [issue #7] Pull out Douglaspeucker Algo to the geom package. Removed the spherical version as we don't need it and there may be an issue with the distance formula that was being used, and got pulled. Improved test scaffolding, with the aim to add more tests. --- planar/simplify/debug.go | 16 +++++++ planar/simplify/douglaspeucker.go | 54 ++++++++++++++++++++++-- planar/simplify/douglaspeucker_test.go | 25 ++++++----- spherical/spherical.go | 58 -------------------------- spherical/spherical_test.go | 19 --------- 5 files changed, 79 insertions(+), 93 deletions(-) create mode 100644 planar/simplify/debug.go delete mode 100644 spherical/spherical.go delete mode 100644 spherical/spherical_test.go diff --git a/planar/simplify/debug.go b/planar/simplify/debug.go new file mode 100644 index 00000000..68b56fa5 --- /dev/null +++ b/planar/simplify/debug.go @@ -0,0 +1,16 @@ +package simplify + +import ( + "log" + "os" +) + +const debug = false + +var logger *log.Logger + +func init() { + if debug { + logger = log.New(os.Stderr, "simplify:", log.Lshortfile|log.LstdFlags) + } +} diff --git a/planar/simplify/douglaspeucker.go b/planar/simplify/douglaspeucker.go index ca071040..abe78ef8 100644 --- a/planar/simplify/douglaspeucker.go +++ b/planar/simplify/douglaspeucker.go @@ -1,6 +1,10 @@ package simplify -import "github.com/go-spatial/geom/planar" +import ( + "strings" + + "github.com/go-spatial/geom/planar" +) type DouglasPeucker struct { @@ -12,11 +16,38 @@ type DouglasPeucker struct { } func (dp DouglasPeucker) Simplify(linestring [][2]float64, isClosed bool) ([][2]float64, error) { + return dp.simplify(0, linestring, isClosed) +} + +func (dp DouglasPeucker) simplify(depth uint8, linestring [][2]float64, isClosed bool) ([][2]float64, error) { + + // helper function for debugging and tracing the code + var printf = func(msg string, depth uint8, params ...interface{}) { + if debug { + ps := make([]interface{}, 1, len(params)+1) + ps[0] = depth + ps = append(ps, params...) + logger.Printf(strings.Repeat(" ", int(depth*2))+"[%v]"+msg, ps...) + } + } if dp.Tolerance <= 0 || len(linestring) <= 2 { + if debug { + if dp.Tolerance <= 0 { + printf("skipping due to Tolerance (%v) ≤ zero:", depth, dp.Tolerance) + + } + if len(linestring) <= 2 { + printf("skipping due to len(linestring) (%v) ≤ two:", depth, len(linestring)) + } + } return linestring, nil } + if debug { + printf("starting linestring: %v ; tolerance: %v", depth, linestring, dp.Tolerance) + } + dmax, idx := 0.0, 0 dist := planar.PerpendicularDistance if dp.Dist != nil { @@ -25,25 +56,42 @@ func (dp DouglasPeucker) Simplify(linestring [][2]float64, isClosed bool) ([][2] line := [2][2]float64{linestring[0], linestring[len(linestring)-1]} + if debug { + printf("starting dmax: %v ; idx %v ; line : %v", depth, dmax, idx, line) + } + // Find the point that is the furthest away. for i := 1; i <= len(linestring)-2; i++ { d := dist(line, linestring[i]) if d > dmax { dmax, idx = d, i } + + if debug { + printf("looking at %v ; d : %v dmax %v ", depth, i, d, dmax) + } } // If the furtherest point is greater then tolerance, we split at that point, and look again at each // subsections. if dmax > dp.Tolerance { if len(linestring) <= 3 { + if debug { + printf("returning linestring %v", depth, linestring) + } return linestring, nil } - rec1, _ := dp.Simplify(linestring[0:idx], isClosed) - rec2, _ := dp.Simplify(linestring[idx:], isClosed) + rec1, _ := dp.simplify(depth+1, linestring[0:idx], isClosed) + rec2, _ := dp.simplify(depth+1, linestring[idx:], isClosed) + if debug { + printf("returning combined lines: %v %v", depth, rec1, rec2) + } return append(rec1, rec2...), nil } // Drop all points between the end points. + if debug { + printf("dropping all points between the end points: %v", depth, line) + } return line[:], nil } diff --git a/planar/simplify/douglaspeucker_test.go b/planar/simplify/douglaspeucker_test.go index 3a1842bb..e3bd4d3e 100644 --- a/planar/simplify/douglaspeucker_test.go +++ b/planar/simplify/douglaspeucker_test.go @@ -1,10 +1,17 @@ package simplify import ( + "flag" "reflect" "testing" ) +var ignoreSanityCheck bool + +func init() { + flag.BoolVar(&ignoreSanityCheck, "ignoreSanityCheck", false, "ignore sanity checks in test cases.") +} + func TestDouglasPeucker(t *testing.T) { type tcase struct { l [][2]float64 @@ -24,6 +31,11 @@ func TestDouglasPeucker(t *testing.T) { t.Errorf("simplified points, expected %v got %v", tc.el, gl) return } + + if ignoreSanityCheck { + return + } + // Let's try it with true, it should not matter, as DP does not care. // More sanity checking. gl, _ = tc.dp.Simplify(tc.l, true) @@ -42,19 +54,6 @@ func TestDouglasPeucker(t *testing.T) { }, el: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, }, - /* - This isn't working: I don't if this is a problem or not; we will not be using Haversine - for things - - "haversine simple box": { - l: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, - dp: DouglasPeucker{ - Tolerance: 0.001, - Dist: spherical.Haversine{}.PerpendicularDistance, - }, - el: [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}}, - }, - */ } for name, tc := range tests { diff --git a/spherical/spherical.go b/spherical/spherical.go deleted file mode 100644 index f4902db5..00000000 --- a/spherical/spherical.go +++ /dev/null @@ -1,58 +0,0 @@ -package spherical - -import ( - "math" -) - -const Rad = math.Pi / 180 - -// https://resources.arcgis.com/en/help/main/10.1/index.html#//003r00000003000000 -// we will default to the semimajor of WGS84 -const wgs84Semimajor = 6378137.0 - -// Angle3Pts returns the angle between three points… -func Angle3Pts(pt1, pt2, pt3 [2]float64) float64 { - m1 := math.Atan2((pt1[1] - pt2[1]), (pt1[0] - pt2[0])) - m2 := math.Atan2((pt1[1] - pt3[1]), (pt1[0] - pt3[0])) - return m2 - m1 -} - -type Haversine struct { - Semimajor *float64 -} - -// semimajor returns a default value if the main one is not set. -func (hs Haversine) semimajor() float64 { - if hs.Semimajor == nil { - return wgs84Semimajor - } - return *hs.Semimajor -} - -// Distance returns the Haversine distance between two points in meters -// ref: https://en.wikipedia.org/wiki/Haversine_formula -func (hs Haversine) Distance(pt1, pt2 [2]float64) float64 { - rpt1x, rpt2x := pt1[0]*Rad, pt2[0]*Rad - distx := rpt1x - rpt2x - disty := (pt1[1] * Rad) - (pt2[1] * Rad) - dist := 2 * math.Asin( - math.Sqrt( - (math.Pow(math.Sin(distx/2), 2)+math.Cos(rpt1x))* - math.Cos(rpt2x)* - (math.Pow(math.Sin(disty/2), 2)), - ), - ) - if math.IsNaN(dist) { - dist = 0 - } - - return dist * hs.semimajor() -} - -// PerpendicularDistance returns the distance between a point and a line in meters, using the Harversine distance. -// cross track distance: https://www.movable-type.co.uk/scripts/latlong.html -func (hs Haversine) PerpendicularDistance(line [2][2]float64, point [2]float64) float64 { - dist := hs.Distance(line[0], point) - angle := Angle3Pts(line[0], line[1], point) - return math.Abs(dist * math.Sin(angle)) -} diff --git a/spherical/spherical_test.go b/spherical/spherical_test.go deleted file mode 100644 index 80166a3e..00000000 --- a/spherical/spherical_test.go +++ /dev/null @@ -1,19 +0,0 @@ -package spherical - -import "testing" - -// TODO(gdey): Need to add real tests. - -func TestSemimajor(t *testing.T) { - hs := Haversine{} - sm := hs.semimajor() - if sm != wgs84Semimajor { - t.Errorf("seimimajor, expected %v got %v", wgs84Semimajor, sm) - } - f := 1.0 - hs.Semimajor = &f - sm = hs.semimajor() - if sm != 1.0 { - t.Errorf("seimimajor, expected %v got %v", 1.0, sm) - } -}