

/*
 * The JTS Topology Suite is a collection of Java classes that
 * implement the fundamental operations required to validate a given
 * geo-spatial data set to a known topological specification.
 *
 * Copyright (C) 2001 Vivid Solutions
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * For more information, contact:
 *
 *     Vivid Solutions
 *     Suite #1A
 *     2328 Government Street
 *     Victoria BC  V8T 5G5
 *     Canada
 *
 *     (250)385-6040
 *     www.vividsolutions.com
 */
package com.vividsolutions.jts.algorithm;

/**
 *@version 1.7
 */

import com.vividsolutions.jts.geom.*;
import com.vividsolutions.jts.util.Assert;

/**
 * A robust version of {@link LineIntersector}.
 *
 * @version 1.7
 * @see RobustDeterminant
 */
public class RobustLineIntersector
    extends LineIntersector
{

  public RobustLineIntersector() {
  }

  public void computeIntersection(Coordinate p, Coordinate p1, Coordinate p2) {
    isProper = false;
    // do between check first, since it is faster than the orientation test
    if (Envelope.intersects(p1, p2, p)) {
      if ((CGAlgorithms.orientationIndex(p1, p2, p) == 0)
          && (CGAlgorithms.orientationIndex(p2, p1, p) == 0)) {
        isProper = true;
        if (p.equals(p1) || p.equals(p2)) {
          isProper = false;
        }
        result = POINT_INTERSECTION;
        return;
      }
    }
    result = NO_INTERSECTION;
  }

  protected int computeIntersect(
                Coordinate p1, Coordinate p2,
                Coordinate q1, Coordinate q2  ) {
    isProper = false;

    // first try a fast test to see if the envelopes of the lines intersect
    if (! Envelope.intersects(p1, p2, q1, q2))
      return NO_INTERSECTION;

    // for each endpoint, compute which side of the other segment it lies
    // if both endpoints lie on the same side of the other segment,
    // the segments do not intersect
    int Pq1 = CGAlgorithms.orientationIndex(p1, p2, q1);
    int Pq2 = CGAlgorithms.orientationIndex(p1, p2, q2);

    if ((Pq1>0 && Pq2>0) || (Pq1<0 && Pq2<0)) {
      return NO_INTERSECTION;
    }

    int Qp1 = CGAlgorithms.orientationIndex(q1, q2, p1);
    int Qp2 = CGAlgorithms.orientationIndex(q1, q2, p2);

    if ((Qp1>0 && Qp2>0) || (Qp1<0 && Qp2<0)) {
        return NO_INTERSECTION;
    }

    boolean collinear = Pq1 == 0
         && Pq2 == 0
         && Qp1 == 0
         && Qp2 == 0;
    if (collinear) {
      return computeCollinearIntersection(p1, p2, q1, q2);
    }
    
    /**
     * At this point we know that there is a single intersection point
     * (since the lines are not collinear).
     */
    
    /**
     *  Check if the intersection is an endpoint. If it is, copy the endpoint as
     *  the intersection point. Copying the point rather than computing it
     *  ensures the point has the exact value, which is important for
     *  robustness. It is sufficient to simply check for an endpoint which is on
     *  the other line, since at this point we know that the inputLines must
     *  intersect.
     */
    if (Pq1 == 0 || Pq2 == 0 || Qp1 == 0 || Qp2 == 0) {
      isProper = false;
      
      /**
       * Check for two equal endpoints.  
       * This is done explicitly rather than by the orientation tests
       * below in order to improve robustness.
       * 
       * [An example where the orientation tests fail to be consistent is
       * the following (where the true intersection is at the shared endpoint
       * POINT (19.850257749638203 46.29709338043669)
       * 
       * LINESTRING ( 19.850257749638203 46.29709338043669, 20.31970698357233 46.76654261437082 ) 
       * and 
       * LINESTRING ( -48.51001596420236 -22.063180333403878, 19.850257749638203 46.29709338043669 )
       * 
       * which used to produce the INCORRECT result: (20.31970698357233, 46.76654261437082, NaN)
       * 
       */
      if (p1.equals2D(q1) 
      		|| p1.equals2D(q2)) {
      	intPt[0] = p1;
      }
      else if (p2.equals2D(q1) 
      		|| p2.equals2D(q2)) {
      	intPt[0] = p2;
      }
      
      /**
       * Now check to see if any endpoint lies on the interior of the other segment.
       */
      else if (Pq1 == 0) {
        intPt[0] = new Coordinate(q1);
      }
      else if (Pq2 == 0) {
        intPt[0] = new Coordinate(q2);
      }
      else if (Qp1 == 0) {
        intPt[0] = new Coordinate(p1);
      }
      else if (Qp2 == 0) {
        intPt[0] = new Coordinate(p2);
      }
    }
    else {
      isProper = true;
      intPt[0] = intersection(p1, p2, q1, q2);
    }
    return POINT_INTERSECTION;
  }

  private int computeCollinearIntersection(Coordinate p1, Coordinate p2,
      Coordinate q1, Coordinate q2) {
    boolean p1q1p2 = Envelope.intersects(p1, p2, q1);
    boolean p1q2p2 = Envelope.intersects(p1, p2, q2);
    boolean q1p1q2 = Envelope.intersects(q1, q2, p1);
    boolean q1p2q2 = Envelope.intersects(q1, q2, p2);

    if (p1q1p2 && p1q2p2) {
      intPt[0] = q1;
      intPt[1] = q2;
      return COLLINEAR_INTERSECTION;
    }
    if (q1p1q2 && q1p2q2) {
      intPt[0] = p1;
      intPt[1] = p2;
      return COLLINEAR_INTERSECTION;
    }
    if (p1q1p2 && q1p1q2) {
      intPt[0] = q1;
      intPt[1] = p1;
      return q1.equals(p1) && !p1q2p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
    }
    if (p1q1p2 && q1p2q2) {
      intPt[0] = q1;
      intPt[1] = p2;
      return q1.equals(p2) && !p1q2p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
    }
    if (p1q2p2 && q1p1q2) {
      intPt[0] = q2;
      intPt[1] = p1;
      return q2.equals(p1) && !p1q1p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
    }
    if (p1q2p2 && q1p2q2) {
      intPt[0] = q2;
      intPt[1] = p2;
      return q2.equals(p2) && !p1q1p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
    }
    return NO_INTERSECTION;
  }

  /**
   * This method computes the actual value of the intersection point.
   * To obtain the maximum precision from the intersection calculation,
   * the coordinates are normalized by subtracting the minimum
   * ordinate values (in absolute value).  This has the effect of
   * removing common significant digits from the calculation to
   * maintain more bits of precision.
   */
  private Coordinate intersection(
    Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2)
  {
    Coordinate intPt = intersectionWithNormalization(p1, p2, q1, q2);
    
    /*
    // TESTING ONLY
    Coordinate intPtDD = CGAlgorithmsDD.intersection(p1, p2, q1, q2);
    double dist = intPt.distance(intPtDD);
    System.out.println(intPt + " - " + intPtDD + " dist = " + dist);
    //intPt = intPtDD;
    */
    
    /**
     * Due to rounding it can happen that the computed intersection is
     * outside the envelopes of the input segments.  Clearly this
     * is inconsistent. 
     * This code checks this condition and forces a more reasonable answer
     * 
     * MD - May 4 2005 - This is still a problem.  Here is a failure case:
     *
     * LINESTRING (2089426.5233462777 1180182.3877339689, 2085646.6891757075 1195618.7333999649)
     * LINESTRING (1889281.8148903656 1997547.0560044837, 2259977.3672235999 483675.17050843034)
     * int point = (2097408.2633752143,1144595.8008114607)
     * 
     * MD - Dec 14 2006 - This does not seem to be a failure case any longer
     */
    if (! isInSegmentEnvelopes(intPt)) {
//      System.out.println("Intersection outside segment envelopes: " + intPt);
//      System.out.println("Segments: " + this);
      // compute a safer result
      intPt = CentralEndpointIntersector.getIntersection(p1, p2, q1, q2);
//      System.out.println("Snapped to " + intPt);
    }

    if (precisionModel != null) {
      precisionModel.makePrecise(intPt);
    }

    return intPt;
  }

  private Coordinate intersectionWithNormalization(
    Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2)
    {
      Coordinate n1 = new Coordinate(p1);
      Coordinate n2 = new Coordinate(p2);
      Coordinate n3 = new Coordinate(q1);
      Coordinate n4 = new Coordinate(q2);
      Coordinate normPt = new Coordinate();
      normalizeToEnvCentre(n1, n2, n3, n4, normPt);

      Coordinate intPt = safeHCoordinateIntersection(n1, n2, n3, n4);

      intPt.x += normPt.x;
      intPt.y += normPt.y;
      
      return intPt;
  }
  
  /**
   * Computes a segment intersection using homogeneous coordinates.
   * Round-off error can cause the raw computation to fail, 
   * (usually due to the segments being approximately parallel).
   * If this happens, a reasonable approximation is computed instead.
   * 
   * @param p1 a segment endpoint
   * @param p2 a segment endpoint
   * @param q1 a segment endpoint
   * @param q2 a segment endpoint
   * @return the computed intersection point
   */
  private Coordinate safeHCoordinateIntersection(Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2)
  {
    Coordinate intPt = null;
    try {
      intPt = HCoordinate.intersection(p1, p2, q1, q2);
    }
    catch (NotRepresentableException e) {
//    	System.out.println("Not calculable: " + this);
      // compute an approximate result
      intPt = CentralEndpointIntersector.getIntersection(p1, p2, q1, q2);
 //     System.out.println("Snapped to " + intPt);
    }
    return intPt;
  }

  /**
   * Normalize the supplied coordinates so that
   * their minimum ordinate values lie at the origin.
   * NOTE: this normalization technique appears to cause
   * large errors in the position of the intersection point for some cases.
   *
   * @param n1
   * @param n2
   * @param n3
   * @param n4
   * @param normPt
   */
  private void normalizeToMinimum(
    Coordinate n1,
    Coordinate n2,
    Coordinate n3,
    Coordinate n4,
    Coordinate normPt)
  {
    normPt.x = smallestInAbsValue(n1.x, n2.x, n3.x, n4.x);
    normPt.y = smallestInAbsValue(n1.y, n2.y, n3.y, n4.y);
    n1.x -= normPt.x;    n1.y -= normPt.y;
    n2.x -= normPt.x;    n2.y -= normPt.y;
    n3.x -= normPt.x;    n3.y -= normPt.y;
    n4.x -= normPt.x;    n4.y -= normPt.y;
  }

  /**
   * Normalize the supplied coordinates to
   * so that the midpoint of their intersection envelope
   * lies at the origin.
   *
   * @param n00
   * @param n01
   * @param n10
   * @param n11
   * @param normPt
   */
  private void normalizeToEnvCentre(
    Coordinate n00,
    Coordinate n01,
    Coordinate n10,
    Coordinate n11,
    Coordinate normPt)
  {
    double minX0 = n00.x < n01.x ? n00.x : n01.x;
    double minY0 = n00.y < n01.y ? n00.y : n01.y;
    double maxX0 = n00.x > n01.x ? n00.x : n01.x;
    double maxY0 = n00.y > n01.y ? n00.y : n01.y;

    double minX1 = n10.x < n11.x ? n10.x : n11.x;
    double minY1 = n10.y < n11.y ? n10.y : n11.y;
    double maxX1 = n10.x > n11.x ? n10.x : n11.x;
    double maxY1 = n10.y > n11.y ? n10.y : n11.y;

    double intMinX = minX0 > minX1 ? minX0 : minX1;
    double intMaxX = maxX0 < maxX1 ? maxX0 : maxX1;
    double intMinY = minY0 > minY1 ? minY0 : minY1;
    double intMaxY = maxY0 < maxY1 ? maxY0 : maxY1;

    double intMidX = (intMinX + intMaxX) / 2.0;
    double intMidY = (intMinY + intMaxY) / 2.0;
    normPt.x = intMidX;
    normPt.y = intMidY;

    /*
    // equilavalent code using more modular but slower method
    Envelope env0 = new Envelope(n00, n01);
    Envelope env1 = new Envelope(n10, n11);
    Envelope intEnv = env0.intersection(env1);
    Coordinate intMidPt = intEnv.centre();

    normPt.x = intMidPt.x;
    normPt.y = intMidPt.y;
    */

    n00.x -= normPt.x;    n00.y -= normPt.y;
    n01.x -= normPt.x;    n01.y -= normPt.y;
    n10.x -= normPt.x;    n10.y -= normPt.y;
    n11.x -= normPt.x;    n11.y -= normPt.y;
  }

  private double smallestInAbsValue(double x1, double x2, double x3, double x4)
  {
    double x = x1;
    double xabs = Math.abs(x);
    if (Math.abs(x2) < xabs) {
      x = x2;
      xabs = Math.abs(x2);
    }
    if (Math.abs(x3) < xabs) {
      x = x3;
      xabs = Math.abs(x3);
    }
    if (Math.abs(x4) < xabs) {
      x = x4;
    }
    return x;
  }

  /**
   * Test whether a point lies in the envelopes of both input segments.
   * A correctly computed intersection point should return <code>true</code>
   * for this test.
   * Since this test is for debugging purposes only, no attempt is
   * made to optimize the envelope test.
   *
   * @return <code>true</code> if the input point lies within both input segment envelopes
   */
  private boolean isInSegmentEnvelopes(Coordinate intPt)
  {
    Envelope env0 = new Envelope(inputLines[0][0], inputLines[0][1]);
    Envelope env1 = new Envelope(inputLines[1][0], inputLines[1][1]);
    return env0.contains(intPt) && env1.contains(intPt);
  }

}
