Fix last merge

This commit is contained in:
tidwall 2021-07-11 10:09:51 -07:00
parent 579a41abae
commit dd4d31ae1b
4 changed files with 138 additions and 134 deletions

View File

@ -1,12 +1,10 @@
package collection package collection
import ( import (
"math"
"runtime" "runtime"
"github.com/tidwall/btree" "github.com/tidwall/btree"
"github.com/tidwall/geoindex" "github.com/tidwall/geoindex"
"github.com/tidwall/geoindex/algo"
"github.com/tidwall/geojson" "github.com/tidwall/geojson"
"github.com/tidwall/geojson/geo" "github.com/tidwall/geojson/geo"
"github.com/tidwall/geojson/geometry" "github.com/tidwall/geojson/geometry"
@ -743,7 +741,7 @@ func (c *Collection) Nearby(
} }
nextStep(count, cursor, deadline) nextStep(count, cursor, deadline)
item := itemv.(*itemT) 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 return alive
}, },
) )
@ -759,123 +757,3 @@ func nextStep(step uint64, cursor Cursor, deadline *deadline.Deadline) {
cursor.Step(1) 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, 146164.
// 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Δλ)
}

View File

@ -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, 146164.
// 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Δλ)
}

View File

@ -66,8 +66,11 @@ type ScanWriterParams struct {
o geojson.Object o geojson.Object
fields []float64 fields []float64
distance float64 distance float64
distOutput bool // query or fence requested distance output
noLock bool noLock bool
ignoreGlobMatch bool
clip geojson.Object clip geojson.Object
skipTesting bool
} }
func (s *Server) newScanWriter( func (s *Server) newScanWriter(

View File

@ -12,6 +12,7 @@ import (
"github.com/tidwall/geojson/geometry" "github.com/tidwall/geojson/geometry"
"github.com/tidwall/resp" "github.com/tidwall/resp"
"github.com/tidwall/tile38/internal/bing" "github.com/tidwall/tile38/internal/bing"
"github.com/tidwall/tile38/internal/clip"
"github.com/tidwall/tile38/internal/glob" "github.com/tidwall/tile38/internal/glob"
) )
@ -433,7 +434,10 @@ func (server *Server) cmdNearby(msg *Message) (res resp.Value, err error) {
o: o, o: o,
fields: fields, fields: fields,
distance: meters, distance: meters,
distOutput: s.distance,
noLock: true, noLock: true,
ignoreGlobMatch: true,
skipTesting: true,
}) })
} }
sw.col.Nearby(s.obj, sw, msg.Deadline, iter) sw.col.Nearby(s.obj, sw, msg.Deadline, iter)