Class CGAlgorithmsDD

1	/*
2	* Copyright (c) 2016 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	import org.locationtech.jts.math.DD;
16
17	/**
18	* Implements basic computational geometry algorithms using {@link DD} arithmetic.
19	*
20	* @author Martin Davis
21	*
22	*/
23	public class CGAlgorithmsDD
24	{
25	private CGAlgorithmsDD() {}
26
27	/**
28	* Returns the index of the direction of the point <code>q</code> relative to
29	* a vector specified by <code>p1-p2</code>.
30	*
31	* @param p1 the origin point of the vector
32	* @param p2 the final point of the vector
33	* @param q the point to compute the direction to
34	*
35	* @return 1 if q is counter-clockwise (left) from p1-p2
36	* @return -1 if q is clockwise (right) from p1-p2
37	* @return 0 if q is collinear with p1-p2
38	*/
39	public static int orientationIndex(Coordinate p1, Coordinate p2, Coordinate q)
40	{
41	// fast filter for orientation index
42	// avoids use of slow extended-precision arithmetic in many cases
43	int index = orientationIndexFilter(p1, p2, q);
44	if (index <= 1) return index;
45
46	// normalize coordinates
47	DD dx1 = DD.valueOf(p2.x).selfAdd(-p1.x);
48	DD dy1 = DD.valueOf(p2.y).selfAdd(-p1.y);
49	DD dx2 = DD.valueOf(q.x).selfAdd(-p2.x);
50	DD dy2 = DD.valueOf(q.y).selfAdd(-p2.y);
51
52	// sign of determinant - unrolled for performance
53	return dx1.selfMultiply(dy2).selfSubtract(dy1.selfMultiply(dx2)).signum();
54	}
55
56	/**
57	* Computes the sign of the determinant of the 2x2 matrix
58	* with the given entries.
59	*
60	* @return -1 if the determinant is negative,
61	* @return 1 if the determinant is positive,
62	* @return 0 if the determinant is 0.
63	*/
64	public static int signOfDet2x2(DD x1, DD y1, DD x2, DD y2)
65	{
66	DD det = x1.multiply(y2).selfSubtract(y1.multiply(x2));
67	return det.signum();
68	}
69
70	/**
71	* Computes the sign of the determinant of the 2x2 matrix
72	* with the given entries.
73	*
74	* @return -1 if the determinant is negative,
75	* @return 1 if the determinant is positive,
76	* @return 0 if the determinant is 0.
77	*/
78	public static int signOfDet2x2(double dx1, double dy1, double dx2, double dy2)
79	{
80	DD x1 = DD.valueOf(dx1);
81	DD y1 = DD.valueOf(dy1);
82	DD x2 = DD.valueOf(dx2);
83	DD y2 = DD.valueOf(dy2);
84
85	DD det = x1.multiply(y2).selfSubtract(y1.multiply(x2));
86	return det.signum();
87	}
88
89	/**
90	* A value which is safely greater than the
91	* relative round-off error in double-precision numbers
92	*/
93	private static final double DP_SAFE_EPSILON = 1e-15;
94
95	/**
96	* A filter for computing the orientation index of three coordinates.
97	* <p>
98	* If the orientation can be computed safely using standard DP
99	* arithmetic, this routine returns the orientation index.
100	* Otherwise, a value i > 1 is returned.
101	* In this case the orientation index must
102	* be computed using some other more robust method.
103	* The filter is fast to compute, so can be used to
104	* avoid the use of slower robust methods except when they are really needed,
105	* thus providing better average performance.
106	* <p>
107	* Uses an approach due to Jonathan Shewchuk, which is in the public domain.
108	*
109	* @param pa a coordinate
110	* @param pb a coordinate
111	* @param pc a coordinate
112	* @return the orientation index if it can be computed safely
113	* @return i > 1 if the orientation index cannot be computed safely
114	*/
115	private static int orientationIndexFilter(Coordinate pa, Coordinate pb, Coordinate pc)
116	{
117	double detsum;
118
119	double detleft = (pa.x - pc.x) * (pb.y - pc.y);
120	double detright = (pa.y - pc.y) * (pb.x - pc.x);
121	double det = detleft - detright;
122
123	if (detleft > 0.0) {
124	if (detright <= 0.0) {
125	return signum(det);
126	}
127	else {
128	detsum = detleft + detright;
129	}
130	}
131	else if (detleft < 0.0) {
132	if (detright >= 0.0) {
133	return signum(det);
134	}
135	else {
136	detsum = -detleft - detright;
137	}
138	}
139	else {
140	return signum(det);
141	}
142
143	double errbound = DP_SAFE_EPSILON * detsum;
144	if ((det >= errbound) \|\| (-det >= errbound)) {
145	return signum(det);
146	}
147
148	return 2;
149	}
150
151	private static int signum(double x)
152	{
153	if (x > 0) return 1;
154	if (x < 0) return -1;
155	return 0;
156	}
157
158	/**
159	* Computes an intersection point between two lines
160	* using DD arithmetic.
161	* If the lines are parallel (either identical
162	* or separate) a null value is returned.
163	*
164	* @param p1 an endpoint of line segment 1
165	* @param p2 an endpoint of line segment 1
166	* @param q1 an endpoint of line segment 2
167	* @param q2 an endpoint of line segment 2
168	* @return an intersection point if one exists, or null if the lines are parallel
169	*/
170	public static Coordinate intersection(
171	Coordinate p1, Coordinate p2,
172	Coordinate q1, Coordinate q2)
173	{
174	DD px = new DD(p1.y).selfSubtract(p2.y);
175	DD py = new DD(p2.x).selfSubtract(p1.x);
176	DD pw = new DD(p1.x).selfMultiply(p2.y).selfSubtract(new DD(p2.x).selfMultiply(p1.y));
177
178	DD qx = new DD(q1.y).selfSubtract(q2.y);
179	DD qy = new DD(q2.x).selfSubtract(q1.x);
180	DD qw = new DD(q1.x).selfMultiply(q2.y).selfSubtract(new DD(q2.x).selfMultiply(q1.y));
181
182	DD x = py.multiply(qw).selfSubtract(qy.multiply(pw));
183	DD y = qx.multiply(pw).selfSubtract(px.multiply(qw));
184	DD w = px.multiply(qy).selfSubtract(qx.multiply(py));
185
186	double xInt = x.selfDivide(w).doubleValue();
187	double yInt = y.selfDivide(w).doubleValue();
188
189	if ((Double.isNaN(xInt)) \|\| (Double.isInfinite(xInt) \|\| Double.isNaN(yInt)) \|\| (Double.isInfinite(yInt))) {
190	return null;
191	}
192
193	return new Coordinate(xInt, yInt);
194	}
195	}
196