/* * 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; import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.math.DD; /** * Implements basic computational geometry algorithms using {@link DD} arithmetic. * * @author Martin Davis * */ public class CGAlgorithmsDD { /** * Returns the index of the direction of the point <code>q</code> relative to * a vector specified by <code>p1-p2</code>. * * @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 1 if q is counter-clockwise (left) from p1-p2 * @return -1 if q is clockwise (right) from p1-p2 * @return 0 if q is collinear with p1-p2 */ public static int orientationIndex(Coordinate p1, Coordinate p2, Coordinate q) { // fast filter for orientation index // avoids use of slow extended-precision arithmetic in many cases int index = orientationIndexFilter(p1, p2, q); if (index <= 1) return index; // normalize coordinates DD dx1 = DD.valueOf(p2.x).selfAdd(-p1.x); DD dy1 = DD.valueOf(p2.y).selfAdd(-p1.y); DD dx2 = DD.valueOf(q.x).selfAdd(-p2.x); DD dy2 = DD.valueOf(q.y).selfAdd(-p2.y); // 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, * @return 1 if the determinant is positive, * @return 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(); } /** * 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 pa a coordinate * @param pb a coordinate * @param pc a 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(Coordinate pa, Coordinate pb, Coordinate pc) { double detsum; double detleft = (pa.x - pc.x) * (pb.y - pc.y); double detright = (pa.y - pc.y) * (pb.x - pc.x); 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. * Currently does not handle case of parallel lines. * * @param p1 * @param p2 * @param q1 * @param q2 * @return */ public static Coordinate intersection( Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2) { DD denom1 = DD.valueOf(q2.y).selfSubtract(q1.y) .selfMultiply(DD.valueOf(p2.x).selfSubtract(p1.x)); DD denom2 = DD.valueOf(q2.x).selfSubtract(q1.x) .selfMultiply(DD.valueOf(p2.y).selfSubtract(p1.y)); DD denom = denom1.subtract(denom2); /** * Cases: * - denom is 0 if lines are parallel * - intersection point lies within line segment p if fracP is between 0 and 1 * - intersection point lies within line segment q if fracQ is between 0 and 1 */ DD numx1 = DD.valueOf(q2.x).selfSubtract(q1.x) .selfMultiply(DD.valueOf(p1.y).selfSubtract(q1.y)); DD numx2 = DD.valueOf(q2.y).selfSubtract(q1.y) .selfMultiply(DD.valueOf(p1.x).selfSubtract(q1.x)); DD numx = numx1.subtract(numx2); double fracP = numx.selfDivide(denom).doubleValue(); double x = DD.valueOf(p1.x).selfAdd(DD.valueOf(p2.x).selfSubtract(p1.x).selfMultiply(fracP)).doubleValue(); DD numy1 = DD.valueOf(p2.x).selfSubtract(p1.x) .selfMultiply(DD.valueOf(p1.y).selfSubtract(q1.y)); DD numy2 = DD.valueOf(p2.y).selfSubtract(p1.y) .selfMultiply(DD.valueOf(p1.x).selfSubtract(q1.x)); DD numy = numy1.subtract(numy2); double fracQ = numy.selfDivide(denom).doubleValue(); double y = DD.valueOf(p1.y).selfAdd(DD.valueOf(p2.y).selfSubtract(p1.y).selfMultiply(fracQ)).doubleValue(); return new Coordinate(x,y); } }