geo.java

来自「OpenMap是一个基于JavaBeansTM的开发工具包。利用OpenMap你」· Java 代码 · 共 1,182 行 · 第 1/3 页

    public Geo crossNormalize(Geo b) {
        return crossNormalize(b, new Geo());
    }

    /**
     * Eqvivalent to <code>this.cross(b).normalize()</code>.
     * 
     * @return ret Do not pass in a null value.
     */
    public Geo crossNormalize(Geo b, Geo ret) {
        double x = this.y() * b.z() - this.z() * b.y();
        double y = this.z() * b.x() - this.x() * b.z();
        double z = this.x() * b.y() - this.y() * b.x();
        double L = Math.sqrt(x * x + y * y + z * z);

        ret.initialize(x / L, y / L, z / L);
        return ret;
    }

    /**
     * Eqvivalent to <code>this.cross(b).normalize()</code>.
     * 
     * @return ret Do not pass in a null value.
     */
    public static Geo crossNormalize(Geo a, Geo b, Geo ret) {
        return a.crossNormalize(b, ret);
    }

    /** Returns this + b. */
    public Geo add(Geo b) {
        return add(b, new Geo());
    }

    /*
     * @return ret Do not pass in a null value.
     */
    public Geo add(Geo b, Geo ret) {
        ret.initialize(this.x() + b.x(), this.y() + b.y(), this.z() + b.z());
        return ret;
    }

    /** Returns this - b. */
    public Geo subtract(Geo b) {
        return subtract(b, new Geo());
    }

    /**
     * Returns this - b. *
     * 
     * @return ret Do not pass in a null value.
     */
    public Geo subtract(Geo b, Geo ret) {
        ret.initialize(this.x() - b.x(), this.y() - b.y(), this.z() - b.z());
        return ret;
    }

    public boolean equals(Geo v2) {
        return this.x == v2.x && this.y == v2.y && this.z == v2.z;
    }

    /** Angular distance, in radians between this and v2. */
    public double distance(Geo v2) {
        return Math.atan2(v2.crossLength(this), v2.dot(this));
    }

    /** Angular distance, in radians between v1 and v2. */
    public static double distance(Geo v1, Geo v2) {
        return v1.distance(v2);
    }

    /** Angular distance, in radians between the two lat lon points. */
    public static double distance(double lat1, double lon1, double lat2,
                                  double lon2) {
        return Geo.distance(new Geo(lat1, lon1), new Geo(lat2, lon2));
    }

    /** Distance in kilometers. * */
    public double distanceKM(Geo v2) {
        return km(distance(v2));
    }

    /** Distance in kilometers. * */
    public static double distanceKM(Geo v1, Geo v2) {
        return v1.distanceKM(v2);
    }

    /** Distance in kilometers. * */
    public static double distanceKM(double lat1, double lon1, double lat2,
                                    double lon2) {
        return Geo.distanceKM(new Geo(lat1, lon1), new Geo(lat2, lon2));
    }

    /** Distance in nautical miles. * */
    public double distanceNM(Geo v2) {
        return nm(distance(v2));
    }

    /** Distance in nautical miles. * */
    public static double distanceNM(Geo v1, Geo v2) {
        return v1.distanceNM(v2);
    }

    /** Distance in nautical miles. * */
    public static double distanceNM(double lat1, double lon1, double lat2,
                                    double lon2) {
        return Geo.distanceNM(new Geo(lat1, lon1), new Geo(lat2, lon2));
    }

    /** Azimuth in radians from this to v2. */
    public double azimuth(Geo v2) {
        /*
         * n1 is the great circle representing the meridian of this. n2 is the
         * great circle between this and v2. The azimuth is the angle between
         * them but we specialized the cross product.
         */
        // Geo n1 = north.cross(this);
        // Geo n2 = v2.cross(this);
        // crossNormalization is needed to geos of different length.
        Geo n1 = north.crossNormalize(this);
        Geo n2 = v2.crossNormalize(this);
        double az = Math.atan2(-north.dot(n2), n1.dot(n2));
        return (az > 0.0) ? az : 2.0 * Math.PI + az;
    }

    /**
     * Given 3 points on a sphere, p0, p1, p2, return the angle between them in
     * radians.
     */
    public static double angle(Geo p0, Geo p1, Geo p2) {
        return Math.PI - p0.cross(p1).distance(p1.cross(p2));
    }

    /**
     * Computes the area of a polygon on the surface of a unit sphere given an
     * enumeration of its point.. For a non unit sphere, multiply this by the
     * radius of sphere squared.
     */
    public static double area(Enumeration vs) {
        int count = 0;
        double area = 0;
        Geo v0 = (Geo) vs.nextElement();
        Geo v1 = (Geo) vs.nextElement();
        Geo p0 = v0;
        Geo p1 = v1;
        Geo p2 = null;
        while (vs.hasMoreElements()) {
            count = count + 1;
            p2 = (Geo) vs.nextElement();
            area = area + angle(p0, p1, p2);
            p0 = p1;
            p1 = p2;
        }

        count = count + 1;
        p2 = v0;
        area = area + angle(p0, p1, p2);
        p0 = p1;
        p1 = p2;

        count = count + 1;
        p2 = v1;
        area = area + angle(p0, p1, p2);

        return area - (count - 2) * Math.PI;
    }

    /**
     * Is the point, p, within radius radians of the great circle segment
     * between this and v2?
     */
    public boolean isInside(Geo v2, double radius, Geo p) {
        // Allocate a Geo to be reused for all of these calculations, instead of
        // creating 3 of them that are just thrown away. There's one more we
        // still need to allocate, for dp below.
        Geo tmp = new Geo();

        /*
         * gc is a unit vector perpendicular to the plane defined by v1 and v2
         */
        Geo gc = this.crossNormalize(v2, tmp);

        /*
         * |gc . p| is the size of the projection of p onto gc (the normal of
         * v1,v2) cos(pi/2-r) is effectively the size of the projection of a
         * vector along gc of the radius length. If the former is larger than
         * the latter, than p is further than radius from arc, so must not be
         * isInside
         */
        if (Math.abs(gc.dot(p)) > Math.cos((Math.PI / 2.0) - radius))
            return false;

        /*
         * If p is within radius of either endpoint, then we know it isInside
         */
        if (this.distance(p) <= radius || v2.distance(p) <= radius)
            return true;

        /* d is the vector from the v2 to v1 */
        Geo d = v2.subtract(this, tmp);

        /* L is the length of the vector d */
        double L = d.length();

        /* n is the d normalized to length=1 */
        Geo n = d.normalize(tmp);

        /* dp is the vector from p to v1 */
        Geo dp = p.subtract(this, new Geo());

        /* size is the size of the projection of dp onto n */
        double size = n.dot(dp);

        /* p is inside iff size>=0 and size <= L */
        return (0 <= size && size <= L);
    }

    /**
     * do the segments v1-v2 and p1-p2 come within radius (radians) of each
     * other?
     */
    public static boolean isInside(Geo v1, Geo v2, double radius, Geo p1, Geo p2) {
        return v1.isInside(v2, radius, p1) || v1.isInside(v2, radius, p2)
                || p1.isInside(p2, radius, v1) || p1.isInside(p2, radius, v2);
    }

    /**
     * Static version of isInside uses conventional (decimal degree)
     * coordinates.
     */
    public static boolean isInside(double lat1, double lon1, double lat2,
                                   double lon2, double radius, double lat3,
                                   double lon3) {
        return (new Geo(lat1, lon1)).isInside(new Geo(lat2, lon2),
                radius,
                new Geo(lat3, lon3));
    }

    /**
     * Is Geo p inside the time bubble along the great circle segment from this
     * to v2 looking forward forwardRadius and backward backwardRadius.
     */
    public boolean inBubble(Geo v2, double forwardRadius, double backRadius,
                            Geo p) {
        return distance(p) <= ((v2.subtract(this)
                .normalize()
                .dot(p.subtract(this)) > 0.0) ? forwardRadius : backRadius);
    }

    /** Returns the point opposite this point on the earth. */
    public Geo antipode() {
        return this.scale(-1.0, new Geo());
    }

    /**
     * Returns the point opposite this point on the earth. *
     * 
     * @return ret Do not pass in a null value.
     */
    public Geo antipode(Geo ret) {
        return this.scale(-1.0, ret);
    }

    /**
     * Find the intersection of the great circle between this and q and the
     * great circle normal to r.
     * <p>
     * 
     * That is, find the point, y, lying between this and q such that
     * 
     * <pre>
     *                              
     *  y = [x*this + (1-x)*q]*c
     *  where c = 1/y.dot(y) is a factor for normalizing y.
     *  y.dot(r) = 0
     *  substituting:
     *  [x*this + (1-x)*q]*c.dot(r) = 0 or
     *  [x*this + (1-x)*q].dot(r) = 0
     *  x*this.dot(r) + (1-x)*q.dot(r) = 0
     *  x*a + (1-x)*b = 0
     *  x = -b/(a - b)
     *                                        
     * </pre>
     * 
     * We assume that this and q are less than 180 degrees appart. When this and
     * q are 180 degrees appart, the point -y is also a valid intersection.
     * <p>
     * Alternatively the intersection point, y, satisfies y.dot(r) = 0
     * y.dot(this.crossNormalize(q)) = 0 which is satisfied by y =
     * r.crossNormalize(this.crossNormalize(q));
     * 
     */
    public Geo intersect(Geo q, Geo r) {
        return intersect(q, r, new Geo());
    }

    /**
     * Find the intersection of the great circle between this and q and the
     * great circle normal to r.
     * <p>
     * 
     * That is, find the point, y, lying between this and q such that
     * 
     * <pre>
     *                              
     *  y = [x*this + (1-x)*q]*c
     *  where c = 1/y.dot(y) is a factor for normalizing y.
     *  y.dot(r) = 0
     *  substituting:
     *  [x*this + (1-x)*q]*c.dot(r) = 0 or
     *  [x*this + (1-x)*q].dot(r) = 0
     *  x*this.dot(r) + (1-x)*q.dot(r) = 0
     *  x*a + (1-x)*b = 0
     *  x = -b/(a - b)
     *                                        
     * </pre>
     * 
     * We assume that this and q are less than 180 degrees appart. When this and
     * q are 180 degrees appart, the point -y is also a valid intersection.
     * <p>
     * Alternatively the intersection point, y, satisfies y.dot(r) = 0
     * y.dot(this.crossNormalize(q)) = 0 which is satisfied by y =
     * r.crossNormalize(this.crossNormalize(q));
     * 
     * @return ret Do not pass in a null value.
     */
    public Geo intersect(Geo q, Geo r, Geo ret) {
        double a = this.dot(r);
        double b = q.dot(r);
        double x = -b / (a - b);
        // This still results in one Geo being allocated and lost, in the
        // q.scale call.
        return this.scale(x, ret).add(q.scale(1.0 - x), ret).normalize(ret);
    }

    /** alias for computeCorridor(path, radius, radians(10), true) * */
    public static Geo[] computeCorridor(Geo[] path, double radius) {
        return computeCorridor(path, radius, radians(10.0), true);
    }

    /**
     * Wrap a fixed-distance corridor around an (open) path, as specified by an
     * array of Geo.
     * 
     * @param path Open path, must not have repeated points or consecutive
     *        antipodes.
     * @param radius Distance from path to widen corridor, in angular radians.
     * @param err maximum angle of rounded edges, in radians. If 0, will
     *        directly cut outside bends.
     * @param capp iff true, will round end caps
     * @return a closed polygon representing the specified corridor around the
     *         path.
     * 
     */
    public static Geo[] computeCorridor(Geo[] path, double radius, double err,
                                        boolean capp) {
        if (path == null || radius <= 0.0) {
            return new Geo[] {};
        }
        // assert path!=null;
        // assert radius > 0.0;

        int pl = path.length;
        if (pl < 2)
            return null;

        // final polygon will be right[0],...,right[n],left[m],...,left[0]
        ArrayList right = new ArrayList((int) (pl * 1.5));
        ArrayList left = new ArrayList((int) (pl * 1.5));

        Geo g0 = null; // previous point
        Geo n0 = null; // previous normal vector
        Geo l0 = null;
        Geo r0 = null;

        Geo g1 = path[0]; // current point

        for (int i = 1; i < pl; i++) {
            Geo g2 = path[i]; // next point
            Geo n1 = g1.crossNormalize(g2); // n is perpendicular to the vector
            // from g1 to g2
            n1 = n1.scale(radius); // normalize to radius
            // these are the offsets on the g2 side at g1
            Geo r1b = g1.add(n1);
            Geo l1b = g1.subtract(n1);

            if (n0 == null) {
                if (capp && err > 0) {
                    // start cap
                    Geo[] arc = approximateArc(g1, l1b, r1b, err);
                    for (int j = arc.length - 1; j >= 0; j--) {
                        right.add(arc[j]);
                    }
                } else {
                    // no previous point - we'll just be square