CGAlgorithms.java

/*
 * Copyright (c) 2016 Vivid Solutions.
 *
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License 2.0
 * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
 * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
 * and the Eclipse Distribution License is available at
 *
 * http://www.eclipse.org/org/documents/edl-v10.php.
 */
package org.locationtech.jts.algorithm;

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.CoordinateSequence;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Location;
import org.locationtech.jts.math.MathUtil;

/**
 * Specifies and implements various fundamental Computational Geometric
 * algorithms. The algorithms supplied in this class are robust for
 * double-precision floating point.
 * 
 * @version 1.7
 * @deprecated See {@link Length}, {@link Area}, {@link Distance},
 *             {@link Orientation}, {@link PointLocation}
 */
public class CGAlgorithms
{

  /**
   * A value that indicates an orientation of clockwise, or a right turn.
   * 
   * @deprecated Use {@link Orientation#CLOCKWISE} instead.
   */
  public static final int CLOCKWISE = -1;

  /**
   * A value that indicates an orientation of clockwise, or a right turn.
   * 
   * @deprecated Use {@link Orientation#RIGHT} instead.
   */
  public static final int RIGHT = CLOCKWISE;

  /**
   * A value that indicates an orientation of counterclockwise, or a left turn.
   * 
   * @deprecated Use {@link Orientation#COUNTERCLOCKWISE} instead.
   */
  public static final int COUNTERCLOCKWISE = 1;

  /**
   * A value that indicates an orientation of counterclockwise, or a left turn.
   * 
   * @deprecated Use {@link Orientation#LEFT} instead.
   */
  public static final int LEFT = COUNTERCLOCKWISE;

  /**
   * A value that indicates an orientation of collinear, or no turn (straight).
   * 
   * @deprecated Use {@link Orientation#COLLINEAR} instead.
   */
  public static final int COLLINEAR = 0;

  /**
   * A value that indicates an orientation of collinear, or no turn (straight).
   * 
   * @deprecated Use {@link Orientation#STRAIGHT} instead.
   */
  public static final int STRAIGHT = COLLINEAR;

  /**
   * Returns the index of the direction of the point {@code q} relative to
   * a vector specified by {@code p1-p2}.
   * 
   * @param p1 the origin point of the vector
   * @param p2 the final point of the vector
   * @param q the point to compute the direction to
   * 
   * @return {@code 1} if q is counter-clockwise (left) from p1-p2
   *        {@code -1} if q is clockwise (right) from p1-p2
   *         {@code 0} if q is collinear with p1-p2
   *
   * @deprecated Use {@link Orientation#index(Coordinate, Coordinate, Coordinate)}
   *             instead.
   */
  public static int orientationIndex(Coordinate p1, Coordinate p2, Coordinate q)
  {
    /*
      MD - 9 Aug 2010 It seems that the basic algorithm is slightly orientation
      dependent, when computing the orientation of a point very close to a
      line. This is possibly due to the arithmetic in the translation to the
      origin.

      For instance, the following situation produces identical results in spite
      of the inverse orientation of the line segment:

      Coordinate p0 = new Coordinate(219.3649559090992, 140.84159161824724);
      Coordinate p1 = new Coordinate(168.9018919682399, -5.713787599646864);

      Coordinate p = new Coordinate(186.80814046338352, 46.28973405831556); int
      orient = orientationIndex(p0, p1, p); int orientInv =
      orientationIndex(p1, p0, p);

      A way to force consistent results is to normalize the orientation of the
      vector using the following code. However, this may make the results of
      orientationIndex inconsistent through the triangle of points, so it's not
      clear this is an appropriate patch.

     */
    return CGAlgorithmsDD.orientationIndex(p1, p2, q);
    // testing only
    //return ShewchuksDeterminant.orientationIndex(p1, p2, q);
    // previous implementation - not quite fully robust
    //return RobustDeterminant.orientationIndex(p1, p2, q);
    
  }

  public CGAlgorithms()
  {
  }

  /**
   * Tests whether a point lies inside or on a ring. The ring may be oriented in
   * either direction. A point lying exactly on the ring boundary is considered
   * to be inside the ring.
   * <p>
   * This method does <i>not</i> first check the point against the envelope of
   * the ring.
   * 
   * @param p
   *          point to check for ring inclusion
   * @param ring
   *          an array of coordinates representing the ring (which must have
   *          first point identical to last point)
   * @return true if p is inside ring
   * 
   * @see #locatePointInRing(Coordinate, Coordinate[])  
   * @deprecated Use {@link PointLocation#isInRing(Coordinate, Coordinate[])}
   *             instead.
   */
  public static boolean isPointInRing(Coordinate p, Coordinate[] ring)
  {
    return locatePointInRing(p, ring) != Location.EXTERIOR;
  }

  /**
   * Determines whether a point lies in the interior, on the boundary, or in the
   * exterior of a ring. The ring may be oriented in either direction.
   * <p>
   * This method does <i>not</i> first check the point against the envelope of
   * the ring.
   * 
   * @param p
   *          point to check for ring inclusion
   * @param ring
   *          an array of coordinates representing the ring (which must have
   *          first point identical to last point)
   * @return the {@link Location} of p relative to the ring
   * @deprecated Use
   *             {@link PointLocation#locateInRing(Coordinate, Coordinate[])}
   *             instead.
   */
  public static int locatePointInRing(Coordinate p, Coordinate[] ring)
  {
    return RayCrossingCounter.locatePointInRing(p, ring);
  }

  /**
   * Tests whether a point lies on the line segments defined by a list of
   * coordinates.
   * 
   * @return true if the point is a vertex of the line or lies in the interior
   *         of a line segment in the linestring
   * @deprecated Use {@link PointLocation#isOnLine(Coordinate, Coordinate[])}
   *             instead.
   */
  public static boolean isOnLine(Coordinate p, Coordinate[] pt)
  {
    for (int i = 1; i < pt.length; i++) {
      Coordinate p0 = pt[i - 1];
      Coordinate p1 = pt[i];
      if (PointLocation.isOnSegment(p, p0, p1)) {
        return true;
      }
    }
    return false;
  }

  /**
   * Computes whether a ring defined by an array of {@link Coordinate}s is
   * oriented counter-clockwise.
   * <ul>
   * <li>The list of points is assumed to have the first and last points equal.
   * <li>This will handle coordinate lists which contain repeated points.
   * </ul>
   * This algorithm is <b>only</b> guaranteed to work with valid rings. If the
   * ring is invalid (e.g. self-crosses or touches), the computed result may not
   * be correct.
   * 
   * @param ring
   *          an array of Coordinates forming a ring
   * @return true if the ring is oriented counter-clockwise.
   * @throws IllegalArgumentException
   *           if there are too few points to determine orientation (&lt; 4)
   * @deprecated Use {@link Orientation#isCCW(Coordinate[])} instead.
   */
  public static boolean isCCW(Coordinate[] ring)
  {
    // # of points without closing endpoint
    int nPts = ring.length - 1;
    // sanity check
    if (nPts < 3)
      throw new IllegalArgumentException(
          "Ring has fewer than 4 points, so orientation cannot be determined");

    // find highest point
    Coordinate hiPt = ring[0];
    int hiIndex = 0;
    for (int i = 1; i <= nPts; i++) {
      Coordinate p = ring[i];
      if (p.y > hiPt.y) {
        hiPt = p;
        hiIndex = i;
      }
    }

    // find distinct point before highest point
    int iPrev = hiIndex;
    do {
      iPrev = iPrev - 1;
      if (iPrev < 0)
        iPrev = nPts;
    } while (ring[iPrev].equals2D(hiPt) && iPrev != hiIndex);

    // find distinct point after highest point
    int iNext = hiIndex;
    do {
      iNext = (iNext + 1) % nPts;
    } while (ring[iNext].equals2D(hiPt) && iNext != hiIndex);

    Coordinate prev = ring[iPrev];
    Coordinate next = ring[iNext];

    /*
      This check catches cases where the ring contains an A-B-A configuration
      of points. This can happen if the ring does not contain 3 distinct points
      (including the case where the input array has fewer than 4 elements), or
      it contains coincident line segments.
     */
    if (prev.equals2D(hiPt) || next.equals2D(hiPt) || prev.equals2D(next))
      return false;

    int disc = computeOrientation(prev, hiPt, next);

    /*
      If disc is exactly 0, lines are collinear. There are two possible cases:
      (1) the lines lie along the x axis in opposite directions (2) the lines
      lie on top of one another

      (1) is handled by checking if next is left of prev ==> CCW (2) will never
      happen if the ring is valid, so don't check for it (Might want to assert
      this)
     */
    boolean isCCW;
    if (disc == 0) {
      // poly is CCW if prev x is right of next x
      isCCW = (prev.x > next.x);
    }
    else {
      // if area is positive, points are ordered CCW
      isCCW = (disc > 0);
    }
    return isCCW;
  }

  /**
   * Computes the orientation of a point q to the directed line segment p1-p2.
   * The orientation of a point relative to a directed line segment indicates
   * which way you turn to get to q after travelling from p1 to p2.
   * 
   * @param p1 the first vertex of the line segment
   * @param p2 the second vertex of the line segment
   * @param q the point to compute the relative orientation of
   * @return 1 if q is counter-clockwise from p1-p2,
   * or -1 if q is clockwise from p1-p2,
   * or 0 if q is collinear with p1-p2
   * @deprecated Use {@link Orientation#index(Coordinate, Coordinate, Coordinate)}
   *             instead.
   */
  public static int computeOrientation(Coordinate p1, Coordinate p2,
      Coordinate q)
  {
    return orientationIndex(p1, p2, q);
  }

  /**
   * Computes the distance from a point p to a line segment AB
   * 
   * Note: NON-ROBUST!
   * 
   * @param p
   *          the point to compute the distance for
   * @param A
   *          one point of the line
   * @param B
   *          another point of the line (must be different to A)
   * @return the distance from p to line segment AB
   * @deprecated Use
   *             {@link Distance#pointToSegment(Coordinate, Coordinate, Coordinate)}
   *             instead.
   */
  public static double distancePointLine(Coordinate p, Coordinate A,
      Coordinate B)
  {
    // if start = end, then just compute distance to one of the endpoints
    if (A.x == B.x && A.y == B.y)
      return p.distance(A);

    // otherwise use comp.graphics.algorithms Frequently Asked Questions method
    /*
     * (1) r = AC dot AB 
     *         --------- 
     *         ||AB||^2 
     *         
     * r has the following meaning: 
     *   r=0 P = A 
     *   r=1 P = B 
     *   r<0 P is on the backward extension of AB 
     *   r>1 P is on the forward extension of AB 
     *   0<r<1 P is interior to AB
     */

    double len2 = (B.x - A.x) * (B.x - A.x) + (B.y - A.y) * (B.y - A.y);
    double r = ((p.x - A.x) * (B.x - A.x) + (p.y - A.y) * (B.y - A.y))
        / len2;

    if (r <= 0.0)
      return p.distance(A);
    if (r >= 1.0)
      return p.distance(B);

    /*
     * (2) s = (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) 
     *         ----------------------------- 
     *                    L^2
     * 
     * Then the distance from C to P = |s|*L.
     * 
     * This is the same calculation as {@link #distancePointLinePerpendicular}.
     * Unrolled here for performance.
     */
    double s = ((A.y - p.y) * (B.x - A.x) - (A.x - p.x) * (B.y - A.y))
        / len2;
    return Math.abs(s) * Math.sqrt(len2);
  }

  /**
   * Computes the perpendicular distance from a point p to the (infinite) line
   * containing the points AB
   * 
   * @param p
   *          the point to compute the distance for
   * @param A
   *          one point of the line
   * @param B
   *          another point of the line (must be different to A)
   * @return the distance from p to line AB
   * @deprecated Use
   *             {@link Distance#pointToLinePerpendicular(Coordinate, Coordinate, Coordinate)}
   *             instead.
   */
  public static double distancePointLinePerpendicular(Coordinate p,
      Coordinate A, Coordinate B)
  {
    // use comp.graphics.algorithms Frequently Asked Questions method
    /*
     * (2) s = (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) 
     *         ----------------------------- 
     *                    L^2
     * 
     * Then the distance from C to P = |s|*L.
     */
    double len2 = (B.x - A.x) * (B.x - A.x) + (B.y - A.y) * (B.y - A.y);
    double s = ((A.y - p.y) * (B.x - A.x) - (A.x - p.x) * (B.y - A.y))
        / len2;

    return Math.abs(s) * Math.sqrt(len2);
  }

  /**
   * Computes the distance from a point to a sequence of line segments.
   * 
   * @param p
   *          a point
   * @param line
   *          a sequence of contiguous line segments defined by their vertices
   * @return the minimum distance between the point and the line segments
   * @deprecated Use
   *             {@link Distance#pointToSegmentString(Coordinate, Coordinate[])}
   *             instead.
   */
  public static double distancePointLine(Coordinate p, Coordinate[] line)
  {
    if (line.length == 0)
      throw new IllegalArgumentException(
          "Line array must contain at least one vertex");
    // this handles the case of length = 1
    double minDistance = p.distance(line[0]);
    for (int i = 0; i < line.length - 1; i++) {
      double dist = distancePointLine(p, line[i], line[i + 1]);
      if (dist < minDistance) {
        minDistance = dist;
      }
    }
    return minDistance;
  }

  /**
   * Computes the distance from a line segment AB to a line segment CD
   * 
   * Note: NON-ROBUST!
   * 
   * @param A
   *          a point of one line
   * @param B
   *          the second point of (must be different to A)
   * @param C
   *          one point of the line
   * @param D
   *          another point of the line (must be different to A)
   * @deprecated Use
   *             {@link Distance#segmentToSegment(Coordinate, Coordinate, Coordinate, Coordinate)}
   *             instead.
   */
  public static double distanceLineLine(Coordinate A, Coordinate B,
      Coordinate C, Coordinate D)
  {
    // check for zero-length segments
    if (A.equals(B))
      return distancePointLine(A, C, D);
    if (C.equals(D))
      return distancePointLine(D, A, B);

    // AB and CD are line segments
    /*
     * from comp.graphics.algo
     * 
     * Solving the above for r and s yields 
     * 
     *     (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy) 
     * r = ----------------------------- (eqn 1) 
     *     (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
     * 
     *     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)  
     * s = ----------------------------- (eqn 2)
     *     (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx) 
     *     
     * Let P be the position vector of the
     * intersection point, then 
     *   P=A+r(B-A) or 
     *   Px=Ax+r(Bx-Ax) 
     *   Py=Ay+r(By-Ay) 
     * By examining the values of r & s, you can also determine some other limiting
     * conditions: 
     *   If 0<=r<=1 & 0<=s<=1, intersection exists 
     *      r<0 or r>1 or s<0 or s>1 line segments do not intersect 
     *   If the denominator in eqn 1 is zero, AB & CD are parallel 
     *   If the numerator in eqn 1 is also zero, AB & CD are collinear.
     */

    boolean noIntersection = false;
    if (! Envelope.intersects(A, B, C, D)) {
      noIntersection = true;
    }
    else {
      double denom = (B.x - A.x) * (D.y - C.y) - (B.y - A.y) * (D.x - C.x);
      
      if (denom == 0) {
        noIntersection = true;
      }
      else {
        double r_num = (A.y - C.y) * (D.x - C.x) - (A.x - C.x) * (D.y - C.y);
        double s_num = (A.y - C.y) * (B.x - A.x) - (A.x - C.x) * (B.y - A.y);
        
        double s = s_num / denom;
        double r = r_num / denom;
  
        if ((r < 0) || (r > 1) || (s < 0) || (s > 1)) {
          noIntersection = true;
        }
      }
    }
    if (noIntersection) {
      return MathUtil.min(
            distancePointLine(A, C, D),
            distancePointLine(B, C, D),
            distancePointLine(C, A, B),
            distancePointLine(D, A, B));
    }
    // segments intersect
    return 0.0; 
  }

  /**
   * Computes the signed area for a ring. The signed area is positive if the
   * ring is oriented CW, negative if the ring is oriented CCW, and zero if the
   * ring is degenerate or flat.
   * 
   * @param ring
   *          the coordinates forming the ring
   * @return the signed area of the ring
   * @deprecated Use {@link Area#ofRing(Coordinate[])} or
   *             {@link Area#ofRingSigned(Coordinate[])} instead.
   */
  public static double signedArea(Coordinate[] ring)
  {
    if (ring.length < 3)
      return 0.0;
    double sum = 0.0;
    /*
      Based on the Shoelace formula.
      http://en.wikipedia.org/wiki/Shoelace_formula
     */
    double x0 = ring[0].x;
    for (int i = 1; i < ring.length - 1; i++) {
      double x = ring[i].x - x0;
      double y1 = ring[i + 1].y;
      double y2 = ring[i - 1].y;
      sum += x * (y2 - y1);
    }
    return sum / 2.0;
  }

  /**
   * Computes the signed area for a ring. The signed area is:
   * <ul>
   * <li>positive if the ring is oriented CW
   * <li>negative if the ring is oriented CCW
   * <li>zero if the ring is degenerate or flat
   * </ul>
   * 
   * @param ring
   *          the coordinates forming the ring
   * @return the signed area of the ring
   * @deprecated Use {@link Area#ofRing(CoordinateSequence)} or
   *             {@link Area#ofRingSigned(CoordinateSequence)} instead.
   */
  public static double signedArea(CoordinateSequence ring)
  {
    int n = ring.size();
    if (n < 3)
      return 0.0;
    /*
      Based on the Shoelace formula.
      http://en.wikipedia.org/wiki/Shoelace_formula
     */
    Coordinate p0 = new Coordinate();
    Coordinate p1 = new Coordinate();
    Coordinate p2 = new Coordinate();
    ring.getCoordinate(0, p1);
    ring.getCoordinate(1, p2);
    double x0 = p1.x;
    p2.x -= x0;
    double sum = 0.0;
    for (int i = 1; i < n - 1; i++) {
      p0.y = p1.y;
      p1.x = p2.x;
      p1.y = p2.y;
      ring.getCoordinate(i + 1, p2);
      p2.x -= x0;
      sum += p1.x * (p0.y - p2.y);
    }
    return sum / 2.0;
  }

  /**
   * Computes the length of a linestring specified by a sequence of points.
   * 
   * @param pts
   *          the points specifying the linestring
   * @return the length of the linestring
   * @deprecated Use {@link Length#ofLine(CoordinateSequence)} instead.
   */
  public static double length(CoordinateSequence pts)
  {
    // optimized for processing CoordinateSequences
    int n = pts.size();
    if (n <= 1)
      return 0.0;

    double len = 0.0;

    Coordinate p = new Coordinate();
    pts.getCoordinate(0, p);
    double x0 = p.x;
    double y0 = p.y;

    for (int i = 1; i < n; i++) {
      pts.getCoordinate(i, p);
      double x1 = p.x;
      double y1 = p.y;
      double dx = x1 - x0;
      double dy = y1 - y0;

      len += Math.hypot(dx, dy);

      x0 = x1;
      y0 = y1;
    }
    return len;
  }

}