|
10 | 10 | // except according to those terms. |
11 | 11 | #![allow(non_snake_case)] |
12 | 12 |
|
13 | | -//! This is a direct transcript of the sourcecode and algorithms provided by |
| 13 | +//! # Adaptive Precision Floating-Point Arithmetic and Fast Robust Predicates for Computational Geometry |
| 14 | +//! This is a direct transcript of the source code and algorithms provided by |
14 | 15 | //! Jonathan Richard Shewchuk ([https://www.cs.cmu.edu/~quake/robust.html](https://www.cs.cmu.edu/~quake/robust.html)) |
15 | 16 | //! See the paper and the source code for more information. |
16 | 17 | //! |
17 | 18 | //! The module offers adaptive and precise calculations for orientation queries |
18 | | -//! (on which side of a line does a point lie?) and in-circle queries |
19 | | -//! (is a given point contained in the circumference of a triangle?) |
| 19 | +//! – "on which side of a line (2d) or plane (3d) does a point lie?" – and in-circle / in-sphere queries |
| 20 | +//! – "is a given point contained in the circumference of a triangle?". |
20 | 21 | //! The "adaptive" nature will increase performance only if a simpler calculation |
21 | 22 | //! cannot be guaranteed to be accurate enough, yielding higher performance on |
22 | 23 | //! average. |
@@ -64,7 +65,7 @@ const ISPERRBOUND_B: f64 = (5.0 + 72.0 * EPSILON) * EPSILON; |
64 | 65 | const ISPERRBOUND_C: f64 = (71.0 + 1408.0 * EPSILON) * EPSILON * EPSILON; |
65 | 66 |
|
66 | 67 | /// Returns a positive value if the coordinates `pa`, `pb`, and `pc` occur in counterclockwise order |
67 | | -/// (pc lies to the **left** of the directed line defined by coordinates pa and pb). |
| 68 | +/// (`pc` lies to the **left** of the directed line defined by coordinates `pa` and `pb`). |
68 | 69 | /// Returns a negative value if they occur in clockwise order (`pc` lies to the **right** of the directed line `pa, pb`). |
69 | 70 | /// Returns `0` if they are **collinear**. |
70 | 71 | pub fn orient2d<T: Into<f64>>(pa: Coord<T>, pb: Coord<T>, pc: Coord<T>) -> f64 { |
@@ -168,8 +169,8 @@ fn orient2dadapt(pa: Coord<f64>, pb: Coord<f64>, pc: Coord<f64>, detsum: f64) -> |
168 | 169 | } |
169 | 170 |
|
170 | 171 | /// Returns a positive value if the point `pd` lies below the plane passing through `pa`, `pb`, and `pc` |
171 | | -/// ("below" is defined so that pa, pb, and pc appear in counterclockwise order when viewed from above the plane). |
172 | | -/// Returns a negative value if `pd` lies above the plane |
| 172 | +/// ("below" is defined so that `pa`, `pb`, and `pc` appear in counterclockwise order when viewed from above the plane). |
| 173 | +/// Returns a negative value if `pd` lies above the plane. |
173 | 174 | /// Returns `0` if they are **coplanar**. |
174 | 175 | pub fn orient3d<T: Into<f64>>( |
175 | 176 | pa: Coord3D<T>, |
@@ -1325,10 +1326,10 @@ fn incircleadapt( |
1325 | 1326 | fin1[finlength - 1] |
1326 | 1327 | } |
1327 | 1328 |
|
1328 | | -/// Return a positive value if the point pe lies inside the sphere passing through `pa`, `pb`, `pc`, and `pd` |
1329 | | -/// Returns a negative value if it lies outside |
1330 | | -/// Returns `0` if the five points are **cospherical** |
1331 | | -/// The points pa, pb, pc, and pd must be ordered so that they have a positive orientation |
| 1329 | +/// Returns a positive value if the point `pe` lies inside the sphere passing through `pa`, `pb`, `pc`, and `pd`. |
| 1330 | +/// Returns a negative value if it lies outside. |
| 1331 | +/// Returns `0` if the five points are **cospherical**. |
| 1332 | +/// **NOTE**: The points `pa`, `pb`, `pc`, and `pd` must be ordered so that they have a positive orientation. |
1332 | 1333 | pub fn insphere<T: Into<f64>>( |
1333 | 1334 | pa: Coord3D<T>, |
1334 | 1335 | pb: Coord3D<T>, |
@@ -2240,8 +2241,237 @@ fn abs(x: f64) -> f64 { |
2240 | 2241 |
|
2241 | 2242 | #[cfg(test)] |
2242 | 2243 | mod test { |
| 2244 | + |
| 2245 | + // Test fixtures are from the Staged Geometric Predicates (2001) paper by Aleksandar Nanevski et al |
| 2246 | + // Fixtures copied from https://github.com/mourner/robust-predicates/tree/main/test/fixtures |
| 2247 | + // Original location: https://www.cs.cmu.edu/afs/cs/project/pscico/pscico/src/arithmetic/compiler1/test/ |
| 2248 | + |
2243 | 2249 | use super::{incircle, insphere, orient2d, orient3d, Coord, Coord3D}; |
2244 | 2250 |
|
| 2251 | + #[cfg(not(feature = "no_std"))] |
| 2252 | + use std::fs::File; |
| 2253 | + #[cfg(not(feature = "no_std"))] |
| 2254 | + use std::io::{self, Read}; |
| 2255 | + |
| 2256 | + #[cfg(not(feature = "no_std"))] |
| 2257 | + fn filename_to_string(s: &str) -> io::Result<String> { |
| 2258 | + let mut file = File::open(s)?; |
| 2259 | + let mut s = String::new(); |
| 2260 | + file.read_to_string(&mut s)?; |
| 2261 | + Ok(s) |
| 2262 | + } |
| 2263 | + |
| 2264 | + #[cfg(not(feature = "no_std"))] |
| 2265 | + fn results_by_line(s: &str) -> Vec<Vec<f64>> { |
| 2266 | + s.lines() |
| 2267 | + .map(|line| { |
| 2268 | + line.split_whitespace() |
| 2269 | + .skip(1) |
| 2270 | + .map(|foo| foo.parse::<f64>().unwrap()) |
| 2271 | + .collect() |
| 2272 | + }) |
| 2273 | + .collect() |
| 2274 | + } |
| 2275 | + |
| 2276 | + #[cfg(not(feature = "no_std"))] |
| 2277 | + #[test] |
| 2278 | + fn test_orient2d_fixtures() { |
| 2279 | + let f = filename_to_string("fixtures/orient2d.txt").unwrap(); |
| 2280 | + let fixtures = results_by_line(&f); |
| 2281 | + fixtures.iter().enumerate().for_each(|(idx, fixture)| { |
| 2282 | + let ax = fixture[0]; |
| 2283 | + let ay = fixture[1]; |
| 2284 | + let bx = fixture[2]; |
| 2285 | + let by = fixture[3]; |
| 2286 | + let cx = fixture[4]; |
| 2287 | + let cy = fixture[5]; |
| 2288 | + let sign = fixture[6]; |
| 2289 | + let c1 = Coord { x: ax, y: ay }; |
| 2290 | + let c2 = Coord { x: bx, y: by }; |
| 2291 | + let c3 = Coord { x: cx, y: cy }; |
| 2292 | + let res = orient2d(c1, c2, c3); |
| 2293 | + // result sign and fixture sign should be equal |
| 2294 | + assert!( |
| 2295 | + res.signum() == sign.signum(), |
| 2296 | + "Line {line}: Result sign ({result}) and fixture sign ({sign}) should match |
| 2297 | + \nCoord 1: {c1:?}\nCoord 2: {c2:?}\nCoord 3: {c3:?}", |
| 2298 | + line = idx + 1, |
| 2299 | + result = res, |
| 2300 | + sign = sign, |
| 2301 | + c1 = c1, |
| 2302 | + c2 = c2, |
| 2303 | + c3 = c3 |
| 2304 | + ); |
| 2305 | + }) |
| 2306 | + } |
| 2307 | + |
| 2308 | + #[cfg(not(feature = "no_std"))] |
| 2309 | + #[test] |
| 2310 | + fn test_incircle_fixtures() { |
| 2311 | + let f = filename_to_string("fixtures/incircle.txt").unwrap(); |
| 2312 | + let fixtures = results_by_line(&f); |
| 2313 | + fixtures.iter().enumerate().for_each(|(idx, fixture)| { |
| 2314 | + let ax = fixture[0]; |
| 2315 | + let ay = fixture[1]; |
| 2316 | + let bx = fixture[2]; |
| 2317 | + let by = fixture[3]; |
| 2318 | + let cx = fixture[4]; |
| 2319 | + let cy = fixture[5]; |
| 2320 | + let dx = fixture[6]; |
| 2321 | + let dy = fixture[7]; |
| 2322 | + let sign = fixture[8]; |
| 2323 | + let c1 = Coord { x: ax, y: ay }; |
| 2324 | + let c2 = Coord { x: bx, y: by }; |
| 2325 | + let c3 = Coord { x: cx, y: cy }; |
| 2326 | + let c4 = Coord { x: dx, y: dy }; |
| 2327 | + let res = incircle(c1, c2, c3, c4); |
| 2328 | + assert!( |
| 2329 | + res.signum() == sign.signum(), |
| 2330 | + "Line {line}: Result sign ({result}) and fixture sign ({sign}) should match |
| 2331 | + \nCoord 1: {c1:?}\nCoord 2: {c2:?}\nCoord 3: {c3:?}\nCoord 4: {c4:?}", |
| 2332 | + line = idx + 1, |
| 2333 | + result = res, |
| 2334 | + sign = sign, |
| 2335 | + c1 = c1, |
| 2336 | + c2 = c2, |
| 2337 | + c3 = c3, |
| 2338 | + c4 = c4 |
| 2339 | + ); |
| 2340 | + }) |
| 2341 | + } |
| 2342 | + |
| 2343 | + #[cfg(not(feature = "no_std"))] |
| 2344 | + #[test] |
| 2345 | + fn test_orient3d_fixtures() { |
| 2346 | + let f = filename_to_string("fixtures/orient3d.txt").unwrap(); |
| 2347 | + let fixtures = results_by_line(&f); |
| 2348 | + fixtures.iter().enumerate().for_each(|(idx, fixture)| { |
| 2349 | + let ax = fixture[0]; |
| 2350 | + let ay = fixture[1]; |
| 2351 | + let az = fixture[2]; |
| 2352 | + |
| 2353 | + let bx = fixture[3]; |
| 2354 | + let by = fixture[4]; |
| 2355 | + let bz = fixture[5]; |
| 2356 | + |
| 2357 | + let cx = fixture[6]; |
| 2358 | + let cy = fixture[7]; |
| 2359 | + let cz = fixture[8]; |
| 2360 | + |
| 2361 | + let dx = fixture[9]; |
| 2362 | + let dy = fixture[10]; |
| 2363 | + let dz = fixture[11]; |
| 2364 | + |
| 2365 | + let sign = fixture[12]; |
| 2366 | + |
| 2367 | + let c1 = Coord3D { |
| 2368 | + x: ax, |
| 2369 | + y: ay, |
| 2370 | + z: az, |
| 2371 | + }; |
| 2372 | + let c2 = Coord3D { |
| 2373 | + x: bx, |
| 2374 | + y: by, |
| 2375 | + z: bz, |
| 2376 | + }; |
| 2377 | + let c3 = Coord3D { |
| 2378 | + x: cx, |
| 2379 | + y: cy, |
| 2380 | + z: cz, |
| 2381 | + }; |
| 2382 | + let c4 = Coord3D { |
| 2383 | + x: dx, |
| 2384 | + y: dy, |
| 2385 | + z: dz, |
| 2386 | + }; |
| 2387 | + let res = orient3d(c1, c2, c3, c4); |
| 2388 | + assert!( |
| 2389 | + res.signum() == sign.signum(), |
| 2390 | + "Line {line}: Result sign ({result}) and fixture sign ({sign}) should match |
| 2391 | + \nCoord 1: {c1:?}\nCoord 2: {c2:?}\nCoord 3: {c3:?}\nCoord 4: {c4:?}", |
| 2392 | + line = idx + 1, |
| 2393 | + result = res, |
| 2394 | + sign = sign, |
| 2395 | + c1 = c1, |
| 2396 | + c2 = c2, |
| 2397 | + c3 = c3, |
| 2398 | + c4 = c4 |
| 2399 | + ); |
| 2400 | + // symmetry |
| 2401 | + assert!(res.signum() == orient3d(c1, c2, c3, c4).signum()); |
| 2402 | + }) |
| 2403 | + } |
| 2404 | + |
| 2405 | + #[cfg(not(feature = "no_std"))] |
| 2406 | + #[test] |
| 2407 | + fn test_insphere_fixtures() { |
| 2408 | + let f = filename_to_string("fixtures/insphere.txt").unwrap(); |
| 2409 | + let fixtures = results_by_line(&f); |
| 2410 | + fixtures.iter().enumerate().for_each(|(idx, fixture)| { |
| 2411 | + let ax = fixture[0]; |
| 2412 | + let ay = fixture[1]; |
| 2413 | + let az = fixture[2]; |
| 2414 | + |
| 2415 | + let bx = fixture[3]; |
| 2416 | + let by = fixture[4]; |
| 2417 | + let bz = fixture[5]; |
| 2418 | + |
| 2419 | + let cx = fixture[6]; |
| 2420 | + let cy = fixture[7]; |
| 2421 | + let cz = fixture[8]; |
| 2422 | + |
| 2423 | + let dx = fixture[9]; |
| 2424 | + let dy = fixture[10]; |
| 2425 | + let dz = fixture[11]; |
| 2426 | + |
| 2427 | + let ex = fixture[12]; |
| 2428 | + let ey = fixture[13]; |
| 2429 | + let ez = fixture[14]; |
| 2430 | + |
| 2431 | + let sign = fixture[15]; |
| 2432 | + |
| 2433 | + let c1 = Coord3D { |
| 2434 | + x: ax, |
| 2435 | + y: ay, |
| 2436 | + z: az, |
| 2437 | + }; |
| 2438 | + let c2 = Coord3D { |
| 2439 | + x: bx, |
| 2440 | + y: by, |
| 2441 | + z: bz, |
| 2442 | + }; |
| 2443 | + let c3 = Coord3D { |
| 2444 | + x: cx, |
| 2445 | + y: cy, |
| 2446 | + z: cz, |
| 2447 | + }; |
| 2448 | + let c4 = Coord3D { |
| 2449 | + x: dx, |
| 2450 | + y: dy, |
| 2451 | + z: dz, |
| 2452 | + }; |
| 2453 | + let c5 = Coord3D { |
| 2454 | + x: ex, |
| 2455 | + y: ey, |
| 2456 | + z: ez, |
| 2457 | + }; |
| 2458 | + let res = insphere(c1, c2, c3, c4, c5); |
| 2459 | + assert!( |
| 2460 | + res.signum() == sign.signum(), |
| 2461 | + "Line {line}: Result sign ({result}) and fixture sign ({sign}) should match |
| 2462 | + \nCoord 1: {c1:?}\nCoord 2: {c2:?}\nCoord 3: {c3:?}\nCoord 4: {c4:?}\nCoord 5: {c5:?}", |
| 2463 | + line = idx + 1, |
| 2464 | + result = res, |
| 2465 | + sign = sign, |
| 2466 | + c1 = c1, |
| 2467 | + c2 = c2, |
| 2468 | + c3 = c3, |
| 2469 | + c4 = c4, |
| 2470 | + c5 = c5 |
| 2471 | + ); |
| 2472 | + }) |
| 2473 | + } |
| 2474 | + |
2245 | 2475 | #[test] |
2246 | 2476 | fn test_orient2d() { |
2247 | 2477 | let from = Coord { x: -1f64, y: -1.0 }; |
|
0 commit comments