CGAlgorithmsDD.java
/*
* Copyright (c) 2016 Martin Davis.
*
* 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.math.DD;
/**
* Implements basic computational geometry algorithms using {@link DD} arithmetic.
*
* @author Martin Davis
*
*/
public class CGAlgorithmsDD
{
private CGAlgorithmsDD() {}
/**
* 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
*/
public static int orientationIndex(Coordinate p1, Coordinate p2, Coordinate q)
{
return orientationIndex(p1.x, p1.y, p2.x, p2.y, q.x, q.y);
}
/**
* Returns the index of the direction of the point {@code q} relative to
* a vector specified by {@code p1-p2}.
*
* @param p1x the x ordinate of the vector origin point
* @param p1y the y ordinate of the vector origin point
* @param p2x the x ordinate of the vector final point
* @param p2y the y ordinate of the vector final point
* @param qx the x ordinate of the query point
* @param qy the y ordinate of the query point
*
* @return 1 if q is counter-clockwise (left) from p1-p2
* -1 if q is clockwise (right) from p1-p2
* 0 if q is collinear with p1-p2
*/
public static int orientationIndex(double p1x, double p1y,
double p2x, double p2y,
double qx, double qy)
{
// fast filter for orientation index
// avoids use of slow extended-precision arithmetic in many cases
int index = orientationIndexFilter(p1x, p1y, p2x, p2y, qx, qy);
if (index <= 1) return index;
// normalize coordinates
DD dx1 = DD.valueOf(p2x).selfAdd(-p1x);
DD dy1 = DD.valueOf(p2y).selfAdd(-p1y);
DD dx2 = DD.valueOf(qx).selfAdd(-p2x);
DD dy2 = DD.valueOf(qy).selfAdd(-p2y);
// sign of determinant - unrolled for performance
return dx1.selfMultiply(dy2).selfSubtract(dy1.selfMultiply(dx2)).signum();
}
/**
* Computes the sign of the determinant of the 2x2 matrix
* with the given entries.
*
* @return -1 if the determinant is negative,
* 1 if the determinant is positive,
* 0 if the determinant is 0.
*/
public static int signOfDet2x2(DD x1, DD y1, DD x2, DD y2)
{
DD det = x1.multiply(y2).selfSubtract(y1.multiply(x2));
return det.signum();
}
/**
* Computes the sign of the determinant of the 2x2 matrix
* with the given entries.
*
* @return -1 if the determinant is negative,
* 1 if the determinant is positive,
* 0 if the determinant is 0.
*/
public static int signOfDet2x2(double dx1, double dy1, double dx2, double dy2)
{
DD x1 = DD.valueOf(dx1);
DD y1 = DD.valueOf(dy1);
DD x2 = DD.valueOf(dx2);
DD y2 = DD.valueOf(dy2);
DD det = x1.multiply(y2).selfSubtract(y1.multiply(x2));
return det.signum();
}
/**
* A value which is safely greater than the
* relative round-off error in double-precision numbers
*/
private static final double DP_SAFE_EPSILON = 1e-15;
/**
* A filter for computing the orientation index of three coordinates.
* <p>
* If the orientation can be computed safely using standard DP
* arithmetic, this routine returns the orientation index.
* Otherwise, a value i > 1 is returned.
* In this case the orientation index must
* be computed using some other more robust method.
* The filter is fast to compute, so can be used to
* avoid the use of slower robust methods except when they are really needed,
* thus providing better average performance.
* <p>
* Uses an approach due to Jonathan Shewchuk, which is in the public domain.
*
* @param pax A coordinate
* @param pay A coordinate
* @param pbx B coordinate
* @param pby B coordinate
* @param pcx C coordinate
* @param pcy C coordinate
* @return the orientation index if it can be computed safely
* @return i > 1 if the orientation index cannot be computed safely
*/
private static int orientationIndexFilter(double pax, double pay,
double pbx, double pby, double pcx, double pcy)
{
double detsum;
double detleft = (pax - pcx) * (pby - pcy);
double detright = (pay - pcy) * (pbx - pcx);
double det = detleft - detright;
if (detleft > 0.0) {
if (detright <= 0.0) {
return signum(det);
}
else {
detsum = detleft + detright;
}
}
else if (detleft < 0.0) {
if (detright >= 0.0) {
return signum(det);
}
else {
detsum = -detleft - detright;
}
}
else {
return signum(det);
}
double errbound = DP_SAFE_EPSILON * detsum;
if ((det >= errbound) || (-det >= errbound)) {
return signum(det);
}
return 2;
}
private static int signum(double x)
{
if (x > 0) return 1;
if (x < 0) return -1;
return 0;
}
/**
* Computes an intersection point between two lines
* using DD arithmetic.
* If the lines are parallel (either identical
* or separate) a null value is returned.
*
* @param p1 an endpoint of line segment 1
* @param p2 an endpoint of line segment 1
* @param q1 an endpoint of line segment 2
* @param q2 an endpoint of line segment 2
* @return an intersection point if one exists, or null if the lines are parallel
*/
public static Coordinate intersection(
Coordinate p1, Coordinate p2,
Coordinate q1, Coordinate q2)
{
DD px = new DD(p1.y).selfSubtract(p2.y);
DD py = new DD(p2.x).selfSubtract(p1.x);
DD pw = new DD(p1.x).selfMultiply(p2.y).selfSubtract(new DD(p2.x).selfMultiply(p1.y));
DD qx = new DD(q1.y).selfSubtract(q2.y);
DD qy = new DD(q2.x).selfSubtract(q1.x);
DD qw = new DD(q1.x).selfMultiply(q2.y).selfSubtract(new DD(q2.x).selfMultiply(q1.y));
DD x = py.multiply(qw).selfSubtract(qy.multiply(pw));
DD y = qx.multiply(pw).selfSubtract(px.multiply(qw));
DD w = px.multiply(qy).selfSubtract(qx.multiply(py));
double xInt = x.selfDivide(w).doubleValue();
double yInt = y.selfDivide(w).doubleValue();
if ((Double.isNaN(xInt)) || (Double.isInfinite(xInt) || Double.isNaN(yInt)) || (Double.isInfinite(yInt))) {
return null;
}
return new Coordinate(xInt, yInt);
}
}