lucene-java-3.4.0/lucene/contrib/spatial/src/java/org/apache/lucene/spatial/DistanceUtils.java

   1 /**
   2  * Licensed to the Apache Software Foundation (ASF) under one or more
   3  * contributor license agreements.  See the NOTICE file distributed with
   4  * this work for additional information regarding copyright ownership.
   5  * The ASF licenses this file to You under the Apache License, Version 2.0
   6  * (the "License"); you may not use this file except in compliance with
   7  * the License.  You may obtain a copy of the License at
   8  *
   9  *     http://www.apache.org/licenses/LICENSE-2.0
  10  *
  11  * Unless required by applicable law or agreed to in writing, software
  12  * distributed under the License is distributed on an "AS IS" BASIS,
  13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14  * See the License for the specific language governing permissions and
  15  * limitations under the License.
  16  */
  17
  18 package org.apache.lucene.spatial;
  19
  20 import org.apache.lucene.spatial.geometry.DistanceUnits;
  21 import org.apache.lucene.spatial.geometry.FloatLatLng;
  22 import org.apache.lucene.spatial.geometry.LatLng;
  23 import org.apache.lucene.spatial.geometry.shape.LLRect;
  24 import org.apache.lucene.spatial.geometry.shape.Rectangle;
  25 import org.apache.lucene.spatial.tier.InvalidGeoException;
  26
  27 /**
  28  * <p><font color="red"><b>NOTE:</b> This API is still in
  29  * flux and might change in incompatible ways in the next
  30  * release.</font>
  31  */
  32
  33 public class DistanceUtils {
  34
  35   public static final double DEGREES_TO_RADIANS = Math.PI / 180.0;
  36   public static final double RADIANS_TO_DEGREES = 180.0 / Math.PI;
  37   //pre-compute some angles that are commonly used
  38   public static final double DEG_45_AS_RADS = Math.PI / 4.0;
  39   public static final double SIN_45_AS_RADS = Math.sin(DEG_45_AS_RADS);
  40   public static final double DEG_90_AS_RADS = Math.PI / 2;
  41   public static final double DEG_180_AS_RADS = Math.PI;
  42   public static final double DEG_225_AS_RADS = 5 * DEG_45_AS_RADS;
  43   public static final double DEG_270_AS_RADS = 3*DEG_90_AS_RADS;
  44
  45
  46   public static final double KM_TO_MILES = 0.621371192;
  47   public static final double MILES_TO_KM = 1.609344;
  48     /**
  49    * The International Union of Geodesy and Geophysics says the Earth's mean radius in KM is:
  50    *
  51    * [1] http://en.wikipedia.org/wiki/Earth_radius
  52    */
  53   public static final double EARTH_MEAN_RADIUS_KM = 6371.009;
  54
  55   public static final double EARTH_MEAN_RADIUS_MI = EARTH_MEAN_RADIUS_KM / MILES_TO_KM;
  56
  57   public static final double EARTH_EQUATORIAL_RADIUS_MI = 3963.205;
  58   public static final double EARTH_EQUATORIAL_RADIUS_KM = EARTH_EQUATORIAL_RADIUS_MI * MILES_TO_KM;
  59
  60
  61   public static double getDistanceMi(double x1, double y1, double x2, double y2) {
  62     return getLLMDistance(x1, y1, x2, y2);
  63   }
  64
  65   /**
  66    *
  67    * @param x1
  68    * @param y1
  69    * @param miles
  70    * @return boundary rectangle where getY/getX is top left, getMinY/getMinX is bottom right
  71    */
  72   public static Rectangle getBoundary (double x1, double y1, double miles) {
  73
  74     LLRect box = LLRect.createBox( new FloatLatLng( x1, y1 ), miles, miles );
  75
  76     //System.out.println("Box: "+maxX+" | "+ maxY+" | "+ minX + " | "+ minY);
  77     return box.toRectangle();
  78
  79   }
  80
  81   public static double getLLMDistance (double x1, double y1, double x2, double y2) {
  82
  83     LatLng p1 = new FloatLatLng( x1, y1 );
  84     LatLng p2 = new FloatLatLng( x2, y2 );
  85     return p1.arcDistance( p2, DistanceUnits.MILES );
  86   }
  87
  88   /**
  89    * distance/radius.
  90    * @param distance The distance travelled
  91    * @param radius The radius of the sphere
  92    * @return The angular distance, in radians
  93    */
  94   public static double angularDistance(double distance, double radius){
  95     return distance/radius;
  96   }
  97
  98   /**
  99    * Calculate the p-norm (i.e. length) beteen two vectors
 100    *
 101    * @param vec1  The first vector
 102    * @param vec2  The second vector
 103    * @param power The power (2 for Euclidean distance, 1 for manhattan, etc.)
 104    * @return The length.
 105    *         <p/>
 106    *         See http://en.wikipedia.org/wiki/Lp_space
 107    * @see #vectorDistance(double[], double[], double, double)
 108    */
 109   public static double vectorDistance(double[] vec1, double[] vec2, double power) {
 110     return vectorDistance(vec1, vec2, power, 1.0 / power);
 111   }
 112
 113   /**
 114    * Calculate the p-norm (i.e. length) between two vectors
 115    *
 116    * @param vec1         The first vector
 117    * @param vec2         The second vector
 118    * @param power        The power (2 for Euclidean distance, 1 for manhattan, etc.)
 119    * @param oneOverPower If you've precalculated oneOverPower and cached it, use this method to save one division operation over {@link #vectorDistance(double[], double[], double)}.
 120    * @return The length.
 121    */
 122   public static double vectorDistance(double[] vec1, double[] vec2, double power, double oneOverPower) {
 123     double result = 0;
 124
 125     if (power == 0) {
 126       for (int i = 0; i < vec1.length; i++) {
 127         result += vec1[i] - vec2[i] == 0 ? 0 : 1;
 128       }
 129
 130     } else if (power == 1.0) {
 131       for (int i = 0; i < vec1.length; i++) {
 132         result += vec1[i] - vec2[i];
 133       }
 134     } else if (power == 2.0) {
 135       result = Math.sqrt(squaredEuclideanDistance(vec1, vec2));
 136     } else if (power == Integer.MAX_VALUE || Double.isInfinite(power)) {//infinite norm?
 137       for (int i = 0; i < vec1.length; i++) {
 138         result = Math.max(result, Math.max(vec1[i], vec2[i]));
 139       }
 140     } else {
 141       for (int i = 0; i < vec1.length; i++) {
 142         result += Math.pow(vec1[i] - vec2[i], power);
 143       }
 144       result = Math.pow(result, oneOverPower);
 145     }
 146     return result;
 147   }
 148
 149   /**
 150    * Return the coordinates of a vector that is the corner of a box (upper right or lower left), assuming a Rectangular
 151    * coordinate system.  Note, this does not apply for points on a sphere or ellipse (although it could be used as an approximatation).
 152    *
 153    * @param center     The center point
 154    * @param result Holds the result, potentially resizing if needed.
 155    * @param distance   The d from the center to the corner
 156    * @param upperRight If true, return the coords for the upper right corner, else return the lower left.
 157    * @return The point, either the upperLeft or the lower right
 158    */
 159   public static double[] vectorBoxCorner(double[] center, double[] result, double distance, boolean upperRight) {
 160     if (result == null || result.length != center.length) {
 161       result = new double[center.length];
 162     }
 163     if (upperRight == false) {
 164       distance = -distance;
 165     }
 166     //We don't care about the power here,
 167     // b/c we are always in a rectangular coordinate system, so any norm can be used by
 168     //using the definition of sine
 169     distance = SIN_45_AS_RADS * distance; // sin(Pi/4) == (2^0.5)/2 == opp/hyp == opp/distance, solve for opp, similarily for cosine
 170     for (int i = 0; i < center.length; i++) {
 171       result[i] = center[i] + distance;
 172     }
 173     return result;
 174   }
 175
 176   /**
 177    * @param latCenter  In degrees
 178    * @param lonCenter  In degrees
 179    * @param distance The distance
 180    * @param result A preallocated array to hold the results.  If null, a new one is constructed.
 181    * @param upperRight If true, calculate the upper right corner, else the lower left
 182    * @param sphereRadius The radius of the sphere to use.
 183    * @return The Lat/Lon in degrees
 184    *
 185    * @see #latLonCorner(double, double, double, double[], boolean, double)
 186    */
 187   public static double[] latLonCornerDegs(double latCenter, double lonCenter,
 188                                           double distance, double [] result,
 189                                           boolean upperRight, double sphereRadius) {
 190     result = latLonCorner(latCenter * DEGREES_TO_RADIANS,
 191             lonCenter * DEGREES_TO_RADIANS, distance, result, upperRight, sphereRadius);
 192     result[0] = result[0] * RADIANS_TO_DEGREES;
 193     result[1] = result[1] * RADIANS_TO_DEGREES;
 194     return result;
 195   }
 196
 197   /**
 198    * Uses Haversine to calculate the corner of a box (upper right or lower left) that is the <i>distance</i> away, given a sphere of the specified <i>radius</i>.
 199    *
 200    * NOTE: This is not the same as calculating a box that transcribes a circle of the given distance.
 201    *
 202    * @param latCenter  In radians
 203    * @param lonCenter  In radians
 204    * @param distance   The distance
 205    * @param result A preallocated array to hold the results.  If null, a new one is constructed.
 206    * @param upperRight If true, give lat/lon for the upper right corner, else lower left
 207    * @param sphereRadius     The radius to use for the calculation
 208    * @return The Lat/Lon in Radians
 209
 210    */
 211   public static double[] latLonCorner(double latCenter, double lonCenter,
 212                                       double distance, double [] result, boolean upperRight, double sphereRadius) {
 213     // Haversine formula
 214     double brng = upperRight ? DEG_45_AS_RADS : DEG_225_AS_RADS;
 215     result = pointOnBearing(latCenter, lonCenter, distance, brng, result, sphereRadius);
 216
 217     return result;
 218   }
 219
 220   /**
 221    * Given a start point (startLat, startLon) and a bearing on a sphere of radius <i>sphereRadius</i>, return the destination point.
 222    * @param startLat The starting point latitude, in radians
 223    * @param startLon The starting point longitude, in radians
 224    * @param distance The distance to travel along the bearing.  The units are assumed to be the same as the sphereRadius units, both of which is up to the caller to know
 225    * @param bearing The bearing, in radians.  North is a 0 deg. bearing, east is 90 deg, south is 180 deg, west is 270 deg.
 226    * @param result A preallocated array to hold the results.  If null, a new one is constructed.
 227    * @param sphereRadius The radius of the sphere to use for the calculation.
 228    * @return The destination point, in radians.  First entry is latitude, second is longitude
 229    */
 230   public static double[] pointOnBearing(double startLat, double startLon, double distance, double bearing, double[] result, double sphereRadius) {
 231     /*
 232         lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
 233         lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))
 234
 235      */
 236     double cosAngDist = Math.cos(distance / sphereRadius);
 237     double cosStartLat = Math.cos(startLat);
 238     double sinAngDist = Math.sin(distance / sphereRadius);
 239     double lat2 = Math.asin(Math.sin(startLat) * cosAngDist +
 240             cosStartLat * sinAngDist * Math.cos(bearing));
 241
 242     double lon2 = startLon + Math.atan2(Math.sin(bearing) * sinAngDist * cosStartLat,
 243             cosAngDist - Math.sin(startLat) * Math.sin(lat2));
 244
 245     /*lat2 = (lat2*180)/Math.PI;
 246     lon2 = (lon2*180)/Math.PI;*/
 247     //From Lucene.  Move back to Lucene when synced
 248     // normalize long first
 249     if (result == null || result.length != 2){
 250       result = new double[2];
 251     }
 252     result[0] = lat2;
 253     result[1] = lon2;
 254     normLng(result);
 255
 256     // normalize lat - could flip poles
 257     normLat(result);
 258     return result;
 259   }
 260
 261   /**
 262    * @param latLng The lat/lon, in radians. lat in position 0, long in position 1
 263    */
 264   public static void normLat(double[] latLng) {
 265
 266     if (latLng[0] > DEG_90_AS_RADS) {
 267       latLng[0] = DEG_90_AS_RADS - (latLng[0] - DEG_90_AS_RADS);
 268       if (latLng[1] < 0) {
 269         latLng[1] = latLng[1] + DEG_180_AS_RADS;
 270       } else {
 271         latLng[1] = latLng[1] - DEG_180_AS_RADS;
 272       }
 273     } else if (latLng[0] < -DEG_90_AS_RADS) {
 274       latLng[0] = -DEG_90_AS_RADS - (latLng[0] + DEG_90_AS_RADS);
 275       if (latLng[1] < 0) {
 276         latLng[1] = latLng[1] + DEG_180_AS_RADS;
 277       } else {
 278         latLng[1] = latLng[1] - DEG_180_AS_RADS;
 279       }
 280     }
 281
 282   }
 283
 284   /**
 285    * Returns a normalized Lng rectangle shape for the bounding box
 286    *
 287    * @param latLng The lat/lon, in radians, lat in position 0, long in position 1
 288    */
 289   public static void normLng(double[] latLng) {
 290     if (latLng[1] > DEG_180_AS_RADS) {
 291       latLng[1] = -1.0 * (DEG_180_AS_RADS - (latLng[1] - DEG_180_AS_RADS));
 292     } else if (latLng[1] < -DEG_180_AS_RADS) {
 293       latLng[1] = (latLng[1] + DEG_180_AS_RADS) + DEG_180_AS_RADS;
 294     }
 295   }
 296
 297   /**
 298    * The square of the Euclidean Distance.  Not really a distance, but useful if all that matters is
 299    * comparing the result to another one.
 300    *
 301    * @param vec1 The first point
 302    * @param vec2 The second point
 303    * @return The squared Euclidean distance
 304    */
 305   public static double squaredEuclideanDistance(double[] vec1, double[] vec2) {
 306     double result = 0;
 307     for (int i = 0; i < vec1.length; i++) {
 308       double v = vec1[i] - vec2[i];
 309       result += v * v;
 310     }
 311     return result;
 312   }
 313
 314   /**
 315    * @param x1     The x coordinate of the first point, in radians
 316    * @param y1     The y coordinate of the first point, in radians
 317    * @param x2     The x coordinate of the second point, in radians
 318    * @param y2     The y coordinate of the second point, in radians
 319    * @param radius The radius of the sphere
 320    * @return The distance between the two points, as determined by the Haversine formula.
 321
 322    */
 323   public static double haversine(double x1, double y1, double x2, double y2, double radius) {
 324     double result = 0;
 325     //make sure they aren't all the same, as then we can just return 0
 326     if ((x1 != x2) || (y1 != y2)) {
 327       double diffX = x1 - x2;
 328       double diffY = y1 - y2;
 329       double hsinX = Math.sin(diffX * 0.5);
 330       double hsinY = Math.sin(diffY * 0.5);
 331       double h = hsinX * hsinX +
 332               (Math.cos(x1) * Math.cos(x2) * hsinY * hsinY);
 333       result = (radius * 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)));
 334     }
 335     return result;
 336   }
 337
 338   /**
 339    * Given a string containing <i>dimension</i> values encoded in it, separated by commas, return a String array of length <i>dimension</i>
 340    * containing the values.
 341    *
 342    * @param out         A preallocated array.  Must be size dimension.  If it is not it will be resized.
 343    * @param externalVal The value to parse
 344    * @param dimension   The expected number of values for the point
 345    * @return An array of the values that make up the point (aka vector)
 346    * @throws org.apache.lucene.spatial.tier.InvalidGeoException if the dimension specified does not match the number of values in the externalValue.
 347    */
 348   public static String[] parsePoint(String[] out, String externalVal, int dimension) throws InvalidGeoException {
 349     //TODO: Should we support sparse vectors?
 350     if (out == null || out.length != dimension) out = new String[dimension];
 351     int idx = externalVal.indexOf(',');
 352     int end = idx;
 353     int start = 0;
 354     int i = 0;
 355     if (idx == -1 && dimension == 1 && externalVal.length() > 0) {//we have a single point, dimension better be 1
 356       out[0] = externalVal.trim();
 357       i = 1;
 358     } else if (idx > 0) {//if it is zero, that is an error
 359       //Parse out a comma separated list of point values, as in: 73.5,89.2,7773.4
 360       for (; i < dimension; i++) {
 361         while (start < end && externalVal.charAt(start) == ' ') start++;
 362         while (end > start && externalVal.charAt(end - 1) == ' ') end--;
 363         if (start == end) {
 364           break;
 365         }
 366         out[i] = externalVal.substring(start, end);
 367         start = idx + 1;
 368         end = externalVal.indexOf(',', start);
 369         idx = end;
 370         if (end == -1) {
 371           end = externalVal.length();
 372         }
 373       }
 374     }
 375     if (i != dimension) {
 376       throw new InvalidGeoException("incompatible dimension (" + dimension +
 377               ") and values (" + externalVal + ").  Only " + i + " values specified");
 378     }
 379     return out;
 380   }
 381
 382   /**
 383    * Given a string containing <i>dimension</i> values encoded in it, separated by commas, return a double array of length <i>dimension</i>
 384    * containing the values.
 385    *
 386    * @param out         A preallocated array.  Must be size dimension.  If it is not it will be resized.
 387    * @param externalVal The value to parse
 388    * @param dimension   The expected number of values for the point
 389    * @return An array of the values that make up the point (aka vector)
 390    * @throws InvalidGeoException if the dimension specified does not match the number of values in the externalValue.
 391    */
 392   public static double[] parsePointDouble(double[] out, String externalVal, int dimension) throws InvalidGeoException{
 393     if (out == null || out.length != dimension) out = new double[dimension];
 394     int idx = externalVal.indexOf(',');
 395     int end = idx;
 396     int start = 0;
 397     int i = 0;
 398     if (idx == -1 && dimension == 1 && externalVal.length() > 0) {//we have a single point, dimension better be 1
 399       out[0] = Double.parseDouble(externalVal.trim());
 400       i = 1;
 401     } else if (idx > 0) {//if it is zero, that is an error
 402       //Parse out a comma separated list of point values, as in: 73.5,89.2,7773.4
 403       for (; i < dimension; i++) {
 404         //TODO: abstract common code with other parsePoint
 405         while (start < end && externalVal.charAt(start) == ' ') start++;
 406         while (end > start && externalVal.charAt(end - 1) == ' ') end--;
 407         if (start == end) {
 408           break;
 409         }
 410         out[i] = Double.parseDouble(externalVal.substring(start, end));
 411         start = idx + 1;
 412         end = externalVal.indexOf(',', start);
 413         idx = end;
 414         if (end == -1) {
 415           end = externalVal.length();
 416         }
 417       }
 418     }
 419     if (i != dimension) {
 420       throw new InvalidGeoException("incompatible dimension (" + dimension +
 421               ") and values (" + externalVal + ").  Only " + i + " values specified");
 422     }
 423     return out;
 424   }
 425
 426   public static final double[] parseLatitudeLongitude(String latLonStr) throws InvalidGeoException {
 427     return parseLatitudeLongitude(null, latLonStr);
 428   }
 429
 430   /**
 431    * extract (by calling {@link #parsePoint(String[], String, int)} and validate the latitude and longitude contained
 432    * in the String by making sure the latitude is between 90 & -90 and longitude is between -180 and 180.
 433    * <p/>
 434    * The latitude is assumed to be the first part of the string and the longitude the second part.
 435    *
 436    * @param latLon    A preallocated array to hold the result
 437    * @param latLonStr The string to parse.  Latitude is the first value, longitude is the second.
 438    * @return The lat long
 439    * @throws InvalidGeoException if there was an error parsing
 440    */
 441   public static final double[] parseLatitudeLongitude(double[] latLon, String latLonStr) throws InvalidGeoException {
 442     if (latLon == null) {
 443       latLon = new double[2];
 444     }
 445     double[] toks = parsePointDouble(null, latLonStr, 2);
 446
 447     if (toks[0] < -90.0 || toks[0] > 90.0) {
 448       throw new InvalidGeoException(
 449               "Invalid latitude: latitudes are range -90 to 90: provided lat: ["
 450                       + toks[0] + "]");
 451     }
 452     latLon[0] = toks[0];
 453
 454
 455     if (toks[1] < -180.0 || toks[1] > 180.0) {
 456
 457       throw new InvalidGeoException(
 458               "Invalid longitude: longitudes are range -180 to 180: provided lon: ["
 459                       + toks[1] + "]");
 460     }
 461     latLon[1] = toks[1];
 462
 463     return latLon;
 464   }
 465 }