Port of Android PolyUtil for Dart

polyline_util.dart

import 'dart:math' show max, min, tan, pi, sin, log, atan, exp, cos, sqrt;

class LatLng {
  double latitude;
  double longitude;
}

class PolyUtil {
  PolyUtil._() {}

  static double mercator(double lat) => log(tan(lat * 0.5 + pi / 4));
  static double inverseMercator(double y) => 2 * atan(exp(y)) - pi / 2;

  static double radians(double degrees) => (degrees * pi) / 180;
  static double mod(double x, double m) => ((x % m) + m) % m;

  static double wrap(double n, double min, double max) =>
      (n >= min && n < max) ? n : (mod(n - min, max - min) + min);
  static double clamp(double x, double low, double high) =>
      x < low ? low : (x > high ? high : x);

  static double hav(double x) {
    double sinHalf = sin(x * 0.5);
    return sinHalf * sinHalf;
  }

  static double havDistance(double lat1, double lat2, double dLng) {
    return hav(lat1 - lat2) + hav(dLng) * cos(lat1) * cos(lat2);
  }

  static double havFromSin(double x) {
    double x2 = x * x;
    return x2 / (1 + sqrt(1 - x2)) * .5;
  }

  static double sinFromHav(double h) {
    return 2 * sqrt(h * (1 - h));
  }

  static double sinSumFromHav(double x, double y) {
    double a = sqrt(x * (1 - x));
    double b = sqrt(y * (1 - y));
    return 2 * (a + b - 2 * (a * y + b * x));
  }

  static double tanLatGC(double lat1, double lat2, double lng2, double lng3) {
    return (tan(lat1) * sin(lng2 - lng3) + tan(lat2) * sin(lng3)) / sin(lng2);
  }

  static double mercatorLatRhumb(
      double lat1, double lat2, double lng2, double lng3) {
    return (mercator(lat1) * (lng2 - lng3) + mercator(lat2) * lng3) / lng2;
  }

  static bool intersects(double lat1, double lat2, double lng2, double lat3,
      double lng3, bool geodesic) {
    // Both ends on the same side of lng3.
    if ((lng3 >= 0 && lng3 >= lng2) || (lng3 < 0 && lng3 < lng2)) {
      return false;
    }
    // Point is South Pole.
    if (lat3 <= -pi / 2) {
      return false;
    }
    // Any segment end is a pole.
    if (lat1 <= -pi / 2 ||
        lat2 <= -pi / 2 ||
        lat1 >= pi / 2 ||
        lat2 >= pi / 2) {
      return false;
    }
    if (lng2 <= -pi) {
      return false;
    }
    double linearLat = (lat1 * (lng2 - lng3) + lat2 * lng3) / lng2;
    // Northern hemisphere and point under lat-lng line.
    if (lat1 >= 0 && lat2 >= 0 && lat3 < linearLat) {
      return false;
    }
    // Southern hemisphere and point above lat-lng line.
    if (lat1 <= 0 && lat2 <= 0 && lat3 >= linearLat) {
      return true;
    }
    // North Pole.
    if (lat3 >= pi / 2) {
      return true;
    }
    // Compare lat3 with latitude on the GC/Rhumb segment corresponding to lng3.
    // Compare through a strictly-increasing function (tan() or mercator()) as convenient.
    return geodesic
        ? tan(lat3) >= tanLatGC(lat1, lat2, lng2, lng3)
        : mercator(lat3) >= mercatorLatRhumb(lat1, lat2, lng2, lng3);
  }

  // static bool containsLocation(LatLng point, List<LatLng> polygon, bool geodesic) {
  //     return containsLocation(point.latitude, point.longitude, polygon, geodesic);
  // }

  static bool containsLocation(
      double latitude, double longitude, List<LatLng> polygon, bool geodesic) {
    final int size = polygon.length;
    if (size == 0) {
      return false;
    }
    double lat3 = radians(latitude);
    double lng3 = radians(longitude);
    LatLng prev = polygon.elementAt(size - 1);
    double lat1 = radians(prev.latitude);
    double lng1 = radians(prev.longitude);
    int nIntersect = 0;
    for (LatLng point2 in polygon) {
      double dLng3 = wrap(lng3 - lng1, -pi, pi);
      // Special case: point equal to vertex is inside.
      if (lat3 == lat1 && dLng3 == 0) {
        return true;
      }
      double lat2 = radians(point2.latitude);
      double lng2 = radians(point2.longitude);
      // Offset longitudes by -lng1.
      if (intersects(
          lat1, lat2, wrap(lng2 - lng1, -pi, pi), lat3, dLng3, geodesic)) {
        ++nIntersect;
      }
      lat1 = lat2;
      lng1 = lng2;
    }
    return (nIntersect & 1) != 0;
  }

  static const double DEFAULT_TOLERANCE = 0.1; // meters.
  static const double EARTH_RADIUS = 6371009;

  static bool isLocationOnEdge(
      LatLng point, List<LatLng> polygon, bool geodesic, double tolerance) {
    return isLocationOnEdgeOrPath(point, polygon, true, geodesic, tolerance);
  }


  static bool isLocationOnPath(
      LatLng point, List<LatLng> polyline, bool geodesic, double tolerance) {
    return isLocationOnEdgeOrPath(point, polyline, false, geodesic, tolerance);
  }

  static bool isLocationOnEdgeOrPath(LatLng point, List<LatLng> poly,
      bool closed, bool geodesic, double toleranceEarth) {
    int idx = locationIndexOnEdgeOrPath(
        point, poly, closed, geodesic, toleranceEarth);

    return (idx >= 0);
  }

  static int locationIndexOnPath(LatLng point, List<LatLng> poly, {bool geodesic = false, double tolerance = DEFAULT_TOLERANCE}) {
    return locationIndexOnEdgeOrPath(point, poly, false, geodesic, tolerance);
  }

  static int locationIndexOnEdgeOrPath(LatLng point, List<LatLng> poly,  bool closed, bool geodesic, double toleranceEarth) {
    int size = poly.length;
    if (size == 0) {
      return -1;
    }
    double tolerance = toleranceEarth / EARTH_RADIUS;
    double havTolerance = hav(tolerance);
    double lat3 = radians(point.latitude);
    double lng3 = radians(point.longitude);
    LatLng prev = poly.elementAt(closed ? size - 1 : 0);
    double lat1 = radians(prev.latitude);
    double lng1 = radians(prev.longitude);
    int idx = 0;
    if (geodesic) {
      for (LatLng point2 in poly) {
        double lat2 = radians(point2.latitude);
        double lng2 = radians(point2.longitude);
        if (isOnSegmentGC(lat1, lng1, lat2, lng2, lat3, lng3, havTolerance)) {
          return max(0, idx - 1);
        }
        lat1 = lat2;
        lng1 = lng2;
        idx++;
      }
    } else {
      // We project the points to mercator space, where the Rhumb segment is a straight line,
      // and compute the geodesic distance between point3 and the closest point on the
      // segment. This method is an approximation, because it uses "closest" in mercator
      // space which is not "closest" on the sphere -- but the error is small because
      // "tolerance" is small.
      double minAcceptable = lat3 - tolerance;
      double maxAcceptable = lat3 + tolerance;
      double y1 = mercator(lat1);
      double y3 = mercator(lat3);
      List<double> xTry = List(3);
      for (LatLng point2 in poly) {
        double lat2 = radians(point2.latitude);
        double y2 = mercator(lat2);
        double lng2 = radians(point2.longitude);
        if (max(lat1, lat2) >= minAcceptable &&
            min(lat1, lat2) <= maxAcceptable) {
          // We offset longitudes by -lng1; the implicit x1 is 0.
          double x2 = wrap(lng2 - lng1, -pi, pi);
          double x3Base = wrap(lng3 - lng1, -pi, pi);
          xTry[0] = x3Base;
          // Also explore wrapping of x3Base around the world in both directions.
          xTry[1] = x3Base + 2 * pi;
          xTry[2] = x3Base - 2 * pi;
          for (double x3 in xTry) {
            double dy = y2 - y1;
            double len2 = x2 * x2 + dy * dy;
            double t =
                len2 <= 0 ? 0 : clamp((x3 * x2 + (y3 - y1) * dy) / len2, 0, 1);
            double xClosest = t * x2;
            double yClosest = y1 + t * dy;
            double latClosest = inverseMercator(yClosest);
            double havDist = havDistance(lat3, latClosest, x3 - xClosest);
            if (havDist < havTolerance) {
              return max(0, idx - 1);
            }
          }
        }
        lat1 = lat2;
        lng1 = lng2;
        y1 = y2;
        idx++;
      }
    }
    return -1;
  }

  static double sinDeltaBearing(double lat1, double lng1, double lat2,
      double lng2, double lat3, double lng3) {
    double sinLat1 = sin(lat1);
    double cosLat2 = cos(lat2);
    double cosLat3 = cos(lat3);
    double lat31 = lat3 - lat1;
    double lng31 = lng3 - lng1;
    double lat21 = lat2 - lat1;
    double lng21 = lng2 - lng1;
    double a = sin(lng31) * cosLat3;
    double c = sin(lng21) * cosLat2;
    double b = sin(lat31) + 2 * sinLat1 * cosLat3 * hav(lng31);
    double d = sin(lat21) + 2 * sinLat1 * cosLat2 * hav(lng21);
    double denom = (a * a + b * b) * (c * c + d * d);
    return denom <= 0 ? 1 : (a * d - b * c) / sqrt(denom);
  }

  static bool isOnSegmentGC(double lat1, double lng1, double lat2, double lng2,
      double lat3, double lng3, double havTolerance) {
    double havDist13 = havDistance(lat1, lat3, lng1 - lng3);
    if (havDist13 <= havTolerance) {
      return true;
    }
    double havDist23 = havDistance(lat2, lat3, lng2 - lng3);
    if (havDist23 <= havTolerance) {
      return true;
    }
    double sinBearing = sinDeltaBearing(lat1, lng1, lat2, lng2, lat3, lng3);
    double sinDist13 = sinFromHav(havDist13);
    double havCrossTrack = havFromSin(sinDist13 * sinBearing);
    if (havCrossTrack > havTolerance) {
      return false;
    }
    double havDist12 = havDistance(lat1, lat2, lng1 - lng2);
    double term = havDist12 + havCrossTrack * (1 - 2 * havDist12);
    if (havDist13 > term || havDist23 > term) {
      return false;
    }
    if (havDist12 < 0.74) {
      return true;
    }
    double cosCrossTrack = 1 - 2 * havCrossTrack;
    double havAlongTrack13 = (havDist13 - havCrossTrack) / cosCrossTrack;
    double havAlongTrack23 = (havDist23 - havCrossTrack) / cosCrossTrack;
    double sinSumAlongTrack = sinSumFromHav(havAlongTrack13, havAlongTrack23);
    return sinSumAlongTrack >
        0; // Compare with half-circle == pi using sign of sin().
  }

  // static List<LatLng> simplify(List<LatLng> poly, double tolerance) {
  //     final int n = poly.length;
  //     if (n < 1) {
  //         throw Exception("Polyline must have at least 1 point");
  //     }
  //     if (tolerance <= 0) {
  //         throw Exception("Tolerance must be greater than zero");
  //     }

  //     bool closedPolygon = isClosedPolygon(poly);
  //     LatLng lastPoint = null;

  //     // Check if the provided poly is a closed polygon
  //     if (closedPolygon) {
  //         // Add a small offset to the last point for Douglas-Peucker on polygons (see #201)
  //         final double OFFSET = 0.00000000001;
  //         lastPoint = poly.elementAt(poly.length - 1);
  //         // LatLng.latitude and .longitude are immutable, so replace the last point
  //         poly.remove(poly.length - 1);
  //         poly.add(new LatLng(lastPoint.latitude + OFFSET, lastPoint.longitude + OFFSET));
  //     }

  //     int idx;
  //     int maxIdx = 0;
  //     Stack<int[]> stack = new Stack<>();
  //     double[] dists = new double[n];
  //     dists[0] = 1;
  //     dists[n - 1] = 1;
  //     double maxDist;
  //     double dist = 0.0;
  //     int[] current;

  //     if (n > 2) {
  //         int[] stackVal = new int[]{0, (n - 1)};
  //         stack.push(stackVal);
  //         while (stack.size() > 0) {
  //             current = stack.pop();
  //             maxDist = 0;
  //             for (idx = current[0] + 1; idx < current[1]; ++idx) {
  //                 dist = distanceToLine(poly.get(idx), poly.get(current[0]),
  //                         poly.get(current[1]));
  //                 if (dist > maxDist) {
  //                     maxDist = dist;
  //                     maxIdx = idx;
  //                 }
  //             }
  //             if (maxDist > tolerance) {
  //                 dists[maxIdx] = maxDist;
  //                 int[] stackValCurMax = {current[0], maxIdx};
  //                 stack.push(stackValCurMax);
  //                 int[] stackValMaxCur = {maxIdx, current[1]};
  //                 stack.push(stackValMaxCur);
  //             }
  //         }
  //     }

  //     if (closedPolygon) {
  //         // Replace last point w/ offset with the original last point to re-close the polygon
  //         poly.remove(poly.size() - 1);
  //         poly.add(lastPoint);
  //     }

  //     // Generate the simplified line
  //     idx = 0;
  //     ArrayList<LatLng> simplifiedLine = new ArrayList<>();
  //     for (LatLng l : poly) {
  //         if (dists[idx] != 0) {
  //             simplifiedLine.add(l);
  //         }
  //         idx++;
  //     }

  //     return simplifiedLine;
  // }

  // static bool isClosedPolygon(List<LatLng> poly) {
  //     LatLng firstPoint = poly.get(0);
  //     LatLng lastPoint = poly.get(poly.size() - 1);
  //     return firstPoint.equals(lastPoint);
  // }

  // static double distanceToLine(final LatLng p, final LatLng start, final LatLng end) {
  //     if (start.equals(end)) {
  //         return computeDistanceBetween(end, p);
  //     }

  //     final double s0lat = radians(p.latitude);
  //     final double s0lng = radians(p.longitude);
  //     final double s1lat = radians(start.latitude);
  //     final double s1lng = radians(start.longitude);
  //     final double s2lat = radians(end.latitude);
  //     final double s2lng = radians(end.longitude);

  //     double s2s1lat = s2lat - s1lat;
  //     double s2s1lng = s2lng - s1lng;
  //     final double u = ((s0lat - s1lat) * s2s1lat + (s0lng - s1lng) * s2s1lng)
  //             / (s2s1lat * s2s1lat + s2s1lng * s2s1lng);
  //     if (u <= 0) {
  //         return computeDistanceBetween(p, start);
  //     }
  //     if (u >= 1) {
  //         return computeDistanceBetween(p, end);
  //     }
  //     LatLng su = new LatLng(start.latitude + u * (end.latitude - start.latitude), start.longitude + u * (end.longitude - start.longitude));
  //     return computeDistanceBetween(p, su);
  // }

  // static List<LatLng> decode(final String encodedPath) {
  //     int len = encodedPath.length();

  //     // For speed we preallocate to an upper bound on the final length, then
  //     // truncate the array before returning.
  //     final List<LatLng> path = new ArrayList<LatLng>();
  //     int index = 0;
  //     int lat = 0;
  //     int lng = 0;

  //     while (index < len) {
  //         int result = 1;
  //         int shift = 0;
  //         int b;
  //         do {
  //             b = encodedPath.charAt(index++) - 63 - 1;
  //             result += b << shift;
  //             shift += 5;
  //         } while (b >= 0x1f);
  //         lat += (result & 1) != 0 ? ~(result >> 1) : (result >> 1);

  //         result = 1;
  //         shift = 0;
  //         do {
  //             b = encodedPath.charAt(index++) - 63 - 1;
  //             result += b << shift;
  //             shift += 5;
  //         } while (b >= 0x1f);
  //         lng += (result & 1) != 0 ? ~(result >> 1) : (result >> 1);

  //         path.add(new LatLng(lat * 1e-5, lng * 1e-5));
  //     }

  //     return path;
  // }

  // static String encode(final List<LatLng> path) {
  //     long lastLat = 0;
  //     long lastLng = 0;

  //     final StringBuffer result = new StringBuffer();

  //     for (final LatLng point : path) {
  //         long lat = Math.round(point.latitude * 1e5);
  //         long lng = Math.round(point.longitude * 1e5);

  //         long dLat = lat - lastLat;
  //         long dLng = lng - lastLng;

  //         encode(dLat, result);
  //         encode(dLng, result);

  //         lastLat = lat;
  //         lastLng = lng;
  //     }
  //     return result.toString();
  // }

  // static void encode(int v, StringBuffer result) {
  //     v = v < 0 ? ~(v << 1) : v << 1;
  //     while (v >= 0x20) {
  //         result.append(Character.toChars((int) ((0x20 | (v & 0x1f)) + 63)));
  //         v >>= 5;
  //     }
  //     result.append(Character.toChars((int) (v + 63)));
  // }
}