gistools/geometry/predicates/
orient2d.rs

1use super::util::{
2    CCWERRBOUND_A, CCWERRBOUND_B, CCWERRBOUND_C, RESULTERRBOUND, SPLITTER, estimate, pred_sum,
3};
4use libm::fabs;
5
6/// a negative value if the points a, b, and c occur in counterclockwise order
7/// (c lies to the left of the directed line defined by points a and b).
8///
9/// ## Returns
10/// - a positive value if they occur in clockwise order (c lies to the right of the directed line ab).
11/// - zero if they are collinear.
12pub fn orient2dadapt(ax: f64, ay: f64, bx: f64, by: f64, cx: f64, cy: f64, detsum: f64) -> f64 {
13    let mut b: [f64; 4] = [0.0; 4];
14    let mut c1: [f64; 8] = [0.0; 8];
15    let mut c2: [f64; 12] = [0.0; 12];
16    let mut d: [f64; 16] = [0.0; 16];
17    let mut u: [f64; 4] = [0.0; 4];
18    let mut bvirt: f64;
19    let mut c: f64;
20    let mut ahi: f64;
21    let mut alo: f64;
22    let mut bhi: f64;
23    let mut blo: f64;
24    let mut _i: f64;
25    let mut _j: f64;
26    let mut _z: f64;
27    let mut s1: f64;
28    let mut s0: f64;
29    let mut t1: f64;
30    let mut t0: f64;
31    let mut u3: f64;
32
33    let acx = ax - cx;
34    let bcx = bx - cx;
35    let acy = ay - cy;
36    let bcy = by - cy;
37
38    s1 = acx * bcy;
39    c = SPLITTER * acx;
40    ahi = c - (c - acx);
41    alo = acx - ahi;
42    c = SPLITTER * bcy;
43    bhi = c - (c - bcy);
44    blo = bcy - bhi;
45    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
46    t1 = acy * bcx;
47    c = SPLITTER * acy;
48    ahi = c - (c - acy);
49    alo = acy - ahi;
50    c = SPLITTER * bcx;
51    bhi = c - (c - bcx);
52    blo = bcx - bhi;
53    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
54    _i = s0 - t0;
55    bvirt = s0 - _i;
56    b[0] = s0 - (_i + bvirt) + (bvirt - t0);
57    _j = s1 + _i;
58    bvirt = _j - s1;
59    _z = s1 - (_j - bvirt) + (_i - bvirt);
60    _i = _z - t1;
61    bvirt = _z - _i;
62    b[1] = _z - (_i + bvirt) + (bvirt - t1);
63    u3 = _j + _i;
64    bvirt = u3 - _j;
65    b[2] = _j - (u3 - bvirt) + (_i - bvirt);
66    b[3] = u3;
67
68    let mut det = estimate(4, &b);
69    let mut errbound = CCWERRBOUND_B * detsum;
70    if det >= errbound || -det >= errbound {
71        return det;
72    }
73
74    bvirt = ax - acx;
75    let acxtail = ax - (acx + bvirt) + (bvirt - cx);
76    bvirt = bx - bcx;
77    let bcxtail = bx - (bcx + bvirt) + (bvirt - cx);
78    bvirt = ay - acy;
79    let acytail = ay - (acy + bvirt) + (bvirt - cy);
80    bvirt = by - bcy;
81    let bcytail = by - (bcy + bvirt) + (bvirt - cy);
82
83    if acxtail == 0. && acytail == 0. && bcxtail == 0. && bcytail == 0. {
84        return det;
85    }
86
87    errbound = CCWERRBOUND_C * detsum + RESULTERRBOUND * fabs(det);
88    det += acx * bcytail + bcy * acxtail - (acy * bcxtail + bcx * acytail);
89    if det >= errbound || -det >= errbound {
90        return det;
91    }
92
93    s1 = acxtail * bcy;
94    c = SPLITTER * acxtail;
95    ahi = c - (c - acxtail);
96    alo = acxtail - ahi;
97    c = SPLITTER * bcy;
98    bhi = c - (c - bcy);
99    blo = bcy - bhi;
100    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
101    t1 = acytail * bcx;
102    c = SPLITTER * acytail;
103    ahi = c - (c - acytail);
104    alo = acytail - ahi;
105    c = SPLITTER * bcx;
106    bhi = c - (c - bcx);
107    blo = bcx - bhi;
108    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
109    _i = s0 - t0;
110    bvirt = s0 - _i;
111    u[0] = s0 - (_i + bvirt) + (bvirt - t0);
112    _j = s1 + _i;
113    bvirt = _j - s1;
114    _z = s1 - (_j - bvirt) + (_i - bvirt);
115    _i = _z - t1;
116    bvirt = _z - _i;
117    u[1] = _z - (_i + bvirt) + (bvirt - t1);
118    u3 = _j + _i;
119    bvirt = u3 - _j;
120    u[2] = _j - (u3 - bvirt) + (_i - bvirt);
121    u[3] = u3;
122    let c1len = pred_sum(4, &b, 4, &u, &mut c1);
123
124    s1 = acx * bcytail;
125    c = SPLITTER * acx;
126    ahi = c - (c - acx);
127    alo = acx - ahi;
128    c = SPLITTER * bcytail;
129    bhi = c - (c - bcytail);
130    blo = bcytail - bhi;
131    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
132    t1 = acy * bcxtail;
133    c = SPLITTER * acy;
134    ahi = c - (c - acy);
135    alo = acy - ahi;
136    c = SPLITTER * bcxtail;
137    bhi = c - (c - bcxtail);
138    blo = bcxtail - bhi;
139    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
140    _i = s0 - t0;
141    bvirt = s0 - _i;
142    u[0] = s0 - (_i + bvirt) + (bvirt - t0);
143    _j = s1 + _i;
144    bvirt = _j - s1;
145    _z = s1 - (_j - bvirt) + (_i - bvirt);
146    _i = _z - t1;
147    bvirt = _z - _i;
148    u[1] = _z - (_i + bvirt) + (bvirt - t1);
149    u3 = _j + _i;
150    bvirt = u3 - _j;
151    u[2] = _j - (u3 - bvirt) + (_i - bvirt);
152    u[3] = u3;
153    let c2len = pred_sum(c1len, &c1, 4, &u, &mut c2);
154
155    s1 = acxtail * bcytail;
156    c = SPLITTER * acxtail;
157    ahi = c - (c - acxtail);
158    alo = acxtail - ahi;
159    c = SPLITTER * bcytail;
160    bhi = c - (c - bcytail);
161    blo = bcytail - bhi;
162    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
163    t1 = acytail * bcxtail;
164    c = SPLITTER * acytail;
165    ahi = c - (c - acytail);
166    alo = acytail - ahi;
167    c = SPLITTER * bcxtail;
168    bhi = c - (c - bcxtail);
169    blo = bcxtail - bhi;
170    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
171    _i = s0 - t0;
172    bvirt = s0 - _i;
173    u[0] = s0 - (_i + bvirt) + (bvirt - t0);
174    _j = s1 + _i;
175    bvirt = _j - s1;
176    _z = s1 - (_j - bvirt) + (_i - bvirt);
177    _i = _z - t1;
178    bvirt = _z - _i;
179    u[1] = _z - (_i + bvirt) + (bvirt - t1);
180    u3 = _j + _i;
181    bvirt = u3 - _j;
182    u[2] = _j - (u3 - bvirt) + (_i - bvirt);
183    u[3] = u3;
184    let dlen = pred_sum(c2len, &c2, 4, &u, &mut d);
185
186    d[dlen - 1]
187}
188
189/// returns a positive value if the points a, b, and c occur in counterclockwise order
190/// (c lies to the left of the directed line defined by points a and b).
191///
192/// ## Returns
193/// - a negative value if they occur in clockwise order (c lies to the right of the directed line ab).
194/// - zero if they are collinear.
195pub fn orient2d(ax: f64, ay: f64, bx: f64, by: f64, cx: f64, cy: f64) -> f64 {
196    let detleft = (ay - cy) * (bx - cx);
197    let detright = (ax - cx) * (by - cy);
198    let det = detleft - detright;
199
200    let detsum = fabs(detleft + detright);
201    if fabs(det) >= CCWERRBOUND_A * detsum {
202        return det;
203    }
204
205    -orient2dadapt(ax, ay, bx, by, cx, cy, detsum)
206}
207
208/// returns a positive value if the points a, b, and c occur in counterclockwise order
209/// (c lies to the left of the directed line defined by points a and b).
210///
211/// ## Returns
212/// - a negative value if they occur in clockwise order (c lies to the right of the directed line ab).
213/// - zero if they are collinear.
214pub fn orient2dfast(ax: f64, ay: f64, bx: f64, by: f64, cx: f64, cy: f64) -> f64 {
215    (ay - cy) * (bx - cx) - (ax - cx) * (by - cy)
216}
gistools/geometry/predicates/orient2d.rs

gistools/geometry/predicates/
orient2d.rs
