diff --git a/internal/collection/collection.go b/internal/collection/collection.go index df1f75b4..ee4ebe74 100644 --- a/internal/collection/collection.go +++ b/internal/collection/collection.go @@ -1,12 +1,10 @@ package collection import ( - "math" "runtime" "github.com/tidwall/btree" "github.com/tidwall/geoindex" - "github.com/tidwall/geoindex/algo" "github.com/tidwall/geojson" "github.com/tidwall/geojson/geo" "github.com/tidwall/geojson/geometry" @@ -743,7 +741,7 @@ func (c *Collection) Nearby( } nextStep(count, cursor, deadline) item := itemv.(*itemT) - alive = iter(item.id, item.obj, c.getFieldValues(item.id), dist) + alive = iter(item.id, item.obj, c.fieldValues.get(item.fieldValuesSlot), dist) return alive }, ) @@ -759,123 +757,3 @@ func nextStep(step uint64, cursor Cursor, deadline *deadline.Deadline) { cursor.Step(1) } } - -func geodeticDistAlgo(center [2]float64) func( - min, max [2]float64, data interface{}, item bool, - add func(min, max [2]float64, data interface{}, item bool, dist float64), -) { - const earthRadius = 6371e3 - return func( - min, max [2]float64, data interface{}, item bool, - add func(min, max [2]float64, data interface{}, item bool, dist float64), - ) { - add(min, max, data, item, earthRadius*pointRectDistGeodeticDeg( - center[1], center[0], - min[1], min[0], - max[1], max[0], - )) - } -} - -func pointRectDistGeodeticDeg(pLat, pLng, minLat, minLng, maxLat, maxLng float64) float64 { - result := pointRectDistGeodeticRad( - pLat*math.Pi/180, pLng*math.Pi/180, - minLat*math.Pi/180, minLng*math.Pi/180, - maxLat*math.Pi/180, maxLng*math.Pi/180, - ) - return result -} - -func pointRectDistGeodeticRad(φq, λq, φl, λl, φh, λh float64) float64 { - // Algorithm from: - // Schubert, E., Zimek, A., & Kriegel, H.-P. (2013). - // Geodetic Distance Queries on R-Trees for Indexing Geographic Data. - // Lecture Notes in Computer Science, 146–164. - // doi:10.1007/978-3-642-40235-7_9 - const ( - twoΠ = 2 * math.Pi - halfΠ = math.Pi / 2 - ) - - // distance on the unit sphere computed using Haversine formula - distRad := func(φa, λa, φb, λb float64) float64 { - if φa == φb && λa == λb { - return 0 - } - - Δφ := φa - φb - Δλ := λa - λb - sinΔφ := math.Sin(Δφ / 2) - sinΔλ := math.Sin(Δλ / 2) - cosφa := math.Cos(φa) - cosφb := math.Cos(φb) - - return 2 * math.Asin(math.Sqrt(sinΔφ*sinΔφ+sinΔλ*sinΔλ*cosφa*cosφb)) - } - - // Simple case, point or invalid rect - if φl >= φh && λl >= λh { - return distRad(φl, λl, φq, λq) - } - - if λl <= λq && λq <= λh { - // q is between the bounding meridians of r - // hence, q is north, south or within r - if φl <= φq && φq <= φh { // Inside - return 0 - } - - if φq < φl { // South - return φl - φq - } - - return φq - φh // North - } - - // determine if q is closer to the east or west edge of r to select edge for - // tests below - Δλe := λl - λq - Δλw := λq - λh - if Δλe < 0 { - Δλe += twoΠ - } - if Δλw < 0 { - Δλw += twoΠ - } - var Δλ float64 // distance to closest edge - var λedge float64 // longitude of closest edge - if Δλe <= Δλw { - Δλ = Δλe - λedge = λl - } else { - Δλ = Δλw - λedge = λh - } - - sinΔλ, cosΔλ := math.Sincos(Δλ) - tanφq := math.Tan(φq) - - if Δλ >= halfΠ { - // If Δλ > 90 degrees (1/2 pi in radians) we're in one of the corners - // (NW/SW or NE/SE depending on the edge selected). Compare against the - // center line to decide which case we fall into - φmid := (φh + φl) / 2 - if tanφq >= math.Tan(φmid)*cosΔλ { - return distRad(φq, λq, φh, λedge) // North corner - } - return distRad(φq, λq, φl, λedge) // South corner - } - - if tanφq >= math.Tan(φh)*cosΔλ { - return distRad(φq, λq, φh, λedge) // North corner - } - - if tanφq <= math.Tan(φl)*cosΔλ { - return distRad(φq, λq, φl, λedge) // South corner - } - - // We're to the East or West of the rect, compute distance using cross-track - // Note that this is a simplification of the cross track distance formula - // valid since the track in question is a meridian. - return math.Asin(math.Cos(φq) * sinΔλ) -} diff --git a/internal/collection/geodesic.go b/internal/collection/geodesic.go new file mode 100644 index 00000000..d531b01c --- /dev/null +++ b/internal/collection/geodesic.go @@ -0,0 +1,119 @@ +package collection + +import "math" + +func geodeticDistAlgo(center [2]float64) ( + algo func(min, max [2]float64, data interface{}, item bool) (dist float64), +) { + const earthRadius = 6371e3 + return func(min, max [2]float64, data interface{}, item bool) (dist float64) { + return earthRadius * pointRectDistGeodeticDeg( + center[1], center[0], + min[1], min[0], + max[1], max[0], + ) + } +} + +func pointRectDistGeodeticDeg(pLat, pLng, minLat, minLng, maxLat, maxLng float64) float64 { + result := pointRectDistGeodeticRad( + pLat*math.Pi/180, pLng*math.Pi/180, + minLat*math.Pi/180, minLng*math.Pi/180, + maxLat*math.Pi/180, maxLng*math.Pi/180, + ) + return result +} + +func pointRectDistGeodeticRad(φq, λq, φl, λl, φh, λh float64) float64 { + // Algorithm from: + // Schubert, E., Zimek, A., & Kriegel, H.-P. (2013). + // Geodetic Distance Queries on R-Trees for Indexing Geographic Data. + // Lecture Notes in Computer Science, 146–164. + // doi:10.1007/978-3-642-40235-7_9 + const ( + twoΠ = 2 * math.Pi + halfΠ = math.Pi / 2 + ) + + // distance on the unit sphere computed using Haversine formula + distRad := func(φa, λa, φb, λb float64) float64 { + if φa == φb && λa == λb { + return 0 + } + + Δφ := φa - φb + Δλ := λa - λb + sinΔφ := math.Sin(Δφ / 2) + sinΔλ := math.Sin(Δλ / 2) + cosφa := math.Cos(φa) + cosφb := math.Cos(φb) + + return 2 * math.Asin(math.Sqrt(sinΔφ*sinΔφ+sinΔλ*sinΔλ*cosφa*cosφb)) + } + + // Simple case, point or invalid rect + if φl >= φh && λl >= λh { + return distRad(φl, λl, φq, λq) + } + + if λl <= λq && λq <= λh { + // q is between the bounding meridians of r + // hence, q is north, south or within r + if φl <= φq && φq <= φh { // Inside + return 0 + } + + if φq < φl { // South + return φl - φq + } + + return φq - φh // North + } + + // determine if q is closer to the east or west edge of r to select edge for + // tests below + Δλe := λl - λq + Δλw := λq - λh + if Δλe < 0 { + Δλe += twoΠ + } + if Δλw < 0 { + Δλw += twoΠ + } + var Δλ float64 // distance to closest edge + var λedge float64 // longitude of closest edge + if Δλe <= Δλw { + Δλ = Δλe + λedge = λl + } else { + Δλ = Δλw + λedge = λh + } + + sinΔλ, cosΔλ := math.Sincos(Δλ) + tanφq := math.Tan(φq) + + if Δλ >= halfΠ { + // If Δλ > 90 degrees (1/2 pi in radians) we're in one of the corners + // (NW/SW or NE/SE depending on the edge selected). Compare against the + // center line to decide which case we fall into + φmid := (φh + φl) / 2 + if tanφq >= math.Tan(φmid)*cosΔλ { + return distRad(φq, λq, φh, λedge) // North corner + } + return distRad(φq, λq, φl, λedge) // South corner + } + + if tanφq >= math.Tan(φh)*cosΔλ { + return distRad(φq, λq, φh, λedge) // North corner + } + + if tanφq <= math.Tan(φl)*cosΔλ { + return distRad(φq, λq, φl, λedge) // South corner + } + + // We're to the East or West of the rect, compute distance using cross-track + // Note that this is a simplification of the cross track distance formula + // valid since the track in question is a meridian. + return math.Asin(math.Cos(φq) * sinΔλ) +} diff --git a/internal/server/scanner.go b/internal/server/scanner.go index 6ec6dc2c..e11c14ad 100644 --- a/internal/server/scanner.go +++ b/internal/server/scanner.go @@ -62,12 +62,15 @@ type scanWriter struct { // ScanWriterParams ... type ScanWriterParams struct { - id string - o geojson.Object - fields []float64 - distance float64 - noLock bool - clip geojson.Object + id string + o geojson.Object + fields []float64 + distance float64 + distOutput bool // query or fence requested distance output + noLock bool + ignoreGlobMatch bool + clip geojson.Object + skipTesting bool } func (s *Server) newScanWriter( diff --git a/internal/server/search.go b/internal/server/search.go index 7d4b3e90..9d647b5b 100644 --- a/internal/server/search.go +++ b/internal/server/search.go @@ -12,6 +12,7 @@ import ( "github.com/tidwall/geojson/geometry" "github.com/tidwall/resp" "github.com/tidwall/tile38/internal/bing" + "github.com/tidwall/tile38/internal/clip" "github.com/tidwall/tile38/internal/glob" ) @@ -429,11 +430,14 @@ func (server *Server) cmdNearby(msg *Message) (res resp.Value, err error) { meters = dist } return sw.writeObject(ScanWriterParams{ - id: id, - o: o, - fields: fields, - distance: meters, - noLock: true, + id: id, + o: o, + fields: fields, + distance: meters, + distOutput: s.distance, + noLock: true, + ignoreGlobMatch: true, + skipTesting: true, }) } sw.col.Nearby(s.obj, sw, msg.Deadline, iter)