-
Notifications
You must be signed in to change notification settings - Fork 3.8k
/
segmentize.go
112 lines (104 loc) · 4.83 KB
/
segmentize.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
// Copyright 2020 The Cockroach Authors.
//
// Use of this software is governed by the Business Source License
// included in the file licenses/BSL.txt.
//
// As of the Change Date specified in that file, in accordance with
// the Business Source License, use of this software will be governed
// by the Apache License, Version 2.0, included in the file
// licenses/APL.txt.
package geogfn
import (
"math"
"github.com/cockroachdb/cockroach/pkg/geo"
"github.com/cockroachdb/cockroach/pkg/geo/geosegmentize"
"github.com/cockroachdb/cockroach/pkg/sql/pgwire/pgcode"
"github.com/cockroachdb/cockroach/pkg/sql/pgwire/pgerror"
"github.com/golang/geo/s2"
"github.com/twpayne/go-geom"
)
// Segmentize return modified Geography having no segment longer
// that given maximum segment length.
// This works by dividing each segment by a power of 2 to find the
// smallest power less than or equal to the segmentMaxLength.
func Segmentize(geography geo.Geography, segmentMaxLength float64) (geo.Geography, error) {
if math.IsNaN(segmentMaxLength) || math.IsInf(segmentMaxLength, 1 /* sign */) {
return geography, nil
}
geometry, err := geography.AsGeomT()
if err != nil {
return geo.Geography{}, err
}
switch geometry := geometry.(type) {
case *geom.Point, *geom.MultiPoint:
return geography, nil
default:
if segmentMaxLength <= 0 {
return geo.Geography{}, pgerror.Newf(pgcode.InvalidParameterValue, "maximum segment length must be positive")
}
spheroid, err := spheroidFromGeography(geography)
if err != nil {
return geo.Geography{}, err
}
// Convert segmentMaxLength to Angle with respect to earth sphere as
// further calculation is done considering segmentMaxLength as Angle.
segmentMaxAngle := segmentMaxLength / spheroid.SphereRadius()
ret, err := geosegmentize.Segmentize(geometry, segmentMaxAngle, segmentizeCoords)
if err != nil {
return geo.Geography{}, err
}
return geo.MakeGeographyFromGeomT(ret)
}
}
// segmentizeCoords inserts multiple points between given two coordinates and
// return resultant point as flat []float64. Such that distance between any two
// points is less than given maximum segment's length, the total number of
// segments is the power of 2, and all the segments are of the same length.
// Note: List of points does not consist of end point.
func segmentizeCoords(a geom.Coord, b geom.Coord, segmentMaxAngle float64) ([]float64, error) {
if len(a) != len(b) {
return nil, pgerror.Newf(pgcode.InvalidParameterValue, "cannot segmentize two coordinates of different dimensions")
}
if segmentMaxAngle <= 0 {
return nil, pgerror.Newf(pgcode.InvalidParameterValue, "maximum segment angle must be positive")
}
// Converted geom.Coord into s2.Point so we can segmentize the coordinates.
pointA := s2.PointFromLatLng(s2.LatLngFromDegrees(a.Y(), a.X()))
pointB := s2.PointFromLatLng(s2.LatLngFromDegrees(b.Y(), b.X()))
chordAngleBetweenPoints := s2.ChordAngleBetweenPoints(pointA, pointB).Angle().Radians()
// PostGIS' behavior appears to involve cutting this down into segments divisible
// by a power of two. As such, we do not use ceil(chordAngleBetweenPoints/segmentMaxAngle).
//
// This calculation determines the smallest power of 2 that is greater
// than ceil(chordAngleBetweenPoints/segmentMaxAngle).
//
// We can write that power as 2^n in the following inequality
// 2^n >= ceil(chordAngleBetweenPoints/segmentMaxLength) > 2^(n-1).
// We can drop the ceil since 2^n must be an int
// 2^n >= chordAngleBetweenPoints/segmentMaxLength > 2^(n-1).
// Then n = ceil(log2(chordAngleBetweenPoints/segmentMaxLength)).
// Hence numberOfSegmentsToCreate = 2^(ceil(log2(chordAngleBetweenPoints/segmentMaxLength))).
doubleNumberOfSegmentsToCreate := math.Pow(2, math.Ceil(math.Log2(chordAngleBetweenPoints/segmentMaxAngle)))
doubleNumPoints := float64(len(a)) * (1 + doubleNumberOfSegmentsToCreate)
if err := geosegmentize.CheckSegmentizeValidNumPoints(doubleNumPoints, a, b); err != nil {
return nil, err
}
numberOfSegmentsToCreate := int(doubleNumberOfSegmentsToCreate)
numPoints := int(doubleNumPoints)
allSegmentizedCoordinates := make([]float64, 0, numPoints)
allSegmentizedCoordinates = append(allSegmentizedCoordinates, a.Clone()...)
segmentFraction := 1.0 / float64(numberOfSegmentsToCreate)
for pointInserted := 1; pointInserted < numberOfSegmentsToCreate; pointInserted++ {
newPoint := s2.Interpolate(float64(pointInserted)/float64(numberOfSegmentsToCreate), pointA, pointB)
latLng := s2.LatLngFromPoint(newPoint)
allSegmentizedCoordinates = append(allSegmentizedCoordinates, latLng.Lng.Degrees(), latLng.Lat.Degrees())
// Linearly interpolate Z and/or M coordinates.
for i := 2; i < len(a); i++ {
allSegmentizedCoordinates = append(
allSegmentizedCoordinates,
a[i]*(1-float64(pointInserted)*segmentFraction)+b[i]*(float64(pointInserted)*segmentFraction),
)
}
}
return allSegmentizedCoordinates, nil
}