計算比給定距離更近的點的數量

Question

我正在開發一個項目，給定一組座標（經度和緯度）需要找到比給定距離（Km）更近的所有點。

我有工作實現遍歷每個點，併爲每個點迭代所有其他點。這是O（n^2）。我想知道我可以採取什麼樣的方法來改善這種情況，考慮到我不想要最近的點，而是比x Km更近的點（我也希望能夠做到相反，找出距離x Km更遠的所有點）。

如果你可以提供一些想法和算法，它會很好。另外代碼示例會很好，特別是在Scala中（我正在開發此項目Scala）。

Answer 1

遲了幾個月回答我自己。那麼我最終使用KD-Tree來實現這一點。我以自己的Java實現爲基礎：http://people.cs.vt.edu/~shaffer/Book/JAVA/progs/KDtree/

但是，我對它進行了一些修改（包括將其移植到Scala），並將其適用於我在Geocaches周圍的特定問題。由此產生的代碼如下：

/* 
* Based on 
* http://people.cs.vt.edu/~shaffer/Book/JAVA/progs/KDtree/ 
* 
* That code is made available as part of this textbook on Data Structures: 
* http://people.cs.vt.edu/~shaffer/Book/ 
*/ 
//A node for a KDTree 
class KDNode[E] { 
    private var key_value : Array[Double] = null //key for this node 
    def key = key_value      
    def key_=(k : Array[Double]) { key_value = k } 

    private var element_value: E = _    //Element for this node 
    def element = element_value 
    def element_=(e : E) {element_value = e} 

    private var left_value : KDNode[E] = null  //Pointer to left child 
    def left = left_value 
    def left_=(l : KDNode[E]) {left_value = l} 

    private var right_value : KDNode[E] = null  //Pointer to right child 
    def right = right_value 
    def right_=(r : KDNode[E]) {right_value = r} 
    /** Constructors */ 
    def this(k : Array[Double], value : E) = 
    { 
    this() 
    key = k 
    element = value 
    } 
    def this(k : Array[Double], value : E, l : KDNode[E], r : KDNode[E]) = 
    { 
    this(k,value) 
    left = l 
    right = r  
    } 
    //Checks if it's a leaf 
    def isLeaf : Boolean = 
    { 
    return (left == null) && (right == null) 
    } 
} 

//Trait for a basic Dictionary that will be used as basis for our KDTree 
trait Dictionary[K,V] { 
    def size: Int 
    def insert(key: K, value: V) 
    def remove(key: K): V 
    def find(key: K): V 
    def clear 
    def isEmpty: Boolean 
} 

/** 
* A extended trait that defines that the key of the dictionary is an array of 
* doubles and defines the need of a method called regionSearchGeo (very 
* important for All Years 5 & 6 stats) 
*/ 
trait GeoDictionary[E] extends Dictionary[Array[Double], E] { 
    def regionSearchGeo(k: Array[Double], radius: Double) : List[E] 
} 

/** 
* Our implementation of a KDTree. It's based on the code provided in the 
* reference above, but had a few key areas modified for geographic searching. 
* For example, this KDtree is a 3D Tree. The keys are Latitude, Longitude and 
* the Year of the cache (that was also included to speed thins even further) 
*/ 
class K3DGeoTree[E] extends GeoDictionary[E]{ 

    private var root : KDNode[E] = null 
    private val D : Int = 3 // Only supporting 2D points in this implementation 
    private var nodecount : Int = 0 // Number of nodes in the KD tree 

    /** 
    * Implementing everything that is required by the Dictionary trait 
    * Some private auxiliary methods are also present 
    */ 
    def clear() = { root = null } 
    def isEmpty() : Boolean = { root == null } 
    def size() : Int = { return nodecount } 

    def insert(pt : Array[Double], value : E) = { 
    root = inserthelp(root, pt, value, 0) 
    nodecount=nodecount+1 
    } 
    def inserthelp(rt : KDNode[E], key : Array[Double], value : E, level : Int) 
           : KDNode[E] = { 
    if (rt == null) return new KDNode[E](key, value) 
    val rtkey : Array[Double] = rt.key 
    if (rtkey(level) > key(level)) 
     rt.left = inserthelp(rt.left, key, value, (level+1)%D) 
    else 
     rt.right = inserthelp(rt.right, key, value, (level+1)%D) 
    return rt 
    } 

    private def findmin(rt : KDNode[E], descrim : Int, level : Int): KDNode[E]= { 
    var temp1 : KDNode[E] = null 
    var temp2 : KDNode[E] = null 
    var key1 : Array[Double] = null 
    var key2 : Array[Double] = null 
    if (rt == null) return null 
    temp1 = findmin(rt.left, descrim, (level+1)%D) 
    if (temp1 != null) key1 = temp1.key 
    if (descrim != level) { 
    temp2 = findmin(rt.right, descrim, (level+1)%D) 
    if (temp2 != null) key2 = temp2.key 
    if ((temp1 == null) || ((temp2 != null) && (key1(descrim) > key2(descrim)))) 
    temp1 = temp2 
    key1 = key2 
    } // Now, temp1 has the smaller value 
    var rtkey : Array[Double] = rt.key 
    if ((temp1 == null) || (key1(descrim) > rtkey(descrim))) 
    return rt 
    else 
    return temp1 
} 

    def find(key : Array[Double]) : E = { return findhelp(root, key, 0) } 

    private def findhelp(rt : KDNode[E], key : Array[Double], level : Int) : E ={ 
    if (rt == null) return null.asInstanceOf[E] 
    val it : E = rt.element 
    val itkey : Array[Double]= rt.key 
    if ((itkey(0) == key(0)) && (itkey(1) == key(1))) 
    return rt.element 
    if (itkey(level) > key(level)) 
    return findhelp(rt.left, key, (level+1)%D) 
    else 
    return findhelp(rt.right, key, (level+1)%D) 
} 

    def remove(key : Array[Double]) : E = { 
    val temp : E = findhelp(root, key, 0) // First find it 
    if (temp != null) { 
      root = removehelp(root, key, 0) // Now remove it 
      nodecount=nodecount-1 
    } 
    return temp 
    } 

    private def removehelp(rt : KDNode[E], key : Array[Double], level : Int) 
              : KDNode[E] = { 
     if (rt == null) return null 
     val rtkey : Array[Double] = rt.key 
     if (key(level) < rtkey(level)) 
     rt.left = removehelp(rt.left, key, (level+1)%D) 
     else if (key(level) > rtkey(level)) 
     rt.right = removehelp(rt.right, key, (level+1)%D) 
     else { // Found it 
     if (rt.right == null) 
     if (rt.left == null) // Just drop element 
     return null 
     else { // Switch subtree to right 
      rt.right = rt.left 
      rt.left = null 
     } 
     val temp : KDNode[E] = findmin(rt.right, level, (level+1)%D) 
     rt.right = removehelp(rt.right, temp.key, (level+1)%D) 
     rt.element = temp.element 
    } 
    return rt 
    } 

    /** 
    * Implementing the GeoDictionary trait 
    */ 
    def regionSearchGeo(point: Array[Double], radius: Double) : List[E] = 
    { 
    val pointGeo : GeoLocation = GeoLocation.fromDegrees(point(0), point(1)) 
    /** 
    * Calculates a bounding rectangle that contains the circle with the given 
    * radius. This will be explained later in the corresponding class 
    */ 
    val boundingRect = pointGeo.boundingCoordinates(radius) 
    //Return the caches found 
    return rsGeoHelp(root, point, radius, boundingRect, 0) 
    } 
    /** 
    * Auxiliary region search function that does all the heavy work 
    */ 
    private def rsGeoHelp(rt : KDNode[E], point : Array[Double], radius : Double, 
               boundingRect : Tuple2[GeoLocation,GeoLocation], 
               lev : Int): List[E] = { 
    if (rt == null) return Nil 
    val rtkey : Array[Double] = rt.key 
    var found : List[E] = Nil 
//Checks if the current node is in the desired radius (and also the year) 
    if (InCircleGeo(point, radius, rtkey)) 
      found = List(rt.element) 
//First Dimension is latitude 
    if(lev % D == 0){ 
     if (rtkey(lev) >= boundingRect._1.degLat) 
     found = found:::rsGeoHelp(rt.left, point, radius, boundingRect, (lev+1)%D) 
     if (rtkey(lev) <= boundingRect._2.degLat) 
     found = found:::rsGeoHelp(rt.right, point, radius, boundingRect, (lev+1)%D) 
    } 
//Second Dimension is Longitude 
    else if(lev % D == 1){ 
     if (rtkey(lev) >= boundingRect._1.degLon) 
     found = found:::rsGeoHelp(rt.left, point, radius, boundingRect, (lev+1)%D) 
     if (rtkey(lev) <= boundingRect._2.degLon) 
    found = found:::rsGeoHelp(rt.right, point, radius, boundingRect, (lev+1)%D) 
    } 
//Third and last dimension is the year 
    else{ 
    found = found:::rsGeoHelp(rt.left, point, radius, boundingRect, (lev+1)%D) 
     if (rtkey(lev) <= point(lev)) 
     found = found:::rsGeoHelp(rt.right, point, radius, boundingRect, (lev+1)%D) 
    } 
    //Return the found nodes (in our case it will be caches) 
    return found 
    } 
    private def InCircleGeo(point : Array[Double], radius : Double, 
                coord : Array[Double]) : Boolean = { 
    //Creates a GeoLocation object for each point 
    val pointGeo : GeoLocation = GeoLocation.fromDegrees(point(0), point(1)) 
    val coordGeo : GeoLocation = GeoLocation.fromDegrees(coord(0), coord(1)) 
    /** 
    * If the year is smaller than the query point and the distance is within 
    * radius return true. Else it's false. 
    */ 
    return (coord(0) != point(0) && coord(1) != point(1) && coord(2) <= point(2) 
     && pointGeo.distanceTo(coordGeo) < radius) 
    } 
} 

/** 
* This class encapsulates a series of utility methods to deal with geographic 
* coordinates. It was based on the information in the link below that gives 
* a very good insight about how to do math with geographic coordinates and 
* also provides some Java samples that we used as an inspiration for this 
* class. 
* Link: http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates 
*/ 
//Companion object of Class GeoLocation to define static methods and variables 
object GeoLocation { 
    //Min maxs in terms of Latitude and Longitude accross the globe 
    private val MIN_LAT : Double = math.toRadians(-90) // -PI/2 
    private val MAX_LAT : Double = math.toRadians(90) // PI/2 
    private val MIN_LON : Double = math.toRadians(-180) // -PI 
    private val MAX_LON : Double = math.toRadians(180) // PI 
    /** 
    * Earth radius. This value is the most used but there are others that may 
    * give slightly different results. 
    */ 
    private val RADIUS : Double = 6372.8 

    /** 
    * A factory method that creates a GeoLocation object from given latitude and 
    * longitude in degrees 
    */ 
    def fromDegrees(latitude : Double, longitude : Double) : GeoLocation = { 
    val result : GeoLocation = new GeoLocation() 
     result.radLat = math.toRadians(latitude) 
     result.radLon = math.toRadians(longitude) 
     result.degLat = latitude 
     result.degLon = longitude 
     result.checkBounds 
     return result 
    } 

    /** 
    * A factory method that creates a GeoLocation object from given latitude and 
    * longitude in radians 
    */ 
    def fromRadians(latitude : Double, longitude : Double) : GeoLocation = { 
     val result : GeoLocation = new GeoLocation() 
     result.radLat = latitude 
     result.radLon = longitude 
     result.degLat = math.toDegrees(latitude) 
     result.degLon = math.toDegrees(longitude) 
     result.checkBounds 
     return result 
    } 
} 
/** 
* The GeoLocation class itself. The constructor is private use the factory 
* methods above. 
*/ 
class GeoLocation private{ 

    /** 
    * Getters and Setters implemented as properties with syntactic sugar 
    * This properties contain the latitude and longitude in degrees and radians 
    */ 
    private var radLat_value : Double = _ 
    def radLat = radLat_value      
    private def radLat_=(k : Double) { radLat_value = k } 

    private var radLon_value : Double = _ 
    def radLon = radLon_value      
    private def radLon_=(k : Double) { radLon_value = k } 

    private var degLat_value : Double = _ 
    def degLat = degLat_value      
    private def degLat_=(k : Double) { degLat_value = k } 

    private var degLon_value : Double = _ 
    def degLon = degLon_value      
    private def degLon_=(k : Double) { degLon_value = k } 

    /** 
    * Check if the vales are valid considering the MIN and MAX for latitude and 
    * longitude. 
    */ 
    private def checkBounds = { 
     if (radLat < GeoLocation.MIN_LAT || radLat > GeoLocation.MAX_LAT || 
       radLon < GeoLocation.MIN_LON || radLon > GeoLocation.MAX_LON) 
      throw new IllegalArgumentException() 
    } 

    /** 
    * Function to calculate the distance between this GeoLocation and the given 
    * GeoLocation. 
    * 
    * Check the reference above and 
    * http://en.wikipedia.org/wiki/Haversine_formula 
    * for more information. 
    */ 
    def distanceTo(location : GeoLocation) : Double = { 
     return math.acos(math.sin(radLat) * math.sin(location.radLat) + 
        math.cos(radLat) * math.cos(location.radLat) * 
        math.cos(radLon - location.radLon)) * GeoLocation.RADIUS 
    } 

    /** 
    * This method is very important for the search made in the K3DTree. 
    * It allows us to make a bouding rectangle with the given distance/radius 
    * that is geometrically correct. Check the reference above to learn more 
    * about the math involved. 
    */ 
    def boundingCoordinates(distance : Double) 
     : Tuple2[GeoLocation, GeoLocation] = { 
    if (distance < 0d) throw new IllegalArgumentException() 
     // Angular distance in radians on a great circle 
     val radDist : Double = distance/GeoLocation.RADIUS 

     //Initialize local variables to check for poles 
     var minLat : Double = radLat - radDist 
     var maxLat : Double = radLat + radDist 

     var minLon : Double = 0 
     var maxLon : Double = 0 

     //Normal case 
     if (minLat > GeoLocation.MIN_LAT && maxLat < GeoLocation.MAX_LAT) { 
       val deltaLon : Double = math.asin(math.sin(radDist)/math.cos(radLat)) 
         minLon = radLon - deltaLon 
       if (minLon < GeoLocation.MIN_LON) minLon += 2d * math.Pi 
         maxLon = radLon + deltaLon 
       if (maxLon > GeoLocation.MAX_LON) maxLon -= 2d * math.Pi 
     } 
    //Special case in which a pole is within the distance 
     else{ 
         minLat = math.max(minLat, GeoLocation.MIN_LAT) 
         maxLat = math.min(maxLat, GeoLocation.MAX_LAT) 
         minLon = GeoLocation.MIN_LON 
         maxLon = GeoLocation.MAX_LON 
       }  
    /** 
    * Each of the bounding points (one in the south-west, bottom-left, 
    * and other in the north-east, top-right) 
    */ 
    val swPoint : GeoLocation = GeoLocation.fromRadians(minLat, minLon) 
    val nePoint : GeoLocation = GeoLocation.fromRadians(maxLat, maxLon) 
    //Return the tuple with the two points 
    return (swPoint, nePoint) 
    } 
} 

整個代碼被記錄，所以我希望這可以幫助有類似問題的人。在這個特定的問題中，除了緯度和經度之外，我還要處理幾年，所以我增加了一個額外的維度。但是對於一個更普遍的地理問題，只用兩個維度（一個用於緯度，一個用於經度）更容易。

Answer 2

對於這個問題，我會使用一個java庫（這可能在斯卡拉）。 計算地球上的距離比您想象的要困難得多。例如，地球不是一個完美的球體。 在java中做'geo stuff'的主要庫是JTS：http://tsusiatsoftware.net/jts/main.html 你也可以看一下Geotools或甚至像PostGIS（建立在Postgresql上）的GIS數據庫，但這可能對你的項目來說過於矯枉過正