Class RobustDeterminant

1
2
3	/*
4	* Copyright (c) 2016 Vivid Solutions.
5	*
6	* All rights reserved. This program and the accompanying materials
7	* are made available under the terms of the Eclipse Public License 2.0
8	* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
9	* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
10	* and the Eclipse Distribution License is available at
11	*
12	* http://www.eclipse.org/org/documents/edl-v10.php.
13	*/
14	package org.locationtech.jts.algorithm;
15
16	import org.locationtech.jts.geom.Coordinate;
17
18	/**
19	* @version 1.7
20	*/
21
22	/**
23	* Implements an algorithm to compute the
24	* sign of a 2x2 determinant for double precision values robustly.
25	* It is a direct translation of code developed by Olivier Devillers.
26	* <p>
27	* The original code carries the following copyright notice:
28	*
29	* <pre>
30	*************************************************************************
31	* Author : Olivier Devillers
32	* Olivier.Devillers@sophia.inria.fr
33	* http:/www.inria.fr:/prisme/personnel/devillers/anglais/determinant.html
34	*
35	* Relicensed under EDL and EPL with Permission from Olivier Devillers
36	*
37	**************************************************************************
38	*
39	**************************************************************************
40	* Copyright (c) 1995 by INRIA Prisme Project
41	* BP 93 06902 Sophia Antipolis Cedex, France.
42	* All rights reserved
43	**************************************************************************
44	* </pre>
45	*
46	* @version 1.7
47	*/
48	public class RobustDeterminant {
49
50	//public static int callCount = 0; // debugging only
51
52	/*
53	// test point to allow injecting test code
54	public static int signOfDet2x2(double x1, double y1, double x2, double y2)
55	{
56	int d1 = originalSignOfDet2x2(x1, y1, x2, y2);
57	int d2 = -originalSignOfDet2x2(y1, x1, x2, y2);
58	assert d1 == -d2;
59	return d1;
60	}
61	*/
62
63	/*
64	* Test code to force a standard ordering of input ordinates.
65	* A possible fix for a rare problem where evaluation is order-dependent.
66	*/
67	/*
68	public static int signOfDet2x2(double x1, double y1, double x2, double y2)
69	{
70	if (x1 > x2) {
71	return -signOfDet2x2ordX(x2, y2, x1, y1);
72	}
73	return signOfDet2x2ordX(x1, y1, x2, y2);
74	}
75
76	private static int signOfDet2x2ordX(double x1, double y1, double x2, double y2)
77	{
78	if (y1 > y2) {
79	return -originalSignOfDet2x2(y1, x1, y2, x2);
80	}
81	return originalSignOfDet2x2(x1, y1, x2, y2);
82	}
83	// */
84
85	/**
86	* Computes the sign of the determinant of the 2x2 matrix
87	* with the given entries, in a robust way.
88	*
89	* @return -1 if the determinant is negative,
90	* @return 1 if the determinant is positive,
91	* @return 0 if the determinant is 0.
92	*/
93	//private static int originalSignOfDet2x2(double x1, double y1, double x2, double y2) {
94	public static int signOfDet2x2(double x1, double y1, double x2, double y2) {
95	// returns -1 if the determinant is negative,
96	// returns 1 if the determinant is positive,
97	// returns 0 if the determinant is null.
98	int sign;
99	double swap;
100	double k;
101	long count = 0;
102
103	//callCount++; // debugging only
104
105	sign = 1;
106
107	/*
108	* testing null entries
109	*/
110	if ((x1 == 0.0) \|\| (y2 == 0.0)) {
111	if ((y1 == 0.0) \|\| (x2 == 0.0)) {
112	return 0;
113	}
114	else if (y1 > 0) {
115	if (x2 > 0) {
116	return -sign;
117	}
118	else {
119	return sign;
120	}
121	}
122	else {
123	if (x2 > 0) {
124	return sign;
125	}
126	else {
127	return -sign;
128	}
129	}
130	}
131	if ((y1 == 0.0) \|\| (x2 == 0.0)) {
132	if (y2 > 0) {
133	if (x1 > 0) {
134	return sign;
135	}
136	else {
137	return -sign;
138	}
139	}
140	else {
141	if (x1 > 0) {
142	return -sign;
143	}
144	else {
145	return sign;
146	}
147	}
148	}
149
150	/*
151	* making y coordinates positive and permuting the entries
152	*/
153	/*
154	* so that y2 is the biggest one
155	*/
156	if (0.0 < y1) {
157	if (0.0 < y2) {
158	if (y1 <= y2) {
159	;
160	}
161	else {
162	sign = -sign;
163	swap = x1;
164	x1 = x2;
165	x2 = swap;
166	swap = y1;
167	y1 = y2;
168	y2 = swap;
169	}
170	}
171	else {
172	if (y1 <= -y2) {
173	sign = -sign;
174	x2 = -x2;
175	y2 = -y2;
176	}
177	else {
178	swap = x1;
179	x1 = -x2;
180	x2 = swap;
181	swap = y1;
182	y1 = -y2;
183	y2 = swap;
184	}
185	}
186	}
187	else {
188	if (0.0 < y2) {
189	if (-y1 <= y2) {
190	sign = -sign;
191	x1 = -x1;
192	y1 = -y1;
193	}
194	else {
195	swap = -x1;
196	x1 = x2;
197	x2 = swap;
198	swap = -y1;
199	y1 = y2;
200	y2 = swap;
201	}
202	}
203	else {
204	if (y1 >= y2) {
205	x1 = -x1;
206	y1 = -y1;
207	x2 = -x2;
208	y2 = -y2;
209	;
210	}
211	else {
212	sign = -sign;
213	swap = -x1;
214	x1 = -x2;
215	x2 = swap;
216	swap = -y1;
217	y1 = -y2;
218	y2 = swap;
219	}
220	}
221	}
222
223	/*
224	* making x coordinates positive
225	*/
226	/*
227	* if \|x2\| < \|x1\| one can conclude
228	*/
229	if (0.0 < x1) {
230	if (0.0 < x2) {
231	if (x1 <= x2) {
232	;
233	}
234	else {
235	return sign;
236	}
237	}
238	else {
239	return sign;
240	}
241	}
242	else {
243	if (0.0 < x2) {
244	return -sign;
245	}
246	else {
247	if (x1 >= x2) {
248	sign = -sign;
249	x1 = -x1;
250	x2 = -x2;
251	;
252	}
253	else {
254	return -sign;
255	}
256	}
257	}
258
259	/*
260	* all entries strictly positive x1 <= x2 and y1 <= y2
261	*/
262	while (true) {
263	count = count + 1;
264	// MD - UNSAFE HACK for testing only!
265	// k = (int) (x2 / x1);
266	k = Math.floor(x2 / x1);
267	x2 = x2 - k * x1;
268	y2 = y2 - k * y1;
269
270	/*
271	* testing if R (new U2) is in U1 rectangle
272	*/
273	if (y2 < 0.0) {
274	return -sign;
275	}
276	if (y2 > y1) {
277	return sign;
278	}
279
280	/*
281	* finding R'
282	*/
283	if (x1 > x2 + x2) {
284	if (y1 < y2 + y2) {
285	return sign;
286	}
287	}
288	else {
289	if (y1 > y2 + y2) {
290	return -sign;
291	}
292	else {
293	x2 = x1 - x2;
294	y2 = y1 - y2;
295	sign = -sign;
296	}
297	}
298	if (y2 == 0.0) {
299	if (x2 == 0.0) {
300	return 0;
301	}
302	else {
303	return -sign;
304	}
305	}
306	if (x2 == 0.0) {
307	return sign;
308	}
309
310	/*
311	* exchange 1 and 2 role.
312	*/
313	// MD - UNSAFE HACK for testing only!
314	// k = (int) (x1 / x2);
315	k = Math.floor(x1 / x2);
316	x1 = x1 - k * x2;
317	y1 = y1 - k * y2;
318
319	/*
320	* testing if R (new U1) is in U2 rectangle
321	*/
322	if (y1 < 0.0) {
323	return sign;
324	}
325	if (y1 > y2) {
326	return -sign;
327	}
328
329	/*
330	* finding R'
331	*/
332	if (x2 > x1 + x1) {
333	if (y2 < y1 + y1) {
334	return -sign;
335	}
336	}
337	else {
338	if (y2 > y1 + y1) {
339	return sign;
340	}
341	else {
342	x1 = x2 - x1;
343	y1 = y2 - y1;
344	sign = -sign;
345	}
346	}
347	if (y1 == 0.0) {
348	if (x1 == 0.0) {
349	return 0;
350	}
351	else {
352	return sign;
353	}
354	}
355	if (x1 == 0.0) {
356	return -sign;
357	}
358	}
359
360	}
361
362	/**
363	* Returns the index of the direction of the point <code>q</code> relative to
364	* a vector specified by <code>p1-p2</code>.
365	*
366	* @param p1 the origin point of the vector
367	* @param p2 the final point of the vector
368	* @param q the point to compute the direction to
369	*
370	* @return 1 if q is counter-clockwise (left) from p1-p2
371	* @return -1 if q is clockwise (right) from p1-p2
372	* @return 0 if q is collinear with p1-p2
373	*/
374	public static int orientationIndex(Coordinate p1, Coordinate p2, Coordinate q)
375	{
376	/**
377	* MD - 9 Aug 2010 It seems that the basic algorithm is slightly orientation
378	* dependent, when computing the orientation of a point very close to a
379	* line. This is possibly due to the arithmetic in the translation to the
380	* origin.
381	*
382	* For instance, the following situation produces identical results in spite
383	* of the inverse orientation of the line segment:
384	*
385	* Coordinate p0 = new Coordinate(219.3649559090992, 140.84159161824724);
386	* Coordinate p1 = new Coordinate(168.9018919682399, -5.713787599646864);
387	*
388	* Coordinate p = new Coordinate(186.80814046338352, 46.28973405831556); int
389	* orient = orientationIndex(p0, p1, p); int orientInv =
390	* orientationIndex(p1, p0, p);
391	*
392	*
393	*/
394
395	double dx1 = p2.x - p1.x;
396	double dy1 = p2.y - p1.y;
397	double dx2 = q.x - p2.x;
398	double dy2 = q.y - p2.y;
399	return signOfDet2x2(dx1, dy1, dx2, dy2);
400	}
401
402	}
403