Class TrianglePredicate

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
13	package org.locationtech.jts.triangulate.quadedge;
14
15	import org.locationtech.jts.geom.Coordinate;
16	import org.locationtech.jts.geom.Triangle;
17	import org.locationtech.jts.geom.impl.CoordinateArraySequence;
18	import org.locationtech.jts.io.WKTWriter;
19	import org.locationtech.jts.math.DD;
20
21	/**
22	* Algorithms for computing values and predicates
23	* associated with triangles.
24	* For some algorithms extended-precision
25	* implementations are provided, which are more robust
26	* (i.e. they produce correct answers in more cases).
27	* Also, some more robust formulations of
28	* some algorithms are provided, which utilize
29	* normalization to the origin.
30	*
31	* @author Martin Davis
32	*
33	*/
34	public class TrianglePredicate
35	{
36	/**
37	* Tests if a point is inside the circle defined by
38	* the triangle with vertices a, b, c (oriented counter-clockwise).
39	* This test uses simple
40	* double-precision arithmetic, and thus may not be robust.
41	*
42	* @param a a vertex of the triangle
43	* @param b a vertex of the triangle
44	* @param c a vertex of the triangle
45	* @param p the point to test
46	* @return true if this point is inside the circle defined by the points a, b, c
47	*/
48	public static boolean isInCircleNonRobust(
49	Coordinate a, Coordinate b, Coordinate c,
50	Coordinate p) {
51	boolean isInCircle =
52	(a.x * a.x + a.y * a.y) * triArea(b, c, p)
53	- (b.x * b.x + b.y * b.y) * triArea(a, c, p)
54	+ (c.x * c.x + c.y * c.y) * triArea(a, b, p)
55	- (p.x * p.x + p.y * p.y) * triArea(a, b, c)
56	> 0;
57	return isInCircle;
58	}
59
60	/**
61	* Tests if a point is inside the circle defined by
62	* the triangle with vertices a, b, c (oriented counter-clockwise).
63	* This test uses simple
64	* double-precision arithmetic, and thus is not 100% robust.
65	* However, by using normalization to the origin
66	* it provides improved robustness and increased performance.
67	* <p>
68	* Based on code by J.R.Shewchuk.
69	*
70	*
71	* @param a a vertex of the triangle
72	* @param b a vertex of the triangle
73	* @param c a vertex of the triangle
74	* @param p the point to test
75	* @return true if this point is inside the circle defined by the points a, b, c
76	*/
77	public static boolean isInCircleNormalized(
78	Coordinate a, Coordinate b, Coordinate c,
79	Coordinate p) {
80	double adx = a.x - p.x;
81	double ady = a.y - p.y;
82	double bdx = b.x - p.x;
83	double bdy = b.y - p.y;
84	double cdx = c.x - p.x;
85	double cdy = c.y - p.y;
86
87	double abdet = adx * bdy - bdx * ady;
88	double bcdet = bdx * cdy - cdx * bdy;
89	double cadet = cdx * ady - adx * cdy;
90	double alift = adx * adx + ady * ady;
91	double blift = bdx * bdx + bdy * bdy;
92	double clift = cdx * cdx + cdy * cdy;
93
94	double disc = alift * bcdet + blift * cadet + clift * abdet;
95	return disc > 0;
96	}
97
98	/**
99	* Computes twice the area of the oriented triangle (a, b, c), i.e., the area is positive if the
100	* triangle is oriented counterclockwise.
101	*
102	* @param a a vertex of the triangle
103	* @param b a vertex of the triangle
104	* @param c a vertex of the triangle
105	*/
106	private static double triArea(Coordinate a, Coordinate b, Coordinate c) {
107	return (b.x - a.x) * (c.y - a.y)
108	- (b.y - a.y) * (c.x - a.x);
109	}
110
111	/**
112	* Tests if a point is inside the circle defined by
113	* the triangle with vertices a, b, c (oriented counter-clockwise).
114	* This method uses more robust computation.
115	*
116	* @param a a vertex of the triangle
117	* @param b a vertex of the triangle
118	* @param c a vertex of the triangle
119	* @param p the point to test
120	* @return true if this point is inside the circle defined by the points a, b, c
121	*/
122	public static boolean isInCircleRobust(
123	Coordinate a, Coordinate b, Coordinate c,
124	Coordinate p)
125	{
126	//checkRobustInCircle(a, b, c, p);
127	// return isInCircleNonRobust(a, b, c, p);
128	return isInCircleNormalized(a, b, c, p);
129	}
130
131	/**
132	* Tests if a point is inside the circle defined by
133	* the triangle with vertices a, b, c (oriented counter-clockwise).
134	* The computation uses {@link DD} arithmetic for robustness.
135	*
136	* @param a a vertex of the triangle
137	* @param b a vertex of the triangle
138	* @param c a vertex of the triangle
139	* @param p the point to test
140	* @return true if this point is inside the circle defined by the points a, b, c
141	*/
142	public static boolean isInCircleDDSlow(
143	Coordinate a, Coordinate b, Coordinate c,
144	Coordinate p) {
145	DD px = DD.valueOf(p.x);
146	DD py = DD.valueOf(p.y);
147	DD ax = DD.valueOf(a.x);
148	DD ay = DD.valueOf(a.y);
149	DD bx = DD.valueOf(b.x);
150	DD by = DD.valueOf(b.y);
151	DD cx = DD.valueOf(c.x);
152	DD cy = DD.valueOf(c.y);
153
154	DD aTerm = (ax.multiply(ax).add(ay.multiply(ay)))
155	.multiply(triAreaDDSlow(bx, by, cx, cy, px, py));
156	DD bTerm = (bx.multiply(bx).add(by.multiply(by)))
157	.multiply(triAreaDDSlow(ax, ay, cx, cy, px, py));
158	DD cTerm = (cx.multiply(cx).add(cy.multiply(cy)))
159	.multiply(triAreaDDSlow(ax, ay, bx, by, px, py));
160	DD pTerm = (px.multiply(px).add(py.multiply(py)))
161	.multiply(triAreaDDSlow(ax, ay, bx, by, cx, cy));
162
163	DD sum = aTerm.subtract(bTerm).add(cTerm).subtract(pTerm);
164	boolean isInCircle = sum.doubleValue() > 0;
165
166	return isInCircle;
167	}
168
169	/**
170	* Computes twice the area of the oriented triangle (a, b, c), i.e., the area
171	* is positive if the triangle is oriented counterclockwise.
172	* The computation uses {@link DD} arithmetic for robustness.
173	*
174	* @param ax the x ordinate of a vertex of the triangle
175	* @param ay the y ordinate of a vertex of the triangle
176	* @param bx the x ordinate of a vertex of the triangle
177	* @param by the y ordinate of a vertex of the triangle
178	* @param cx the x ordinate of a vertex of the triangle
179	* @param cy the y ordinate of a vertex of the triangle
180	*/
181	public static DD triAreaDDSlow(DD ax, DD ay,
182	DD bx, DD by, DD cx, DD cy) {
183	return (bx.subtract(ax).multiply(cy.subtract(ay)).subtract(by.subtract(ay)
184	.multiply(cx.subtract(ax))));
185	}
186
187	public static boolean isInCircleDDFast(
188	Coordinate a, Coordinate b, Coordinate c,
189	Coordinate p) {
190	DD aTerm = (DD.sqr(a.x).selfAdd(DD.sqr(a.y)))
191	.selfMultiply(triAreaDDFast(b, c, p));
192	DD bTerm = (DD.sqr(b.x).selfAdd(DD.sqr(b.y)))
193	.selfMultiply(triAreaDDFast(a, c, p));
194	DD cTerm = (DD.sqr(c.x).selfAdd(DD.sqr(c.y)))
195	.selfMultiply(triAreaDDFast(a, b, p));
196	DD pTerm = (DD.sqr(p.x).selfAdd(DD.sqr(p.y)))
197	.selfMultiply(triAreaDDFast(a, b, c));
198
199	DD sum = aTerm.selfSubtract(bTerm).selfAdd(cTerm).selfSubtract(pTerm);
200	boolean isInCircle = sum.doubleValue() > 0;
201
202	return isInCircle;
203	}
204
205	public static DD triAreaDDFast(
206	Coordinate a, Coordinate b, Coordinate c) {
207
208	DD t1 = DD.valueOf(b.x).selfSubtract(a.x)
209	.selfMultiply(
210	DD.valueOf(c.y).selfSubtract(a.y));
211
212	DD t2 = DD.valueOf(b.y).selfSubtract(a.y)
213	.selfMultiply(
214	DD.valueOf(c.x).selfSubtract(a.x));
215
216	return t1.selfSubtract(t2);
217	}
218
219	public static boolean isInCircleDDNormalized(
220	Coordinate a, Coordinate b, Coordinate c,
221	Coordinate p) {
222	DD adx = DD.valueOf(a.x).selfSubtract(p.x);
223	DD ady = DD.valueOf(a.y).selfSubtract(p.y);
224	DD bdx = DD.valueOf(b.x).selfSubtract(p.x);
225	DD bdy = DD.valueOf(b.y).selfSubtract(p.y);
226	DD cdx = DD.valueOf(c.x).selfSubtract(p.x);
227	DD cdy = DD.valueOf(c.y).selfSubtract(p.y);
228
229	DD abdet = adx.multiply(bdy).selfSubtract(bdx.multiply(ady));
230	DD bcdet = bdx.multiply(cdy).selfSubtract(cdx.multiply(bdy));
231	DD cadet = cdx.multiply(ady).selfSubtract(adx.multiply(cdy));
232	DD alift = adx.multiply(adx).selfAdd(ady.multiply(ady));
233	DD blift = bdx.multiply(bdx).selfAdd(bdy.multiply(bdy));
234	DD clift = cdx.multiply(cdx).selfAdd(cdy.multiply(cdy));
235
236	DD sum = alift.selfMultiply(bcdet)
237	.selfAdd(blift.selfMultiply(cadet))
238	.selfAdd(clift.selfMultiply(abdet));
239
240	boolean isInCircle = sum.doubleValue() > 0;
241
242	return isInCircle;
243	}
244
245	/**
246	* Computes the inCircle test using distance from the circumcentre.
247	* Uses standard double-precision arithmetic.
248	* <p>
249	* In general this doesn't
250	* appear to be any more robust than the standard calculation. However, there
251	* is at least one case where the test point is far enough from the
252	* circumcircle that this test gives the correct answer.
253	* <pre>
254	* LINESTRING
255	* (1507029.9878 518325.7547, 1507022.1120341457 518332.8225183258,
256	* 1507029.9833 518325.7458, 1507029.9896965567 518325.744909031)
257	* </pre>
258	*
259	* @param a a vertex of the triangle
260	* @param b a vertex of the triangle
261	* @param c a vertex of the triangle
262	* @param p the point to test
263	* @return true if this point is inside the circle defined by the points a, b, c
264	*/
265	public static boolean isInCircleCC(Coordinate a, Coordinate b, Coordinate c,
266	Coordinate p) {
267	Coordinate cc = Triangle.circumcentre(a, b, c);
268	double ccRadius = a.distance(cc);
269	double pRadiusDiff = p.distance(cc) - ccRadius;
270	return pRadiusDiff <= 0;
271	}
272
273	/**
274	* Checks if the computed value for isInCircle is correct, using
275	* double-double precision arithmetic.
276	*
277	* @param a a vertex of the triangle
278	* @param b a vertex of the triangle
279	* @param c a vertex of the triangle
280	* @param p the point to test
281	*/
282	private static void checkRobustInCircle(Coordinate a, Coordinate b, Coordinate c,
283	Coordinate p)
284	{
285	boolean nonRobustInCircle = isInCircleNonRobust(a, b, c, p);
286	boolean isInCircleDD = TrianglePredicate.isInCircleDDSlow(a, b, c, p);
287	boolean isInCircleCC = TrianglePredicate.isInCircleCC(a, b, c, p);
288
289	Coordinate circumCentre = Triangle.circumcentre(a, b, c);
290	System.out.println("p radius diff a = "
291	+ Math.abs(p.distance(circumCentre) - a.distance(circumCentre))
292	/ a.distance(circumCentre));
293
294	if (nonRobustInCircle != isInCircleDD \|\| nonRobustInCircle != isInCircleCC) {
295	System.out.println("inCircle robustness failure (double result = "
296	+ nonRobustInCircle
297	+ ", DD result = " + isInCircleDD
298	+ ", CC result = " + isInCircleCC + ")");
299	System.out.println(WKTWriter.toLineString(new CoordinateArraySequence(
300	new Coordinate[] { a, b, c, p })));
301	System.out.println("Circumcentre = " + WKTWriter.toPoint(circumCentre)
302	+ " radius = " + a.distance(circumCentre));
303	System.out.println("p radius diff a = "
304	+ Math.abs(p.distance(circumCentre)/a.distance(circumCentre) - 1));
305	System.out.println("p radius diff b = "
306	+ Math.abs(p.distance(circumCentre)/b.distance(circumCentre) - 1));
307	System.out.println("p radius diff c = "
308	+ Math.abs(p.distance(circumCentre)/c.distance(circumCentre) - 1));
309	System.out.println();
310	}
311	}
312
313
314	}
315