Doc Classic

Class Intersection

  1 /*
  2  * Copyright (c) 2019 martin Davis
  3  *
  4  * All rights reserved. This program and the accompanying materials
  5  * are made available under the terms of the Eclipse Public License 2.0
  6  * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
  7  * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
  8  * and the Eclipse Distribution License is available at
  9  *
 10  * http://www.eclipse.org/org/documents/edl-v10.php.
 11  */
 12 package org.locationtech.jts.algorithm;
 13  
 14 import org.locationtech.jts.geom.Coordinate;
 15  
 16 /**
 17  * Contains functions to compute intersections between lines.
 18  * 
 19  * @author Martin Davis
 20  *
 21  */
 22 public class Intersection {
 23   
 24   /**
 25    * Computes the intersection point of two lines.
 26    * If the lines are parallel or collinear this case is detected 
 27    * and <code>null</code> is returned.
 28    * <p>
 29    * In general it is not possible to accurately compute
 30    * the intersection point of two lines, due to 
 31    * numerical roundoff.
 32    * This is particularly true when the input lines are nearly parallel.
 33    * This routine uses numerical conditioning on the input values
 34    * to ensure that the computed value should be very close to the correct value.
 35    * 
 36    * @param p1 an endpoint of line 1
 37    * @param p2 an endpoint of line 1
 38    * @param q1 an endpoint of line 2
 39    * @param q2 an endpoint of line 2
 40    * @return the intersection point between the lines, if there is one,
 41    * or null if the lines are parallel or collinear
 42    * 
 43    * @see CGAlgorithmsDD#intersection(Coordinate, Coordinate, Coordinate, Coordinate)
 44    */
 45   public static Coordinate intersection(Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2) {
 46     // compute midpoint of "kernel envelope"
 47     double minX0 = p1.x < p2.x ? p1.x : p2.x;
 48     double minY0 = p1.y < p2.y ? p1.y : p2.y;
 49     double maxX0 = p1.x > p2.x ? p1.x : p2.x;
 50     double maxY0 = p1.y > p2.y ? p1.y : p2.y;
 51  
 52     double minX1 = q1.x < q2.x ? q1.x : q2.x;
 53     double minY1 = q1.y < q2.y ? q1.y : q2.y;
 54     double maxX1 = q1.x > q2.x ? q1.x : q2.x;
 55     double maxY1 = q1.y > q2.y ? q1.y : q2.y;
 56  
 57     double intMinX = minX0 > minX1 ? minX0 : minX1;
 58     double intMaxX = maxX0 < maxX1 ? maxX0 : maxX1;
 59     double intMinY = minY0 > minY1 ? minY0 : minY1;
 60     double intMaxY = maxY0 < maxY1 ? maxY0 : maxY1;
 61  
 62     double midx = (intMinX + intMaxX) / 2.0;
 63     double midy = (intMinY + intMaxY) / 2.0;
 64     
 65     // condition ordinate values by subtracting midpoint
 66     double p1x = p1.x - midx;
 67     double p1y = p1.y - midy;
 68     double p2x = p2.x - midx;
 69     double p2y = p2.y - midy;
 70     double q1x = q1.x - midx;
 71     double q1y = q1.y - midy;
 72     double q2x = q2.x - midx;
 73     double q2y = q2.y - midy;
 74      
 75     // unrolled computation using homogeneous coordinates eqn
 76     double px = p1y - p2y;
 77     double py = p2x - p1x;
 78     double pw = p1x * p2y - p2x * p1y;
 79     
 80     double qx = q1y - q2y;
 81     double qy = q2x - q1x;
 82     double qw = q1x * q2y - q2x * q1y;
 83     
 84     double x = py * qw - qy * pw;
 85     double y = qx * pw - px * qw;
 86     double w = px * qy - qx * py;
 87     
 88     double xInt = x/w;
 89     double yInt = y/w;
 90     
 91     // check for parallel lines
 92     if ((Double.isNaN(xInt)) || (Double.isInfinite(xInt)
 93         || Double.isNaN(yInt)) || (Double.isInfinite(yInt))) {
 94       return null;
 95     }
 96     // de-condition intersection point
 97     return new Coordinate(xInt + midx, yInt + midy);
 98   }
 99  
100 }
101  
102  
103