搜索附近
最近开发中用到位置搜索附近功能，看到网上有流行的Geohash实现，也有KD树实现。
分别写了实现算法进行验证，感觉既复杂又低效（尤其是难以理解）。
于是想能不能更直观高效的算法实现呢。
假设地球表面展开是一个正方形 -180~180，-90~90，（可以修正）
把这个正方向分割成许许多多的小正方型（一个正方形认为是一个点，不考虑立体经纬度重合问题）
当计算附近的时候，首先计算中心点所在的正方形，再根据距离计算出要找的周围正方形的范围（绿色区域），
然后遍历计算一个正方形区域（蓝色区域），获得距离内的待搜索目标（红色区域内）
图一
地球周长： 6371*2*PI*1000 （m）
用30bit存贮： 2^30=1342177279
计算 6371*2*PI*1000/2^30=0.030（m）
共需要2^30*2^30=2^60= 1152921504606846975 个正方形
这个原理简单明了，但是有两个致命缺陷：
1 正方形太多，没有哪台计算机能存的下这么多正方形
2 确定中心后，周围待查询的小正方形太多，效率太低。
针对问题1 ，可以通过稀疏矩阵来解决
图二
针对问题2，暂时没有办法。
改进一下：
1 30bit分两级存贮，0~14一级，15~29二级 (这里要根据情况需要调整)
2 用正方形左下角的点表示这个正方形。（由于所有正方形的宽、高都是相同的）
3 只存贮有对象的正方形，没对象的不存
图三
如上图 一级正方行里面保存二级正方形，二级正方形保存实际的对象（m，n可以根据实际情况选合适的值）。
1 一级，二级均通过稀疏矩阵存贮。
2 一级矩阵由于命中概率较高，采用for循环行列（如图一的蓝色区域）
3 二级矩阵命中率较低，采用稀疏矩阵遍历。
4 把稀疏矩阵的行列拼接成一个整数（long）。
算法实现：
假定 m=15，n=15；
计算 6371*2*PI/2^15=1.22（km）（一级表示的精度是1.22km）
6371*2*PI*1000/2^30=0.030（m） （二级表示的精度是0.030m）
把所有待搜索目标计算经纬度，放置到这些正方形里面。（计算方法同Geohash）
地球纬度区间是[-90,90]，中关村大厦的纬度是39.97694，可以通过下面算法对纬度39.97694进行逼近编码:
1）区间[-90,90]进行二分为[-90,0),[0,90]，称为左右区间，可以确定39.97694属于右区间[0,90]，给标记为1；
2）接着将区间[0,90]进行二分为 [0,45),[45,90]，可以确定39.97694属于左区间 [0,45)，给标记为0；
3）递归上述过程39.928167总是属于某个区间[a,b]。随着每次迭代区间[a,b]总在缩小，并越来越逼近39.97694；
4）如果给定的纬度x（39.97694）属于左区间，则记录0，如果属于右区间则记录1，这样随着算法的进行会产生一个序列101110001101101，记作：lat。
int lat=getLatCode(39.97694)=775342647; (101110001101101100101000110111)
地球经度区间是[-90,90]，中关村大厦的经度是116.31，同样的原理计算出编码为：110100101011010，记作：lng
int lng=getLngCode(116.31)=883778999; (110100101011010110010110110111)
int MAXIMIZE=0x40000000;
int DATA_BITS=30;
int DATA_BITS_FIRST=15;
LONGITUDE_RATIO=MAXIMIZE/360.0;
double LATITUDE_RATIO=MAXIMIZE/180.0;
//Geohash 经度编码
int getLngCode(double lng){
int result=0;
double start=-180,end=180,midd;
for(int i=DATA_BITS;i>0;i--){
result<<=1;
midd=(start+end)/2;
if(lng>midd){
result|=1;
start=midd;
}
else{
end=midd;
}
}
return result;
}
//Geohash 纬度编码
int getLatCode(double lat){
int result=0;
double start=-90,end=90,midd;
for(int i=DATA_BITS-1;i>0;i--){
result<<=1;
midd=(start+end)/2;
if(lat>midd){
result|=1;
start=midd;
}
else{
end=midd;
}
}
return result;
}
分别取经度、纬度的高15位作为行、列构造一级稀疏矩阵，然后将行列合成。
经纬度合成： 把经度高15bit左移16bit | 纬度高15bit ，记作：key1 （一级正方形左下角的点）
int key1=（((lng>>15)&0x7fff)<<16)|((lat>>15)&0x7fff); （一级KEY）
把经度左移32bit | 纬度 ，记作：key 2（二级正方形左下角的点）
long key2=（lng<<32)|lat; （二级KEY）
这样稀疏矩阵就可以通过Hash表来存贮：
Hashtable<key1,二级稀疏矩阵>
二级稀疏矩阵: Hashtable<Object,key2>
详细算法实现。
import java.util.HashMap;
import java.util.Map;
public class Position {
//2^30
protected static final int MAXIMIZE=0x40000000;
protected static final int DATA_BITS=30;
public static final int LONGITUDE_MAXIMIZE=MAXIMIZE;
public static final int LATITUDE_MAXIMIZE=(MAXIMIZE>>1);
protected static final double EARTH_PERIMETER=6378.137*2000*Math.PI;
protected static final double EARTH_DISTANCE=EARTH_PERIMETER/360;
protected static final double DATA_RATIO=MAXIMIZE/360.0;
protected static final double EARTH_RATIO=EARTH_PERIMETER/MAXIMIZE;
protected static final int DATA_MASK=0x3FFFFFFF;
protected static final int HIGH_MASK=0xFFFF;
protected static final int HIGH_BITS=16;
protected static final int LOW_BITS=14;
//浮点型经纬度转换为整型
public static int getLngCode(double pos){
return (int)(DATA_RATIO*(pos+180));
}
//整型经纬度转换浮点型
public static double getLongitude(int code){
return code/DATA_RATIO-180;
}
//根据浮点纬度计算编码
public static int getLatCode(double lat){
return (int)(DATA_RATIO*(lat+90));
}
//根据编码计算浮点纬度
public static double getLatitude(int code){
return code/DATA_RATIO-90;
}
//计算两点间距离(单位：米，勾股定理)
public static int getDistance(double lng1,double lat1,double lng2,double lat2){
lng1-=lng2;
//修正曲面引起的误差 r=R*cos(纬度)
lng1*=Math.cos((lat1+lat2)*Math.PI/360);
lat1-=lat2;
return (int)(Math.sqrt(lng1*lng1+lat1*lat1)*EARTH_DISTANCE);
}
//计算两点间距离(单位：米，勾股定理)
public static int getDistance(int lng1,int lat1,int lng2,int lat2){
long lg=lng1-lng2;
//修正曲面引起的误差 r=R*cos(纬度-90)=R*sin(纬度)
lg*=Math.sin(((lat1+lat2)*Math.PI/LATITUDE_MAXIMIZE/2));
long lt=lat1-lat2;
return (int)(Math.sqrt((lg*lg)+((lt*lt)))*EARTH_RATIO);
}
private Map<Integer,Map<Object,Long>> table=new HashMap<Integer,Map<Object,Long>>();
//插入对象
public void insert(double lng,double lat,Object code){
insert(getLngCode(lng),getLatCode(lat),code);
}
public void insert(int lg,int lt,Object code){
int key=((lt>>LOW_BITS)<<HIGH_BITS)|(lg>>LOW_BITS);
Map<Object,Long> value=table.get(key);
if(value==null){
value=new HashMap<Object,Long>();
table.put(key, value);
}
long l=lt;
//这里value是 ((纬度<<32)|经度)
value.put(code,((l<<32)|lg));
}
//删除对象
public void remove(double lng,double lat,Object code){
remove(getLngCode(lng),getLatCode(lat),code);
}
public void remove(int lg,int lt,Object code){
int key=((lt>>LOW_BITS)<<HIGH_BITS)|(lg>>LOW_BITS);
Map<Object,Long> value=table.get(key);
if(value!=null){
value.remove(code);
//不再有任何对象，直接删除一级点
if(value.size()==0){
table.remove(key);
}
}
}
//查询对象
public Map<Object,Long> search(double lng,double lat,int distance){
return search(getLngCode(lng),getLatCode(lat),distance);
}
public Map<Object,Long> search(int cntX,int cntY,int distance){
int dist=(int)(distance/EARTH_RATIO);
//计算纬度开始，结束值
int startY=((cntY-dist)>>LOW_BITS),endY=((cntY+dist)>>LOW_BITS);
dist/=Math.sin((cntY*Math.PI/LATITUDE_MAXIMIZE));
//计算经度开始，结束值
int startX=((cntX-dist)>>LOW_BITS),endX=((cntX+dist)>>LOW_BITS);
Map<Object,Long> map=new HashMap<Object,Long>();
//循环遍历蓝色格子
for(int lt=startY;lt<=endY;lt++){
for(int lg=startX;lg<=endX;lg++){
int key=((lt<<HIGH_BITS)|lg);
Map<Object,Long> value=table.get(key);
if(value!=null){
//遍历二级点
for(Map.Entry<Object, Long> it:value.entrySet()){
long t=it.getValue();
if(getDistance((int)(t&DATA_MASK),(int)(t>>32),cntX,cntY)<distance){
map.put(it.getKey(), t);
}
}
}
}
}
return map;
}
}
注意事项：
1 一级、二级的划分，根据实际情况合理划分一二级，比如二级的对象太多（一般几十个点比较合适），就应该放大一级的位数，反之缩小一级的位数
浙公网安备 33010602011771号