intersection.java

openmap java写的开源数字地图程序. 用applet实现,可以像google map 那样放大缩小地图.

     * <code>poly<code> is an array of latitude/longitude points where:
     * <br>
     * <pre>
     *                                  
     *                                                                         
     *                                                                          poly[0] = latitude 1
     *                                                                          poly[1] = longitude 1
     *                                                                          poly[2] = latitude 2
     *                                                                          poly[3] = longitude 2
     *                                                                          .
     *                                                                          .
     *                                                                          .
     *                                                                          poly[n-1] = latitude 1
     *                                                                          poly[n] = longitude 1
     *                                                                          
     * </pre>
     *
     * @param x a geographic coordinate
     * @param poly an array of lat/lons describing a closed polygon
     * @return true iff <code>x</code> or <code>antipode(x)</code> is
     * inside <code>poly</code>
     */
    public static boolean isPointInPolygon(Geo x, Geo[] poly) {
        Geo c = center(poly);

        // bail out if the point is more than 90 degrees off the
        // centroid
        double d = x.distance(c);
        if (d >= (Math.PI / 2)) {
            return false;
        }
        // ray is normal to the great circle from c to x.
        Geo ray = c.crossNormalize(x);
        /*
         * side is a point on the great circle between c and x. It is
         * used to choose a direction.
         */
        Geo side = Geo.crossNormalize(x, ray);
        boolean in = false;
        Geo p1 = new Geo(poly[0]);
        Geo p2 = new Geo(poly[0]);
        for (int i = 1; i < poly.length; i++) {
            p2.initialize(poly[i]);
            /*
             * p1 and p2 are on different sides of the ray, and the
             * great acircle between p1 and p2 is on the side that
             * counts;
             */
            if ((p1.dot(ray) < 0.0) != (p2.dot(ray) < 0.0)
                    && p1.intersect(p2, ray).dot(side) > 0.0)
                in = !in;
            Geo temp = p1;
            p1 = p2;
            p2 = temp;
        }

        // Check for unclosed polygons, if the polygon isn't closed,
        // do the calculation for the last point to the starting
        // point.
        if (!p1.equals(poly[0])) {
            p2.initialize(poly[0]);
            if ((p1.dot(ray) < 0.0) != (p2.dot(ray) < 0.0)
                    && p1.intersect(p2, ray).dot(side) > 0.0) {
                in = !in;
            }
        }

        return in;
    }

    /**
     * Ask if a Geo point is in a polygon, with the poly coordinates
     * specified in radians.
     * 
     * @param x
     * @param poly float array where [lat, lon, lat, lon,...] are in
     *        radians
     * @return true for Geo in poly
     */
    public static boolean isPointInPolygonRadians(Geo x, float[] poly) {
        return isPointInPolygon(x, poly, false);
    }

    /**
     * Ask if a Geo point is in a polygon, with the poly coordinates
     * specified in decimal degrees.
     * 
     * @param x
     * @param poly float array where [lat, lon, lat, lon,...] are in
     *        decimal degrees
     * @return true for Geo in poly
     */
    public static boolean isPointInPolygon(Geo x, float[] poly) {
        return isPointInPolygon(x, poly, true);
    }

    /**
     * Ask if a Geo point is in a polygon.
     * 
     * @param x
     * @param poly float array where [lat, lon, lat, lon,...]
     * @param polyInDegrees true of poly floats represent decimal
     *        degrees.
     * @return true for Geo in poly
     */
    public static boolean isPointInPolygon(Geo x, float[] poly,
                                           boolean polyInDegrees) {
        Geo[] rg = new Geo[poly.length / 2];
        for (int j = 0; j < poly.length / 2; j++) {
            rg[j] = new Geo(poly[j * 2], poly[j * 2 + 1], polyInDegrees);
        }
        return isPointInPolygon(x, rg);
    }

    /**
     * return true IFF some point of the first argument is inside the
     * region specified by the closed polygon specified by the second
     * argument
     */
    public static boolean isPolylineInsidePolygon(float[] poly, float[] region) {
        int l = poly.length / 2;
        Geo[] rg = new Geo[region.length / 2];
        {
            for (int j = 0; j < region.length / 2; j++) {
                rg[j] = new Geo(region[j * 2], region[j * 2 + 1]);
            }
        }

        for (int i = 0; i < l; i++) {
            if (isPointInPolygon(new Geo(poly[i * 2], poly[i * 2 + 1]), rg)) {
                return true;
            }
        }
        return false;
    }

    /**
     * Returns the point of intersection of two great circle segments
     * defined by the segments. (lat1, lon1) to (lat2, lon2) and
     * (lat2, lon2) to (lat4, lon4). All lat-lon values are in
     * degrees.
     * 
     * @return a float array of length 4 containing upto 2 valid
     *         lat-lon points of intersection that lie on both
     *         segments. Positions in the array not containing a valid
     *         lat/lon value are initialized to Float.MAX_VALUE.
     */
    public static float[] getSegIntersection(float lat1, float lon1,
                                             float lat2, float lon2,
                                             float lat3, float lon3,
                                             float lat4, float lon4) {
        // KRA 03SEP03: The original version of this consed 26+ Geo's.
        // This one conses 8+.

        Geo p1 = new Geo(lat1, lon1);
        Geo p2 = new Geo(lat2, lon2);
        Geo p3 = new Geo(lat3, lon3);
        Geo p4 = new Geo(lat4, lon4);

        Geo geoCross1 = p1.crossNormalize(p2);
        Geo geoCross2 = p3.crossNormalize(p4);

        Geo i1 = geoCross1.crossNormalize(geoCross2);
        Geo i2 = i1.antipode();

        // check if the point of intersection lies on both segs
        // length of seg1
        // double d1 = Geo.distance(lat1, lon1, lat2, lon2);
        double d1 = p1.distance(p2);
        // length of seg2
        // double d2 = Geo.distance(lat3, lon3, lat4, lon4);
        double d2 = p3.distance(p4);

        // between seg1 endpoints and first point of intersection
        // double d111 = Geo.distance(lat1, lon1, ll[0], ll[1]);
        double d111 = p1.distance(i1);
        // double d121 = Geo.distance(lat2, lon2, ll[0], ll[1]);
        double d121 = p2.distance(i1);

        // between seg1 endpoints and second point of intersection
        // double d112 = Geo.distance(lat1, lon1, ll[2], ll[3]);
        double d112 = p1.distance(i2);
        // double d122 = Geo.distance(lat2, lon2, ll[2], ll[3]);
        double d122 = p2.distance(i2);

        // between seg2 endpoints and first point of intersection
        // double d211 = Geo.distance(lat3, lon3, ll[0], ll[1]);
        double d211 = p3.distance(i1);
        // double d221 = Geo.distance(lat4, lon4, ll[0], ll[1]);
        double d221 = p4.distance(i1);

        // between seg2 endpoints and second point of intersection
        // double d212 = Geo.distance(lat3, lon3, ll[2], ll[3]);
        double d212 = p3.distance(i2);
        // double d222 = Geo.distance(lat4, lon4, ll[2], ll[3]);
        double d222 = p4.distance(i2);

        float[] llp = new float[] { Float.MAX_VALUE, Float.MAX_VALUE,
                Float.MAX_VALUE, Float.MAX_VALUE };

        // check if first point of intersection lies on both segments
        if (d1 >= d111 && d1 >= d121 && d2 >= d211 && d2 >= d221) {
            llp[0] = ((float) i1.getLatitude());
            llp[1] = ((float) i1.getLongitude());
        }
        // check if second point of intersection lies on both segments
        if (d1 >= d112 && d1 >= d122 && d2 >= d212 && d2 >= d222) {
            llp[2] = ((float) i2.getLatitude());
            llp[3] = ((float) i2.getLongitude());
        }
        return llp;
    }

    // /**
    // * returns the point of interection of two great circle segments
    // * defined by the segments. (lat1, lon1) to (lat2, lon2) and
    // * (lat2, lon2) to (lat4, lon4). All lat-lon values are in
    // * degrees.
    // *
    // * @return a float array of length 4 containing upto 2 valid
    // * lat-lon points of intersection that lie on both
    // * segments. Positions in the array not containing a valid
    // * lat/lon value are initialized to Float.MAX_VALUE.
    // */
    // public static float[] getSegIntersectionOrig(float lat1, float
    // lon1,
    // float lat2, float lon2,
    // float lat3, float lon3,
    // float lat4, float lon4) {
    // // KRA 03SEP03: We can do better than this.
    //
    // float[] ll = getIntersection(lat1,
    // lon1,
    // lat2,
    // lon2,
    // lat3,
    // lon3,
    // lat4,
    // lon4);
    //
    // // check if the point of intersection lies on both segs
    //
    // // length of seg1
    // double d1 = Geo.distance(lat1, lon1, lat2, lon2);
    // // length of seg2
    // double d2 = Geo.distance(lat3, lon3, lat4, lon4);
    //
    // // between seg1 endpoints and first point of intersection
    // double d111 = Geo.distance(lat1, lon1, ll[0], ll[1]);
    // double d121 = Geo.distance(lat2, lon2, ll[0], ll[1]);
    //
    // // between seg1 endpoints and second point of intersection
    // double d112 = Geo.distance(lat1, lon1, ll[2], ll[3]);
    // double d122 = Geo.distance(lat2, lon2, ll[2], ll[3]);
    //
    // // between seg2 endpoints and first point of intersection
    // double d211 = Geo.distance(lat3, lon3, ll[0], ll[1]);
    // double d221 = Geo.distance(lat4, lon4, ll[0], ll[1]);
    //
    // // between seg2 endpoints and second point of intersection
    // double d212 = Geo.distance(lat3, lon3, ll[2], ll[3]);
    // double d222 = Geo.distance(lat4, lon4, ll[2], ll[3]);
    //
    // float[] llp = new float[] { Float.MAX_VALUE, Float.MAX_VALUE,
    // Float.MAX_VALUE, Float.MAX_VALUE };
    //
    // // check if first point of intersection lies on both segments
    // if (d1 >= d111 && d1 >= d121 && d2 >= d211 && d2 >= d221) {
    // llp[0] = ll[0];
    // llp[1] = ll[1];
    // }
    //
    // // check if second point of intersection lies on both segments
    // if (d1 >= d112 && d1 >= d122 && d2 >= d212 && d2 >= d222) {
    // llp[2] = ll[2];
    // llp[3] = ll[3];
    // }
    //
    // return llp;
    // }

    /**
     * Does the segment come within near radians of the region defined
     * by rCenter at rRadius?
     */
    public static final boolean isSegmentNearRadialRegion(GeoSegment segment,
                                                          Geo rCenter,
                                                          double rRadius,
                                                          double near) {
        Geo[] s = segment.getSeg();
        if (s != null && s.length == 2) {
            return isSegmentNearRadialRegion(s[0], s[1], rCenter, rRadius, near);
        }
        return false;
    }

    /**
     * Does the segment come within near radians of the region defined
     * by rCenter at rRadius?
     */
    public static final boolean isSegmentNearRadialRegion(Geo s1, Geo s2,
                                                          Geo rCenter,
                                                          double rRadius,
                                                          double near) {
        return s1.isInside(s2, near + rRadius, rCenter);
    }

    /** Is a segment horizontally within range of a Region region? */
    public static final boolean isSegmentNearRegion(GeoSegment segment,
                                                    double hrange,
                                                    GeoRegion region) {
        // Need to be careful here - calling
        // region.isSegmentNear(segment, hrange) can result in
        // circular code if the region just calls this method, which
        // may seem reasonable, if you look at the API.
        return isSegmentNearPolyRegion(segment, region.toPointArray(), hrange);
    }

    /**
     * Does the segment come within near radians of the region defined
     * by the polygon in r[*]? Catches segments within poly region and
     * returns after first hit, which is why it returns boolean.
     */
    public static final boolean isSegmentNearPolyRegion(GeoSegment segment,
                                                        Geo[] r, double near) {
        Geo[] s = segment.getSeg();
        if (s != null && s.length == 2) {
            return isSegmentNearPolyRegion(s[0], s[1], r, near);
        }
        return false;
    }

    /**
     * Does the segment s1-s2 come within near radians of the region
     * defined by the polygon in r[*]? Catches segments within poly