package org.apache.lucene.util;

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

import java.util.ArrayList;

/**
 * Basic reusable geo-spatial utility methods
 *
 * @lucene.experimental
 */
public final class GeoUtils {
  public static final short BITS = 31;
  private static final double LON_SCALE = (0x1L<<BITS)/360.0D;
  private static final double LAT_SCALE = (0x1L<<BITS)/180.0D;
  public static final double TOLERANCE = 1E-6;

  /** Minimum longitude value. */
  public static final double MIN_LON_INCL = -180.0D;

  /** Maximum longitude value. */
  public static final double MAX_LON_INCL = 180.0D;

  /** Minimum latitude value. */
  public static final double MIN_LAT_INCL = -90.0D;

  /** Maximum latitude value. */
  public static final double MAX_LAT_INCL = 90.0D;

  // No instance:
  private GeoUtils() {
  }

  public static final Long mortonHash(final double lon, final double lat) {
    return BitUtil.interleave(scaleLon(lon), scaleLat(lat));
  }

  public static final double mortonUnhashLon(final long hash) {
    return unscaleLon(BitUtil.deinterleave(hash));
  }

  public static final double mortonUnhashLat(final long hash) {
    return unscaleLat(BitUtil.deinterleave(hash >>> 1));
  }

  private static final long scaleLon(final double val) {
    return (long) ((val-MIN_LON_INCL) * LON_SCALE);
  }

  private static final long scaleLat(final double val) {
    return (long) ((val-MIN_LAT_INCL) * LAT_SCALE);
  }

  private static final double unscaleLon(final long val) {
    return (val / LON_SCALE) + MIN_LON_INCL;
  }

  private static final double unscaleLat(final long val) {
    return (val / LAT_SCALE) + MIN_LAT_INCL;
  }

  public static double compare(final double v1, final double v2) {
    final double delta = v1-v2;
    return Math.abs(delta) <= TOLERANCE ? 0 : delta;
  }

  /**
   * Puts longitude in range of -180 to +180.
   */
  public static double normalizeLon(double lon_deg) {
    if (lon_deg >= -180 && lon_deg <= 180) {
      return lon_deg; //common case, and avoids slight double precision shifting
    }
    double off = (lon_deg + 180) % 360;
    if (off < 0) {
      return 180 + off;
    } else if (off == 0 && lon_deg > 0) {
      return 180;
    } else {
      return -180 + off;
    }
  }

  /**
   * Puts latitude in range of -90 to 90.
   */
  public static double normalizeLat(double lat_deg) {
    if (lat_deg >= -90 && lat_deg <= 90) {
      return lat_deg; //common case, and avoids slight double precision shifting
    }
    double off = Math.abs((lat_deg + 90) % 360);
    return (off <= 180 ? off : 360-off) - 90;
  }

  /**
   * Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
   * NOTE: this is a basic method that does not handle dateline or pole crossing. Unwrapping must be done before
   * calling this method.
   */
  public static boolean bboxContains(final double lon, final double lat, final double minLon,
                                     final double minLat, final double maxLon, final double maxLat) {
    return (compare(lon, minLon) >= 0 && compare(lon, maxLon) <= 0
          && compare(lat, minLat) >= 0 && compare(lat, maxLat) <= 0);
  }

  /**
   * simple even-odd point in polygon computation
   *    1.  Determine if point is contained in the longitudinal range
   *    2.  Determine whether point crosses the edge by computing the latitudinal delta
   *        between the end-point of a parallel vector (originating at the point) and the
   *        y-component of the edge sink
   *
   * NOTE: Requires polygon point (x,y) order either clockwise or counter-clockwise
   */
  public static boolean pointInPolygon(double[] x, double[] y, double lat, double lon) {
    assert x.length == y.length;
    boolean inPoly = false;
    /**
     * Note: This is using a euclidean coordinate system which could result in
     * upwards of 110KM error at the equator.
     * TODO convert coordinates to cylindrical projection (e.g. mercator)
     */
    for (int i = 1; i < x.length; i++) {
      if (x[i] < lon && x[i-1] >= lon || x[i-1] < lon && x[i] >= lon) {
        if (y[i] + (lon - x[i]) / (x[i-1] - x[i]) * (y[i-1] - y[i]) < lat) {
          inPoly = !inPoly;
        }
      }
    }
    return inPoly;
  }

  public static String geoTermToString(long term) {
    StringBuilder s = new StringBuilder(64);
    final int numberOfLeadingZeros = Long.numberOfLeadingZeros(term);
    for (int i = 0; i < numberOfLeadingZeros; i++) {
      s.append('0');
    }
    if (term != 0) {
      s.append(Long.toBinaryString(term));
    }
    return s.toString();
  }


  public static boolean rectDisjoint(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
                                     final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
    return (aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY);
  }

  /**
   * Computes whether the first (a) rectangle is wholly within another (b) rectangle (shared boundaries allowed)
   */
  public static boolean rectWithin(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
                                   final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
    return !(aMinX < bMinX || aMinY < bMinY || aMaxX > bMaxX || aMaxY > bMaxY);
  }

  public static boolean rectCrosses(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
                                    final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
    return !(rectDisjoint(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY) ||
        rectWithin(aMinX, aMinY, aMaxX, aMaxY, bMinX, bMinY, bMaxX, bMaxY));
  }

  /**
   * Computes whether rectangle a contains rectangle b (touching allowed)
   */
  public static boolean rectContains(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
                                     final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
    return !(bMinX < aMinX || bMinY < aMinY || bMaxX > aMaxX || bMaxY > aMaxY);
  }

  /**
   * Computes whether a rectangle intersects another rectangle (crosses, within, touching, etc)
   */
  public static boolean rectIntersects(final double aMinX, final double aMinY, final double aMaxX, final double aMaxY,
                                       final double bMinX, final double bMinY, final double bMaxX, final double bMaxY) {
    return !((aMaxX < bMinX || aMinX > bMaxX || aMaxY < bMinY || aMinY > bMaxY) );
  }

  /**
   * Computes whether a rectangle crosses a shape. (touching not allowed)
   */
  public static boolean rectCrossesPoly(final double rMinX, final double rMinY, final double rMaxX,
                                        final double rMaxY, final double[] shapeX, final double[] shapeY,
                                        final double sMinX, final double sMinY, final double sMaxX,
                                        final double sMaxY) {
    // short-circuit: if the bounding boxes are disjoint then the shape does not cross
    if (rectDisjoint(rMinX, rMinY, rMaxX, rMaxY, sMinX, sMinY, sMaxX, sMaxY)) {
      return false;
    }

    final double[][] bbox = new double[][] { {rMinX, rMinY}, {rMaxX, rMinY}, {rMaxX, rMaxY}, {rMinX, rMaxY}, {rMinX, rMinY} };
    final int polyLength = shapeX.length-1;
    double d, s, t, a1, b1, c1, a2, b2, c2;
    double x00, y00, x01, y01, x10, y10, x11, y11;

    // computes the intersection point between each bbox edge and the polygon edge
    for (short b=0; b<4; ++b) {
      a1 = bbox[b+1][1]-bbox[b][1];
      b1 = bbox[b][0]-bbox[b+1][0];
      c1 = a1*bbox[b+1][0] + b1*bbox[b+1][1];
      for (int p=0; p<polyLength; ++p) {
        a2 = shapeY[p+1]-shapeY[p];
        b2 = shapeX[p]-shapeX[p+1];
        // compute determinant
        d = a1*b2 - a2*b1;
        if (d != 0) {
          // lines are not parallel, check intersecting points
          c2 = a2*shapeX[p+1] + b2*shapeY[p+1];
          s = (1/d)*(b2*c1 - b1*c2);
          t = (1/d)*(a1*c2 - a2*c1);
          x00 = StrictMath.min(bbox[b][0], bbox[b+1][0]) - TOLERANCE;
          x01 = StrictMath.max(bbox[b][0], bbox[b+1][0]) + TOLERANCE;
          y00 = StrictMath.min(bbox[b][1], bbox[b+1][1]) - TOLERANCE;
          y01 = StrictMath.max(bbox[b][1], bbox[b+1][1]) + TOLERANCE;
          x10 = StrictMath.min(shapeX[p], shapeX[p+1]) - TOLERANCE;
          x11 = StrictMath.max(shapeX[p], shapeX[p+1]) + TOLERANCE;
          y10 = StrictMath.min(shapeY[p], shapeY[p+1]) - TOLERANCE;
          y11 = StrictMath.max(shapeY[p], shapeY[p+1]) + TOLERANCE;
          // check whether the intersection point is touching one of the line segments
          boolean touching = ((x00 == s && y00 == t) || (x01 == s && y01 == t))
              || ((x10 == s && y10 == t) || (x11 == s && y11 == t));
          // if line segments are not touching and the intersection point is within the range of either segment
          if (!(touching || x00 > s || x01 < s || y00 > t || y01 < t || x10 > s || x11 < s || y10 > t || y11 < t)) {
            return true;
          }
        }
      } // for each poly edge
    } // for each bbox edge
    return false;
  }

  public static boolean lineCrossesPoly(double x1, double y1, double x2, double y2, double[] shapeX, double[] shapeY, final double sMinX, final double sMinY, final double sMaxX,
                                        final double sMaxY) {
    final double a1 = y2 - y1;
    final double b1 = x2 - x1;
    final double c1 = a1*x2 + b1*y2;
    final int polyLength = shapeX.length;
    double a2, b2, c2, s, t, d;
    double x00, x01, y00, y01, x10, x11, y10, y11;
    for (int p=0; p<polyLength; ++p) {
      a2 = shapeY[p+1]-shapeY[p];
      b2 = shapeX[p]-shapeX[p+1];
      // compute determinant
      d = a1*b2 - a2*b1;
      if (d != 0) {
        // lines are not parallel, check intersecting points
        c2 = a2*shapeX[p+1] + b2*shapeY[p+1];
        s = (1/d)*(b2*c1 - b1*c2);
        t = (1/d)*(a1*c2 - a2*c1);
        x00 = StrictMath.min(x1, x2) - TOLERANCE;
        x01 = StrictMath.max(x1, x2) + TOLERANCE;
        y00 = StrictMath.min(y1, y2) - TOLERANCE;
        y01 = StrictMath.max(y1, y2) + TOLERANCE;
        x10 = StrictMath.min(shapeX[p], shapeX[p+1]) - TOLERANCE;
        x11 = StrictMath.max(shapeX[p], shapeX[p+1]) + TOLERANCE;
        y10 = StrictMath.min(shapeY[p], shapeY[p+1]) - TOLERANCE;
        y11 = StrictMath.max(shapeY[p], shapeY[p+1]) + TOLERANCE;
        // check whether the intersection point is touching one of the line segments
        boolean touching = ((x00 == s && y00 == t) || (x01 == s && y01 == t))
            || ((x10 == s && y10 == t) || (x11 == s && y11 == t));
        // if line segments are not touching and the intersection point is within the range of either segment
        if (!(touching || x00 > s || x01 < s || y00 > t || y01 < t || x10 > s || x11 < s || y10 > t || y11 < t)) {
          return true;
        }
      }
    } // for each poly edge
    return false;
  }

  /**
   * Converts a given circle (defined as a point/radius) to an approximated line-segment polygon
   *
   * @param lon longitudinal center of circle (in degrees)
   * @param lat latitudinal center of circle (in degrees)
   * @param radiusMeters distance radius of circle (in meters)
   * @return a list of lon/lat points representing the circle
   */
  @SuppressWarnings({"unchecked","rawtypes"})
  public static ArrayList<double[]> circleToPoly(final double lon, final double lat, final double radiusMeters) {
    double angle;
    // a little under-sampling (to limit the number of polygonal points): using archimedes estimation of pi
    final int sides = 25;
    ArrayList<double[]> geometry = new ArrayList();
    double[] lons = new double[sides];
    double[] lats = new double[sides];

    double[] pt = new double[2];
    final int sidesLen = sides-1;
    for (int i=0; i<sidesLen; ++i) {
      angle = (i*360/sides);
      pt = GeoProjectionUtils.pointFromLonLatBearing(lon, lat, angle, radiusMeters, pt);
      lons[i] = pt[0];
      lats[i] = pt[1];
    }
    // close the poly
    lons[sidesLen] = lons[0];
    lats[sidesLen] = lats[0];
    geometry.add(lons);
    geometry.add(lats);

    return geometry;
  }

  /**
   * Computes whether a rectangle is within a given polygon (shared boundaries allowed)
   */
  public static boolean rectWithinPoly(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
                                       final double[] shapeX, final double[] shapeY, final double sMinX,
                                       final double sMinY, final double sMaxX, final double sMaxY) {
    // check if rectangle crosses poly (to handle concave/pacman polys), then check that all 4 corners
    // are contained
    return !(rectCrossesPoly(rMinX, rMinY, rMaxX, rMaxY, shapeX, shapeY, sMinX, sMinY, sMaxX, sMaxY) ||
        !pointInPolygon(shapeX, shapeY, rMinY, rMinX) || !pointInPolygon(shapeX, shapeY, rMinY, rMaxX) ||
        !pointInPolygon(shapeX, shapeY, rMaxY, rMaxX) || !pointInPolygon(shapeX, shapeY, rMaxY, rMinX));
  }

  private static boolean rectAnyCornersOutsideCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
                                                     final double centerLon, final double centerLat, final double radiusMeters) {
    return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 > radiusMeters
      || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 > radiusMeters
      || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 > radiusMeters
      || SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 > radiusMeters;
  }

  private static boolean rectAnyCornersInCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
                                                final double centerLon, final double centerLat, final double radiusMeters) {
    return SloppyMath.haversin(centerLat, centerLon, rMinY, rMinX)*1000.0 <= radiusMeters
        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMinX)*1000.0 <= radiusMeters
        || SloppyMath.haversin(centerLat, centerLon, rMaxY, rMaxX)*1000.0 <= radiusMeters
        || SloppyMath.haversin(centerLat, centerLon, rMinY, rMaxX)*1000.0 <= radiusMeters;
  }

  public static boolean rectWithinCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
                                         final double centerLon, final double centerLat, final double radiusMeters) {
    return rectAnyCornersOutsideCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters) == false;
  }

  /**
   * Determine if a bbox (defined by minLon, minLat, maxLon, maxLat) contains the provided point (defined by lon, lat)
   * NOTE: this is basic method that does not handle dateline or pole crossing. Unwrapping must be done before
   * calling this method.
   */
  public static boolean rectCrossesCircle(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
                                          final double centerLon, final double centerLat, final double radiusMeters) {
    return rectAnyCornersInCircle(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters)
        || isClosestPointOnRectWithinRange(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, radiusMeters);
  }

  private static boolean isClosestPointOnRectWithinRange(final double rMinX, final double rMinY, final double rMaxX, final double rMaxY,
                                                         final double centerLon, final double centerLat, final double radiusMeters) {
    double[] closestPt = {0, 0};
    GeoDistanceUtils.closestPointOnBBox(rMinX, rMinY, rMaxX, rMaxY, centerLon, centerLat, closestPt);
    return SloppyMath.haversin(centerLat, centerLon, closestPt[1], closestPt[0])*1000.0 <= radiusMeters;
  }

  /**
   * Compute Bounding Box for a circle using WGS-84 parameters
   */
  public static GeoRect circleToBBox(final double centerLon, final double centerLat, final double radiusMeters) {
    final double radLat = StrictMath.toRadians(centerLat);
    final double radLon = StrictMath.toRadians(centerLon);
    double radDistance = radiusMeters / GeoProjectionUtils.SEMIMAJOR_AXIS;
    double minLat = radLat - radDistance;
    double maxLat = radLat + radDistance;
    double minLon;
    double maxLon;

    if (minLat > GeoProjectionUtils.MIN_LAT_RADIANS && maxLat < GeoProjectionUtils.MAX_LAT_RADIANS) {
      double deltaLon = StrictMath.asin(StrictMath.sin(radDistance) / StrictMath.cos(radLat));
      minLon = radLon - deltaLon;
      if (minLon < GeoProjectionUtils.MIN_LON_RADIANS) {
        minLon += 2d * StrictMath.PI;
      }
      maxLon = radLon + deltaLon;
      if (maxLon > GeoProjectionUtils.MAX_LON_RADIANS) {
        maxLon -= 2d * StrictMath.PI;
      }
    } else {
      // a pole is within the distance
      minLat = StrictMath.max(minLat, GeoProjectionUtils.MIN_LAT_RADIANS);
      maxLat = StrictMath.min(maxLat, GeoProjectionUtils.MAX_LAT_RADIANS);
      minLon = GeoProjectionUtils.MIN_LON_RADIANS;
      maxLon = GeoProjectionUtils.MAX_LON_RADIANS;
    }

    return new GeoRect(StrictMath.toDegrees(minLon), StrictMath.toDegrees(maxLon),
        StrictMath.toDegrees(minLat), StrictMath.toDegrees(maxLat));
  }

  public static GeoRect polyToBBox(double[] polyLons, double[] polyLats) {
    if (polyLons.length != polyLats.length) {
      throw new IllegalArgumentException("polyLons and polyLats must be equal length");
    }

    double minLon = Double.POSITIVE_INFINITY;
    double maxLon = Double.NEGATIVE_INFINITY;
    double minLat = Double.POSITIVE_INFINITY;
    double maxLat = Double.NEGATIVE_INFINITY;

    for (int i=0;i<polyLats.length;i++) {
      if (GeoUtils.isValidLon(polyLons[i]) == false) {
        throw new IllegalArgumentException("invalid polyLons[" + i + "]=" + polyLons[i]);
      }
      if (GeoUtils.isValidLat(polyLats[i]) == false) {
        throw new IllegalArgumentException("invalid polyLats[" + i + "]=" + polyLats[i]);
      }
      minLon = Math.min(polyLons[i], minLon);
      maxLon = Math.max(polyLons[i], maxLon);
      minLat = Math.min(polyLats[i], minLat);
      maxLat = Math.max(polyLats[i], maxLat);
    }

    return new GeoRect(GeoUtils.unscaleLon(GeoUtils.scaleLon(minLon)), GeoUtils.unscaleLon(GeoUtils.scaleLon(maxLon)),
        GeoUtils.unscaleLat(GeoUtils.scaleLat(minLat)), GeoUtils.unscaleLat(GeoUtils.scaleLat(maxLat)));
  }

  /*
  /**
   * Computes whether or a 3dimensional line segment intersects or crosses a sphere
   *
   * @param lon1 longitudinal location of the line segment start point (in degrees)
   * @param lat1 latitudinal location of the line segment start point (in degrees)
   * @param alt1 altitude of the line segment start point (in degrees)
   * @param lon2 longitudinal location of the line segment end point (in degrees)
   * @param lat2 latitudinal location of the line segment end point (in degrees)
   * @param alt2 altitude of the line segment end point (in degrees)
   * @param centerLon longitudinal location of center search point (in degrees)
   * @param centerLat latitudinal location of center search point (in degrees)
   * @param centerAlt altitude of the center point (in meters)
   * @param radiusMeters search sphere radius (in meters)
   * @return whether the provided line segment is a secant of the
   * /
  // NOTE: not used for 2d at the moment. used for 3d w/ altitude (we can keep or add back)
  private static boolean lineCrossesSphere(double lon1, double lat1, double alt1, double lon2,
                                           double lat2, double alt2, double centerLon, double centerLat,
                                           double centerAlt, double radiusMeters) {
    // convert to cartesian 3d (in meters)
    double[] ecf1 = GeoProjectionUtils.llaToECF(lon1, lat1, alt1, null);
    double[] ecf2 = GeoProjectionUtils.llaToECF(lon2, lat2, alt2, null);
    double[] cntr = GeoProjectionUtils.llaToECF(centerLon, centerLat, centerAlt, null);

    final double dX = ecf2[0] - ecf1[0];
    final double dY = ecf2[1] - ecf1[1];
    final double dZ = ecf2[2] - ecf1[2];
    final double fX = ecf1[0] - cntr[0];
    final double fY = ecf1[1] - cntr[1];
    final double fZ = ecf1[2] - cntr[2];

    final double a = dX*dX + dY*dY + dZ*dZ;
    final double b = 2 * (fX*dX + fY*dY + fZ*dZ);
    final double c = (fX*fX + fY*fY + fZ*fZ) - (radiusMeters*radiusMeters);

    double discrim = (b*b)-(4*a*c);
    if (discrim < 0) {
      return false;
    }

    discrim = StrictMath.sqrt(discrim);
    final double a2 = 2*a;
    final double t1 = (-b - discrim)/a2;
    final double t2 = (-b + discrim)/a2;

    if ( (t1 < 0 || t1 > 1) ) {
      return !(t2 < 0 || t2 > 1);
    }

    return true;
  }
  */

  public static boolean isValidLat(double lat) {
    return Double.isNaN(lat) == false && lat >= MIN_LAT_INCL && lat <= MAX_LAT_INCL;
  }

  public static boolean isValidLon(double lon) {
    return Double.isNaN(lon) == false && lon >= MIN_LON_INCL && lon <= MAX_LON_INCL;
  }
}
