From c690361fbf24f63dc8c27a5a395c7f09a2406b2e Mon Sep 17 00:00:00 2001 From: ringsaturn Date: Fri, 23 Jun 2023 00:48:12 +0800 Subject: [PATCH] Exp: Add RTree index for polygon (#4) --- Cargo.lock | 18 +++- Cargo.toml | 3 +- README.md | 52 +---------- benches/bench.rs | 1 + benches/bench_az.rs | 24 ++++- examples/demo.rs | 44 --------- src/lib.rs | 221 ++++++++++++++++++++++++++++++++++++++------ 7 files changed, 240 insertions(+), 123 deletions(-) delete mode 100644 examples/demo.rs diff --git a/Cargo.lock b/Cargo.lock index bae05fb..373aa09 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -10,10 +10,11 @@ checksum = "8bf7cc16383c4b8d58b9905a8509f02926ce3058053c056376248d958c9df1e8" [[package]] name = "geometry-rs" -version = "0.1.2" +version = "0.2.0" dependencies = [ "float_next_after", "lazy_static", + "rtree_rs", "serde", "serde_json", ] @@ -30,6 +31,12 @@ version = "1.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" +[[package]] +name = "pqueue" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2145d14f09d5fc7fe7134b556146599c2929af5b1f1e3a0ecc9d582a9b85e4d" + [[package]] name = "proc-macro2" version = "1.0.60" @@ -48,6 +55,15 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "rtree_rs" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7787960e978c1c675fd6b8eb7d56b9c7162aec55c41d368add6f3ff39251b96e" +dependencies = [ + "pqueue", +] + [[package]] name = "ryu" version = "1.0.13" diff --git a/Cargo.toml b/Cargo.toml index 4d70372..11fff14 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-rs" -version = "0.1.2" +version = "0.2.0" edition = "2021" license-file = "LICENSE" description = "geometry utils" @@ -15,6 +15,7 @@ keywords = ["geometry"] [dependencies] float_next_after = "1.0.0" +rtree_rs = "0.1.4" [dev-dependencies] serde = { version = "1.0", features = ["derive"] } diff --git a/README.md b/README.md index bb4d6d9..d023d2e 100644 --- a/README.md +++ b/README.md @@ -3,54 +3,10 @@ Rewrite parts of [tidwall/geometry](https://github.com/tidwall/geometry) to Rust for [ringsaturn/tzf-rs](https://github.com/ringsaturn/tzf-rs). -```toml -[dependencies] -geometry-rs = "0.1.2" +```bash +cargo add geometry-rs ``` -```rust -use std::vec; +## TODO -use geometry_rs; - -fn main() { - let poly = geometry_rs::Polygon::new( - vec![ - geometry_rs::Point { - x: 90.48826291293898, - y: 45.951129815858565, - }, - geometry_rs::Point { - x: 90.48826291293898, - y: 27.99437617512571, - }, - geometry_rs::Point { - x: 122.83201291294, - y: 27.99437617512571, - }, - geometry_rs::Point { - x: 122.83201291294, - y: 45.951129815858565, - }, - geometry_rs::Point { - x: 90.48826291293898, - y: 45.951129815858565, - }, - ], - vec![], - ); - - let p_out = geometry_rs::Point { - x: 130.74216916294148, - y: 37.649011392900306, - }; - - print!("{:?}\n", poly.contains_point(p_out)); - - let p_in = geometry_rs::Point { - x: 99.9804504129416, - y: 39.70716466970461, - }; - print!("{:?}\n", poly.contains_point(p_in)); -} -``` +- [ ] Use compressed data of RTree index, since too many memory costs. diff --git a/benches/bench.rs b/benches/bench.rs index d380562..e696bcc 100644 --- a/benches/bench.rs +++ b/benches/bench.rs @@ -31,6 +31,7 @@ mod benches_polygon { }, ], vec![], + false, ); return poly; } diff --git a/benches/bench_az.rs b/benches/bench_az.rs index 5aeb713..f189761 100644 --- a/benches/bench_az.rs +++ b/benches/bench_az.rs @@ -5,7 +5,28 @@ mod benches_az_polygon { extern crate test; use test::Bencher; - // test data copy from https://github.com/unitedstates/districts/blob/gh-pages/states/AZ/shape.geojson + /* Test data copy from https://github.com/unitedstates/districts/blob/gh-pages/states/AZ/shape.geojson + + Python code to generate the data: + ```python + import json + + TPL = """geometry_rs::Point {x: {lng}, y: {lat}},""" + + with open("./az.geojson") as f: + data = json.loads(f.read()) + + coordinates = data["coordinates"][0][0] + # print(len(coordinates)) + gens = [] + for coord in coordinates: + # gens.append(TPL.format(lng=coord[0], lat=coord[1])) + gens.append(TPL.replace("{lng}", str(coord[0])).replace("{lat}", str(coord[1]))) + gen_text = "\n".join(gens) + print(gen_text) + ``` + */ + fn load_poly() -> geometry_rs::Polygon{ let poly = geometry_rs::Polygon::new( vec![ @@ -1594,6 +1615,7 @@ mod benches_az_polygon { ], vec![], + true, ); return poly; } diff --git a/examples/demo.rs b/examples/demo.rs deleted file mode 100644 index 4164add..0000000 --- a/examples/demo.rs +++ /dev/null @@ -1,44 +0,0 @@ -use std::vec; - -use geometry_rs; - -fn main() { - let poly = geometry_rs::Polygon::new( - vec![ - geometry_rs::Point { - x: 90.48826291293898, - y: 45.951129815858565, - }, - geometry_rs::Point { - x: 90.48826291293898, - y: 27.99437617512571, - }, - geometry_rs::Point { - x: 122.83201291294, - y: 27.99437617512571, - }, - geometry_rs::Point { - x: 122.83201291294, - y: 45.951129815858565, - }, - geometry_rs::Point { - x: 90.48826291293898, - y: 45.951129815858565, - }, - ], - vec![], - ); - - let p_out = geometry_rs::Point { - x: 130.74216916294148, - y: 37.649011392900306, - }; - - print!("{:?}\n", poly.contains_point(p_out)); - - let p_in = geometry_rs::Point { - x: 99.9804504129416, - y: 39.70716466970461, - }; - print!("{:?}\n", poly.contains_point(p_in)); -} diff --git a/src/lib.rs b/src/lib.rs index f6f9472..7d57c73 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,8 @@ +#![doc = include_str!("../README.md")] + use float_next_after::NextAfter; +use rtree_rs::RTree; +use rtree_rs::Rect as RTreeRect; #[derive(Copy, Clone, Debug)] pub struct Point { @@ -106,7 +110,7 @@ fn segment_at_for_vec_point(exterior: &Vec, index: i64) -> Segment { return Segment { a: seg_a, b: seg_b }; } -fn rins_contains_point(ring: &Vec, point: Point, allow_on_edge: bool) -> bool { +fn rings_contains_point(ring: &Vec, point: Point, allow_on_edge: bool) -> bool { let rect = Rect { min: Point { x: std::f64::NEG_INFINITY, @@ -118,11 +122,12 @@ fn rins_contains_point(ring: &Vec, point: Point, allow_on_edge: bool) -> }, }; let mut inside: bool = false; - let n = (ring.len() - 1) as i64; + let n: i64 = (ring.len() - 1) as i64; for i in 0..n { - let seg = segment_at_for_vec_point(&ring, i); + let seg: Segment = segment_at_for_vec_point(&ring, i); + if seg.rect().intersects_rect(rect) { - let res = raycast(&seg, point); + let res: RaycastResult = raycast(&seg, point); // print!("res= inside:{:?} on:{:?}\n", res.inside, res.on); if res.on { inside = allow_on_edge; @@ -136,26 +141,84 @@ fn rins_contains_point(ring: &Vec, point: Point, allow_on_edge: bool) -> return inside; } -#[derive(Clone, Debug)] +fn rings_contains_point_by_rtree_index( + ring: &Vec, + ring_rtree: &rtree_rs::RTree<2, f64, i64>, + point: Point, + allow_on_edge: bool, +) -> bool { + let rect = Rect { + min: Point { + x: std::f64::NEG_INFINITY, + y: point.y, + }, + max: Point { + x: std::f64::INFINITY, + y: point.y, + }, + }; + for item in ring_rtree.search(RTreeRect::new( + [std::f64::NEG_INFINITY, point.y], + [std::f64::INFINITY, point.y], + )) { + let seg: Segment = segment_at_for_vec_point(&ring, *item.data); + let irect = seg.rect(); + if irect.intersects_rect(rect) { + let res: RaycastResult = raycast(&seg, point); + if res.on { + return allow_on_edge; + } + if res.inside { + return true; + } + } + } + return false; +} + pub struct Polygon { exterior: Vec, + exterior_rtree: rtree_rs::RTree<2, f64, i64>, holes: Vec>, + holes_rtree: Vec>, rect: Rect, + with_index: bool, } impl Polygon { - pub fn contains_point(&self, p: Point) -> bool { - if !self.rect.contains_point(p) { + /// Point-In-Polygon check with RTree index. + /// `with_index` param must true for this query. + /// Query speed is faster compare with [contains_point_normal]. + fn contains_point_with_index(&self, p: Point) -> bool { + if !rings_contains_point_by_rtree_index(&self.exterior, &self.exterior_rtree, p, false) { return false; } - if !rins_contains_point(&self.exterior, p, false) { - return false; + let mut contains: bool = true; + let mut i: usize = 0; + for hole in self.holes.iter() { + let tr = self.holes_rtree.get(i).unwrap(); + if rings_contains_point_by_rtree_index(&hole, &tr, p, false) { + contains = false; + break; + } + + i += 1; } + return contains; + } + /// Point-In-Polygon check, the normal way. + /// It's most used algorithm implementation, port from Go's [geojson] + /// + /// [geojson]: https://github.com/tidwall/geojson + fn contains_point_normal(&self, p: Point) -> bool { + if !rings_contains_point(&self.exterior, p, false) { + return false; + } let mut contains: bool = true; for hole in self.holes.iter() { - if rins_contains_point(&hole, p, false) { + if rings_contains_point(&hole, p, false) { contains = false; break; } @@ -163,7 +226,70 @@ impl Polygon { return contains; } - pub fn new(exterior: Vec, holes: Vec>) -> Polygon { + /// Do point-in-polygon search. + pub fn contains_point(&self, p: Point) -> bool { + if !self.rect.contains_point(p) { + return false; + } + if self.with_index { + return self.contains_point_with_index(p); + } + return self.contains_point_normal(p); + } + + /// Create a new Polygon instance from exterior and holes. + /// + /// Please note that set `with_index` to true will increase performance, but requires more memory. + /// See [#4] for more details. + /// + /// [#4]: https://github.com/ringsaturn/geometry-rs/pull/4 + /// + /// Example: + /// + /// ```rust + /// use std::vec; + /// use geometry_rs; + /// let poly = geometry_rs::Polygon::new( + /// vec![ + /// geometry_rs::Point { + /// x: 90.48826291293898, + /// y: 45.951129815858565, + /// }, + /// geometry_rs::Point { + /// x: 90.48826291293898, + /// y: 27.99437617512571, + /// }, + /// geometry_rs::Point { + /// x: 122.83201291294, + /// y: 27.99437617512571, + /// }, + /// geometry_rs::Point { + /// x: 122.83201291294, + /// y: 45.951129815858565, + /// }, + /// geometry_rs::Point { + /// x: 90.48826291293898, + /// y: 45.951129815858565, + /// }, + /// ], + /// vec![], + /// false, + /// ); + /// + /// let p_out = geometry_rs::Point { + /// x: 130.74216916294148, + /// y: 37.649011392900306, + /// }; + /// + /// print!("{:?}\n", poly.contains_point(p_out)); + /// + /// let p_in = geometry_rs::Point { + /// x: 99.9804504129416, + /// y: 39.70716466970461, + /// }; + /// print!("{:?}\n", poly.contains_point(p_in)); + /// ``` + pub fn new(exterior: Vec, holes: Vec>, with_index: bool) -> Polygon { let mut minx: f64 = exterior.get(0).unwrap().x; let mut miny: f64 = exterior.get(0).unwrap().y; let mut maxx: f64 = exterior.get(0).unwrap().x; @@ -191,10 +317,49 @@ impl Polygon { max: Point { x: maxx, y: maxy }, }; + let mut exterior_rtree = RTree::new(); + let n = (exterior.len() - 1) as i64; + for i in 0..n { + let segrect = segment_at_for_vec_point(&exterior, i).rect(); + if with_index { + exterior_rtree.insert( + RTreeRect::new( + [segrect.min.x, segrect.min.y], + [segrect.max.x, segrect.max.y], + ), + i as i64, + ); + } + } + + let mut holes_rtree = vec![]; + for hole_poly in holes.iter() { + let mut hole_rtre = RTree::new(); + let n = (hole_poly.len() - 1) as i64; + for i in 0..n { + let segrect = segment_at_for_vec_point(&hole_poly, i).rect(); + if with_index { + hole_rtre.insert( + RTreeRect::new( + [segrect.min.x, segrect.min.y], + [segrect.max.x, segrect.max.y], + ), + i as i64, + ); + } + } + if with_index { + holes_rtree.push(hole_rtre); + } + } + return Polygon { exterior, + exterior_rtree, holes, + holes_rtree, rect, + with_index, }; } } @@ -207,28 +372,28 @@ pub struct Segment { impl Segment { pub fn rect(&self) -> Rect { - let mut minx: f64 = self.a.x; - let mut miny: f64 = self.a.y; - let mut maxx: f64 = self.b.x; - let mut maxy: f64 = self.b.y; - - if minx > maxx { - let actualminx = maxx; - let actualmaxx = minx; - minx = actualminx; - maxx = actualmaxx; + let mut min_x: f64 = self.a.x; + let mut min_y: f64 = self.a.y; + let mut max_x: f64 = self.b.x; + let mut max_y: f64 = self.b.y; + + if min_x > max_x { + let actual_min_x = max_x; + let actual_max_x = min_x; + min_x = actual_min_x; + max_x = actual_max_x; } - if miny > maxy { - let actualminy = maxy; - let actualmaxy = miny; - miny = actualminy; - maxy = actualmaxy; + if min_y > max_y { + let actual_min_y = max_y; + let actual_max_y = min_y; + min_y = actual_min_y; + max_y = actual_max_y; } return Rect { - min: Point { x: minx, y: miny }, - max: Point { x: maxx, y: maxy }, + min: Point { x: min_x, y: min_y }, + max: Point { x: max_x, y: max_y }, }; } }