From 4ebddfff71ddff381764dbc9da56aa9d20419bbf Mon Sep 17 00:00:00 2001 From: Timo Millenaar Date: Sat, 23 Mar 2024 22:02:22 +0100 Subject: [PATCH 1/5] Dirty version of line intersection on rect grids --- Cargo.toml | 1 + src/interpolate.rs | 4 +- src/lib.rs | 11 +++++ src/rect_grid.rs | 109 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 123 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index ae621c312..6934f35f2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,6 +9,7 @@ name = "gridkit_rs" crate-type = ["cdylib"] [dependencies] +geo = "0.28.0" geo-types = "0.7.12" numpy = "0.19.0" pyo3 = "0.19.0" diff --git a/src/interpolate.rs b/src/interpolate.rs index 2a17f0f41..d71a02b98 100644 --- a/src/interpolate.rs +++ b/src/interpolate.rs @@ -10,7 +10,7 @@ pub fn vec_norm_1d( return squared_sum.powf(0.5) } -fn _project( +pub fn project_point_on_line( point: &Array1, line_point_1: &Array1, line_point_2: &Array1 @@ -44,7 +44,7 @@ pub fn linear_interp_weights_single_triangle( median = &midpoint_opposite_side - &p1; } - let projected: ArrayBase, Dim<[usize; 1]>> = _project( + let projected: ArrayBase, Dim<[usize; 1]>> = project_point_on_line( &(sample_point - &p1), &(&p2 - &p1), &(&p3 - &p1), diff --git a/src/lib.rs b/src/lib.rs index e9ba05d0c..0b3b3b067 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -253,6 +253,17 @@ impl PyRectGrid { .cells_near_point(&point.as_array()) .into_pyarray(py) } + + fn cells_intersecting_line<'py>( + &self, + py: Python<'py>, + p1: PyReadonlyArray1<'py, f64>, + p2: PyReadonlyArray1<'py, f64>, + ) -> &'py PyArray2 { + self._grid + .cells_intersecting_line(&p1.as_array(), &p2.as_array()) + .into_pyarray(py) + } } #[pyclass] diff --git a/src/rect_grid.rs b/src/rect_grid.rs index dc4dd940c..86c31d566 100644 --- a/src/rect_grid.rs +++ b/src/rect_grid.rs @@ -1,5 +1,9 @@ use numpy::ndarray::*; use crate::utils::*; +use crate::interpolate::*; + +use geo_types::{LineString, Point, Coord, Geometry}; +use geo::Intersects; pub struct RectGrid { pub _dx: f64, @@ -172,4 +176,109 @@ impl RectGrid { } nearby_cells } + + pub fn cells_intersecting_line(&self, p1: &ArrayView1, p2: &ArrayView1) -> Array2 { + let mut ids = Array2::::zeros((0, 2)); + let mut cell_id = self.cell_at_point(&p1.into_shape((1, p1.len())).unwrap()); + let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)]].view()); + // let sides = array![ + // [[corners[Ix3(0,0,0)], corners[Ix3(0,0,1)]],[corners[Ix3(0,1,0)],corners[Ix3(0,1,1)]]], + // [[corners[Ix3(0,1,0)], corners[Ix3(0,1,1)]],[corners[Ix3(0,2,0)],corners[Ix3(0,2,1)]]], + // [[corners[Ix3(0,2,0)], corners[Ix3(0,2,1)]],[corners[Ix3(0,3,0)],corners[Ix3(0,3,1)]]], + // [[corners[Ix3(0,3,0)], corners[Ix3(0,3,1)]],[corners[Ix3(0,0,0)],corners[Ix3(0,0,1)]]], + // ]; + // for i in 0..corners.shape()[0] { + + let point1 = Coord:: {x:p1[Ix1(0)], y:p1[Ix1(1)]}; + let point2 = Coord:: {x:p2[Ix1(0)], y:p2[Ix1(1)]}; + let line = LineString::new(vec![point1, point2]); + + let mut intersection_id: i64 = -1; + let mut counter = 0; + loop { + let corners = self.cell_corners(&cell_id.view()); + let sides = vec![ + LineString::new(vec![Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]},Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]},Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]},Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]},Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]}]), + ]; + + let mut entered_through_side_id: i64; + match intersection_id { // convert intersection_id of previous cell to id of same side for new cell_id + 0 => {entered_through_side_id = 2;} + 1 => {entered_through_side_id = 3;} + 2 => {entered_through_side_id = 0;} + 3 => {entered_through_side_id = 1;} + _ => {entered_through_side_id = -1;} + } + + intersection_id = -1; + + for i in 0..sides.len() { + // TODO: Handle case where line intersects cell corner + println!("{},{}", entered_through_side_id, intersection_id); + if i as i64 == entered_through_side_id { // Discount the intersection towards previous cell + continue; + } + let intersects = sides[i].intersects(&line); + if intersects { + intersection_id = i as i64; + break; + } + + } + println!("{:?}, {}", cell_id, intersection_id); + match intersection_id { + 0 => { // Bottom side + cell_id[Ix2(0,1)] -= 1; + // cell_id = array![[cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] - 1]]; + // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] - 1].view()); + } + 1 => { // Right side + cell_id[Ix2(0,0)] += 1; + // cell_id = array![[cell_id[Ix2(0,0)] + 1, cell_id[Ix2(0,1)]]]; + // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)] + 1, cell_id[Ix2(0,1)]].view()); + } + 2 => { // Top side + cell_id[Ix2(0,1)] += 1; + // cell_id = array![[cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] + 1]]; + // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] + 1].view()); + } + 3 => { // Left side + cell_id[Ix2(0,0)] -= 1; + // cell_id = array![[cell_id[Ix2(0,0)] - 1, cell_id[Ix2(0,1)]]]; + // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)] - 1, cell_id[Ix2(0,1)]].view()); + } + _ => { // No intersection + // Reached the end of the line, break infinite loop and return from function + break; + } + } + let _ = ids.push(Axis(0), cell_id.slice(s![0, ..]).view()); + counter += 1; + if counter == 50 { + break + } + } + + // TODO: Remove previous intersection from sides before checking new intersection + // TODO: Loop until no intersection + + + + // let projected_vec_1 = projected_1 - &sides.slice(s![i, 0, ..]); + // let projected_vec_2 = projected_2 - &sides.slice(s![i, 1, ..]); + // // Check if the direction of the projected arrow is the same + // let same_x = (projected_vec_1[Ix1(0)] > 0.) == (projected_vec_2[Ix1(0)] > 0.); + // let same_y = (projected_vec_1[Ix1(1)] > 0.) == (projected_vec_2[Ix1(1)] > 0.); + // if (!same_x || !same_y) { + // // line crossing cell_side + // println!("{:?}, {:?}", projected_vec_1, projected_vec_2); + // } + // for points in sides.axis_iter(Axis(0)){} + // ids.into() + ids + // let mut ids_arr = Array2::::zeros((ids.shape()[0], 2)); + } } From e7110133e7a25321e3a7501ab5d8a3bca0ee9c3d Mon Sep 17 00:00:00 2001 From: Timo Millenaar Date: Sat, 30 Mar 2024 13:48:44 +0100 Subject: [PATCH 2/5] Fix RectGrid.cells_intersecting_line in rect_grid.rs when line is crossing corners as well as sides of cells --- src/rect_grid.rs | 168 ++++++++++++++++++++++++++++++----------------- 1 file changed, 108 insertions(+), 60 deletions(-) diff --git a/src/rect_grid.rs b/src/rect_grid.rs index 86c31d566..ce69ecf31 100644 --- a/src/rect_grid.rs +++ b/src/rect_grid.rs @@ -178,107 +178,155 @@ impl RectGrid { } pub fn cells_intersecting_line(&self, p1: &ArrayView1, p2: &ArrayView1) -> Array2 { + // This function returns the ids of the cells that intersect the line defined by the outer points p1 and p2. + // This is done through an infinite loop where we start at the cell that contains p1, + // check if any of the line intersects any of the corners or sides of the cell and depending + // on which side or corner is intersected, we find the next cell and repeat the process. + // Since the line also has an intersection where it 'entered' the new cell, + // we ignore this intersection and look for anohter. + // The loop is terminated when the line does not intersect any new corners or lines. + // Since the corners are also part of the lines as far as the line.intersects function is concerned, + // we check the corners first and skip the lines check if a corner is intersected. + // We then also have to ignore the lines that contain the corner when checking intersections for the next cell. + // + // The layout of the corners and lines with their indices are counter-clockwise starting at the bottom-left: + // + // C3 -- L2 -- C2 + // | | + // L3 L1 + // | | + // C0 -- L0 -- C1 + // + // Where C0 represents the first corner (corner index 0), C1 represents corner index 1 and so forth. + // The sampe applies to the lines, where L0 is the first line. + // let mut ids = Array2::::zeros((0, 2)); let mut cell_id = self.cell_at_point(&p1.into_shape((1, p1.len())).unwrap()); let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)]].view()); - // let sides = array![ - // [[corners[Ix3(0,0,0)], corners[Ix3(0,0,1)]],[corners[Ix3(0,1,0)],corners[Ix3(0,1,1)]]], - // [[corners[Ix3(0,1,0)], corners[Ix3(0,1,1)]],[corners[Ix3(0,2,0)],corners[Ix3(0,2,1)]]], - // [[corners[Ix3(0,2,0)], corners[Ix3(0,2,1)]],[corners[Ix3(0,3,0)],corners[Ix3(0,3,1)]]], - // [[corners[Ix3(0,3,0)], corners[Ix3(0,3,1)]],[corners[Ix3(0,0,0)],corners[Ix3(0,0,1)]]], - // ]; - // for i in 0..corners.shape()[0] { + // Create a LineString from the supplied endpoints let point1 = Coord:: {x:p1[Ix1(0)], y:p1[Ix1(1)]}; let point2 = Coord:: {x:p2[Ix1(0)], y:p2[Ix1(1)]}; let line = LineString::new(vec![point1, point2]); - let mut intersection_id: i64 = -1; - let mut counter = 0; + // TODO: Check if the line starts on a cell corner. If so, find out which of the connecting cells is the true starting cell. + + let mut side_intersection_id: i64 = -1; + let mut corner_intersection_id: i64 = -1; + let mut skip_corners = array![false, false, false, false]; + let mut skip_sides = array![false, false, false, false]; loop { let corners = self.cell_corners(&cell_id.view()); + + corner_intersection_id = -1; + for i in 0..4 { // Loop over corners + if skip_corners[i] { // Discount the intersection towards previous cell + println!("Skipping corner {}", i); + continue; + } + let intersects = line.intersects(&Coord:: {x:corners[Ix3(0,i,0)], y:corners[Ix3(0,i,1)]}); + if intersects { + corner_intersection_id = i as i64; + side_intersection_id = -1; // Reset + break; + } + } + + match corner_intersection_id { + // Adjust the cell-id to reflect the next cell and mark the oppisite corner and connecting sides to be skipped. + // To demonstrate what I mean with skipping the opposite corner and connecting sides: + // If the line now intersects the top-right corner, from the perspective of the next cell + // the line intersects the bottom-left corner, which is the one we want to ignore in the next iteration. + // This corner also belongs to the bottom and left sides, which we also want to ignore in the next iteration. + 0 => { // Bottom-left corner + cell_id[Ix2(0,0)] -= 1; + cell_id[Ix2(0,1)] -= 1; + skip_corners = array![false, false, true, false]; + skip_sides = array![false, true, true, false]; + } + 1 => { // Bottom-right corner + cell_id[Ix2(0,0)] += 1; + cell_id[Ix2(0,1)] -= 1; + skip_corners = array![false, false, false, true]; + skip_sides = array![false, false, true, true]; + } + 2 => { // Top-right corner + cell_id[Ix2(0,0)] += 1; + cell_id[Ix2(0,1)] += 1; + skip_corners = array![true, false, false, false]; + skip_sides = array![true, false, false, true]; + } + 3 => { // Top-left corner + cell_id[Ix2(0,0)] -= 1; + cell_id[Ix2(0,1)] += 1; + skip_corners = array![false, true, false, false]; + skip_sides = array![true, true, false, false]; + } + _ => {} // No intersection, check sides next + } + + // Add previous cell to vec and don't bother checking side intersections if we have a corner intersection + if corner_intersection_id != -1 { + println!("Crossed corner {}", corner_intersection_id); + let _ = ids.push(Axis(0), cell_id.slice(s![0, ..]).view()); + continue + } + + // Since there is no corner intersection, reset skip_corners + skip_corners = array![false, false, false, false]; + + // Check insersection on sides let sides = vec![ LineString::new(vec![Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]},Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]}]), LineString::new(vec![Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]},Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]}]), LineString::new(vec![Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]},Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]}]), LineString::new(vec![Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]},Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]}]), ]; - - let mut entered_through_side_id: i64; - match intersection_id { // convert intersection_id of previous cell to id of same side for new cell_id - 0 => {entered_through_side_id = 2;} - 1 => {entered_through_side_id = 3;} - 2 => {entered_through_side_id = 0;} - 3 => {entered_through_side_id = 1;} - _ => {entered_through_side_id = -1;} - } - - intersection_id = -1; + side_intersection_id = -1; for i in 0..sides.len() { - // TODO: Handle case where line intersects cell corner - println!("{},{}", entered_through_side_id, intersection_id); - if i as i64 == entered_through_side_id { // Discount the intersection towards previous cell + if skip_sides[i] { // Discount the intersection towards previous cell + println!("Skipping side {}", i); continue; } let intersects = sides[i].intersects(&line); if intersects { - intersection_id = i as i64; + side_intersection_id = i as i64; + corner_intersection_id = -1; // Also reset the corner break; } - } - println!("{:?}, {}", cell_id, intersection_id); - match intersection_id { + + match side_intersection_id { + // Adjust the cell-id to reflect the next cell and mark the oppisite side to be skipped. + // To demonstrate what I mean with skipping the opposite side: + // If the line now intersects the top side, from the perspective of the next cell + // the line intersects the bottom side, which is the one we want to ignore in the next iteration. 0 => { // Bottom side cell_id[Ix2(0,1)] -= 1; - // cell_id = array![[cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] - 1]]; - // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] - 1].view()); + skip_sides = array![false, false, true, false]; } 1 => { // Right side cell_id[Ix2(0,0)] += 1; - // cell_id = array![[cell_id[Ix2(0,0)] + 1, cell_id[Ix2(0,1)]]]; - // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)] + 1, cell_id[Ix2(0,1)]].view()); + skip_sides = array![false, false, false, true]; } 2 => { // Top side cell_id[Ix2(0,1)] += 1; - // cell_id = array![[cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] + 1]]; - // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)] + 1].view()); + skip_sides = array![true, false, false, false]; } 3 => { // Left side cell_id[Ix2(0,0)] -= 1; - // cell_id = array![[cell_id[Ix2(0,0)] - 1, cell_id[Ix2(0,1)]]]; - // let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)] - 1, cell_id[Ix2(0,1)]].view()); + skip_sides = array![false, true, false, false]; } _ => { // No intersection - // Reached the end of the line, break infinite loop and return from function + // Reached the end point of the line, break infinite loop and return from function break; } } + println!("Crossed side {}", side_intersection_id); let _ = ids.push(Axis(0), cell_id.slice(s![0, ..]).view()); - counter += 1; - if counter == 50 { - break - } } - - // TODO: Remove previous intersection from sides before checking new intersection - // TODO: Loop until no intersection - - - - // let projected_vec_1 = projected_1 - &sides.slice(s![i, 0, ..]); - // let projected_vec_2 = projected_2 - &sides.slice(s![i, 1, ..]); - // // Check if the direction of the projected arrow is the same - // let same_x = (projected_vec_1[Ix1(0)] > 0.) == (projected_vec_2[Ix1(0)] > 0.); - // let same_y = (projected_vec_1[Ix1(1)] > 0.) == (projected_vec_2[Ix1(1)] > 0.); - // if (!same_x || !same_y) { - // // line crossing cell_side - // println!("{:?}, {:?}", projected_vec_1, projected_vec_2); - // } - // for points in sides.axis_iter(Axis(0)){} - // ids.into() - ids - // let mut ids_arr = Array2::::zeros((ids.shape()[0], 2)); - } + return ids + } + } From 182fb119d2bb558ee405f6dab9a4b23c92a6c929 Mon Sep 17 00:00:00 2001 From: Timo Millenaar Date: Sat, 30 Mar 2024 15:26:33 +0100 Subject: [PATCH 3/5] Fix RectGrid.cells_intersecting_line in rect_grid.rs when line is starting or ending on a cell corner or cell side --- src/rect_grid.rs | 123 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 111 insertions(+), 12 deletions(-) diff --git a/src/rect_grid.rs b/src/rect_grid.rs index ce69ecf31..c130f38b1 100644 --- a/src/rect_grid.rs +++ b/src/rect_grid.rs @@ -1,3 +1,4 @@ +use geo::intersects; use numpy::ndarray::*; use crate::utils::*; use crate::interpolate::*; @@ -184,7 +185,8 @@ impl RectGrid { // on which side or corner is intersected, we find the next cell and repeat the process. // Since the line also has an intersection where it 'entered' the new cell, // we ignore this intersection and look for anohter. - // The loop is terminated when the line does not intersect any new corners or lines. + // The loop is terminated when the line does not intersect any new corners or lines, + // or when corner or cell that intersects the line also intersects point2, the end of the line. // Since the corners are also part of the lines as far as the line.intersects function is concerned, // we check the corners first and skip the lines check if a corner is intersected. // We then also have to ignore the lines that contain the corner when checking intersections for the next cell. @@ -202,36 +204,129 @@ impl RectGrid { // let mut ids = Array2::::zeros((0, 2)); let mut cell_id = self.cell_at_point(&p1.into_shape((1, p1.len())).unwrap()); - let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)]].view()); // Create a LineString from the supplied endpoints let point1 = Coord:: {x:p1[Ix1(0)], y:p1[Ix1(1)]}; let point2 = Coord:: {x:p2[Ix1(0)], y:p2[Ix1(1)]}; let line = LineString::new(vec![point1, point2]); - // TODO: Check if the line starts on a cell corner. If so, find out which of the connecting cells is the true starting cell. - let mut side_intersection_id: i64 = -1; let mut corner_intersection_id: i64 = -1; let mut skip_corners = array![false, false, false, false]; let mut skip_sides = array![false, false, false, false]; + + // Check if the line starts on a cell corner. If so, + // find out which of the connecting cells is the true starting cell and + // add the starting corner to 'skip_corners' when determining the next cell. + let mut p1_on_corner: bool = false; + let corners = self.cell_corners(&cell_id.view()); + for i in 0..4 { // loop over corners + p1_on_corner = point1.intersects(&Coord:: {x:corners[Ix3(0,i,0)], y:corners[Ix3(0,i,1)]}); + if p1_on_corner { + break + } + } + if p1_on_corner { + let direction_x = p2[Ix1(0)] - p1[Ix1(0)]; + let direction_y = p2[Ix1(1)] - p1[Ix1(1)]; + let mut azimuth = direction_x.atan2(direction_y) * 180. / std::f64::consts::PI; + azimuth += 360. * (azimuth <= 0.) as i64 as f64; + azimuth -= 360. * (azimuth > 360.) as i64 as f64; + match azimuth { // Line continues in top-right direction + az if (az >= 0.0 && az <= 90.0) => { + // Already in correct starting cell + skip_corners = array![true, false, false, false]; + skip_sides = array![true, false, false, true]; + } // Line continues in bottom-right direction + az if (az > 90.0 && az <= 180.0) => { + cell_id[Ix2(0,1)] -= 1; + skip_corners = array![false, false, false, true]; + skip_sides = array![false, false, true, true]; + } // Line continues in bottom-left direction + az if (az > 180.0 && az <= 270.0) => { + cell_id[Ix2(0,0)] -= 1; + cell_id[Ix2(0,1)] -= 1; + skip_corners = array![false, false, true, false]; + skip_sides = array![false, true, true, false]; + } + _ => { // Line continues in top-left direction + cell_id[Ix2(0,0)] -= 1; + skip_corners = array![false, true, false, false]; + skip_sides = array![true, true, false, false]; + } + } + } + + // Also check if the line starts on a cell side. If so, + // find out which of the connecting cells is the true starting cell and + // add the starting side to 'skip_sides' when determining the next cell. + let mut p1_on_side: i8 = -1; + let sides = vec![ + LineString::new(vec![Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]},Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]},Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]},Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]},Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]}]), + ]; + for i in 0..4 { // loop over sides + let intersects = sides[i].intersects(&point1); + if intersects { + p1_on_side = i as i8; + break + } + } + match p1_on_side { + // Because of the way cell_at_point handles the boundary conditions, + // we only need to check for the point being on the left or bottom sides. + // Points on the right or top sides are considered to be in the next cell over. + 0 => { // Bottom-left corner + if p2[Ix1(1)] > p1[Ix1(1)] { // Line continues towards the top + // cell_id already correct + skip_sides = array![true, false, false, false]; + } else { // Line continues towards the bottom + cell_id[Ix2(0,1)] -= 1; + skip_sides = array![false, false, true, false]; + } + } + 3 => { // Top-left corner + if p2[Ix1(0)] > p1[Ix1(0)] { // Line continues towards the right + // cell_id already correct + skip_sides = array![false, false, false, true]; + } else { // Line continues towards the left + cell_id[Ix2(0,0)] -= 1; + skip_sides = array![false, true, false, false]; + } + } + _ => {} // No intersection, check sides next + } + + // Add starting cell to return ids + let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)]].view()); + + // Loop and continuousely find the next cell until the end of the line is reached. loop { let corners = self.cell_corners(&cell_id.view()); - + let mut reached_point2 = false; corner_intersection_id = -1; for i in 0..4 { // Loop over corners if skip_corners[i] { // Discount the intersection towards previous cell - println!("Skipping corner {}", i); continue; } - let intersects = line.intersects(&Coord:: {x:corners[Ix3(0,i,0)], y:corners[Ix3(0,i,1)]}); + let corner_coord = Coord:: {x:corners[Ix3(0,i,0)], y:corners[Ix3(0,i,1)]}; + let intersects = line.intersects(&corner_coord); if intersects { + if corner_coord.intersects(&point2) { + reached_point2 = true; + } corner_intersection_id = i as i64; - side_intersection_id = -1; // Reset break; } } + if reached_point2 { + // Line ends in this cell, break infinite loop and return + break + } + match corner_intersection_id { // Adjust the cell-id to reflect the next cell and mark the oppisite corner and connecting sides to be skipped. // To demonstrate what I mean with skipping the opposite corner and connecting sides: @@ -267,7 +362,6 @@ impl RectGrid { // Add previous cell to vec and don't bother checking side intersections if we have a corner intersection if corner_intersection_id != -1 { - println!("Crossed corner {}", corner_intersection_id); let _ = ids.push(Axis(0), cell_id.slice(s![0, ..]).view()); continue } @@ -286,17 +380,23 @@ impl RectGrid { side_intersection_id = -1; for i in 0..sides.len() { if skip_sides[i] { // Discount the intersection towards previous cell - println!("Skipping side {}", i); continue; } let intersects = sides[i].intersects(&line); if intersects { + if sides[i].intersects(&point2) { + reached_point2 = true; + } side_intersection_id = i as i64; - corner_intersection_id = -1; // Also reset the corner break; } } + if reached_point2 { + // Line ends in this cell, break infinite loop and return + break + } + match side_intersection_id { // Adjust the cell-id to reflect the next cell and mark the oppisite side to be skipped. // To demonstrate what I mean with skipping the opposite side: @@ -323,7 +423,6 @@ impl RectGrid { break; } } - println!("Crossed side {}", side_intersection_id); let _ = ids.push(Axis(0), cell_id.slice(s![0, ..]).view()); } return ids From 293eda330932a9137813889c54dc246320f3d94b Mon Sep 17 00:00:00 2001 From: Timo Millenaar Date: Sat, 30 Mar 2024 19:22:37 +0100 Subject: [PATCH 4/5] Only check line starting on side if not already starting on corner --- src/rect_grid.rs | 72 ++++++++++++++++++++++++------------------------ 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/src/rect_grid.rs b/src/rect_grid.rs index c130f38b1..4fe7a178c 100644 --- a/src/rect_grid.rs +++ b/src/rect_grid.rs @@ -255,49 +255,49 @@ impl RectGrid { skip_sides = array![true, true, false, false]; } } - } - + } else { // Starting point not on corner // Also check if the line starts on a cell side. If so, // find out which of the connecting cells is the true starting cell and // add the starting side to 'skip_sides' when determining the next cell. - let mut p1_on_side: i8 = -1; - let sides = vec![ - LineString::new(vec![Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]},Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]}]), - LineString::new(vec![Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]},Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]}]), - LineString::new(vec![Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]},Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]}]), - LineString::new(vec![Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]},Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]}]), - ]; - for i in 0..4 { // loop over sides - let intersects = sides[i].intersects(&point1); - if intersects { - p1_on_side = i as i8; - break - } - } - match p1_on_side { - // Because of the way cell_at_point handles the boundary conditions, - // we only need to check for the point being on the left or bottom sides. - // Points on the right or top sides are considered to be in the next cell over. - 0 => { // Bottom-left corner - if p2[Ix1(1)] > p1[Ix1(1)] { // Line continues towards the top - // cell_id already correct - skip_sides = array![true, false, false, false]; - } else { // Line continues towards the bottom - cell_id[Ix2(0,1)] -= 1; - skip_sides = array![false, false, true, false]; + let mut p1_on_side: i8 = -1; + let sides = vec![ + LineString::new(vec![Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]},Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,1,0)], y:corners[Ix3(0,1,1)]},Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,2,0)], y:corners[Ix3(0,2,1)]},Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]}]), + LineString::new(vec![Coord:: {x:corners[Ix3(0,3,0)], y:corners[Ix3(0,3,1)]},Coord:: {x:corners[Ix3(0,0,0)], y:corners[Ix3(0,0,1)]}]), + ]; + for i in 0..4 { // loop over sides + let intersects = sides[i].intersects(&point1); + if intersects { + p1_on_side = i as i8; + break } } - 3 => { // Top-left corner - if p2[Ix1(0)] > p1[Ix1(0)] { // Line continues towards the right - // cell_id already correct - skip_sides = array![false, false, false, true]; - } else { // Line continues towards the left - cell_id[Ix2(0,0)] -= 1; - skip_sides = array![false, true, false, false]; + match p1_on_side { + // Because of the way cell_at_point handles the boundary conditions, + // we only need to check for the point being on the left or bottom sides. + // Points on the right or top sides are considered to be in the next cell over. + 0 => { // Bottom-left corner + if p2[Ix1(1)] > p1[Ix1(1)] { // Line continues towards the top + // cell_id already correct + skip_sides = array![true, false, false, false]; + } else { // Line continues towards the bottom + cell_id[Ix2(0,1)] -= 1; + skip_sides = array![false, false, true, false]; + } + } + 3 => { // Top-left corner + if p2[Ix1(0)] > p1[Ix1(0)] { // Line continues towards the right + // cell_id already correct + skip_sides = array![false, false, false, true]; + } else { // Line continues towards the left + cell_id[Ix2(0,0)] -= 1; + skip_sides = array![false, true, false, false]; + } } + _ => {} // No intersection, check sides next } - _ => {} // No intersection, check sides next - } + } // End checking for starting point on corner or side // Add starting cell to return ids let _ = ids.push(Axis(0), array![cell_id[Ix2(0,0)], cell_id[Ix2(0,1)]].view()); From 39c955226a6e01ae6b5ad24bdd65677e4e522d07 Mon Sep 17 00:00:00 2001 From: Timo Millenaar Date: Sat, 30 Mar 2024 19:41:01 +0100 Subject: [PATCH 5/5] Expose cells_intersecting_line in BaseGrid and RectGrid and add docstring and tests --- gridkit/base_grid.py | 26 ++++++ gridkit/rect_grid.py | 7 ++ tests/test_gridkit/test_rect_grid.py | 117 +++++++++++++++++++++++++++ 3 files changed, 150 insertions(+) diff --git a/gridkit/base_grid.py b/gridkit/base_grid.py index f30f06bd7..161837923 100644 --- a/gridkit/base_grid.py +++ b/gridkit/base_grid.py @@ -503,6 +503,32 @@ def _geom_iterator(): intersecting_cells.extend(cells_in_bounds[mask]) return GridIndex(intersecting_cells).unique() + @abc.abstractmethod + def cells_intersecting_line(self, line): + """Obtain the cell ids of the cells that are intersected the a supplied line. + + Similar functionality is found in :meth:`.BaseGrid.intersect_geometries`. + That method is more general and also handles point and polygon geometries. + However, :meth:`cells_intersecting_line` is more performant than :meth:`.BaseGrid.intersect_geometries` + when dealing with only line geometries. + + Parameters + ---------- + line: `numpy.ndarray` + A line as described by two points in the form [[x1,y1], [x2,y2]] + + Returns + ------- + :class:`GridIndex` + The ids of the cells intersected by the supplied line + + See also + -------- + :meth:`.BaseGrid.intersect_geometries` + The intersect_geometries method is more general, but less performant. + """ + pass + @validate_index def to_shapely(self, index, as_multipolygon: bool = False): """Represent the cells as Shapely Polygons diff --git a/gridkit/rect_grid.py b/gridkit/rect_grid.py index 8988b1916..0815f9731 100644 --- a/gridkit/rect_grid.py +++ b/gridkit/rect_grid.py @@ -531,6 +531,13 @@ def cells_in_bounds(self, bounds, return_cell_count: bool = False): ids = GridIndex(ids.T.reshape((*shape, 2))) return (ids, shape) if return_cell_count else ids + def cells_intersecting_line(self, line): + line = numpy.array(line, dtype=float) + if not line.shape == (2,2): + raise ValueError(f"Expected a line to be supplied in the form [[x1,y1], [x2,y2]]. Got shape {line.shape}") + cell_ids = self._grid.cells_intersecting_line(*line) + return GridIndex(cell_ids) + @property def parent_grid_class(self): return RectGrid diff --git a/tests/test_gridkit/test_rect_grid.py b/tests/test_gridkit/test_rect_grid.py index 75c359e45..d9f5fdbbb 100644 --- a/tests/test_gridkit/test_rect_grid.py +++ b/tests/test_gridkit/test_rect_grid.py @@ -459,3 +459,120 @@ def test_dx_dy_setter(): grid.dy = 0 with pytest.raises(ValueError): grid.dy = -1 + + +@pytest.mark.parametrize( + "line, expected_ids", + ( + ( + numpy.array( + [ + [-2.5, -1], + [2.5, 3], + ] + ), + numpy.array( + [[-3, -1], [-2, -1], [-2, 0], [-1, 0], [0, 1], [1, 1], [1, 2], [2, 2]] + ), + ), + ( + numpy.array( + [ + [-2.5, 1], + [2.5, -3], + ] + ), + numpy.array( + [ + [-3, 0], + [-2, 0], + [-2, -1], + [-1, -1], + [0, -2], + [1, -2], + [1, -3], + [2, -3], + ] + ), + ), + ( + numpy.array( + [ + [3, 2.5], + [-1, -2.5], + ] + ), + numpy.array( + [[2, 2], [2, 1], [1, 1], [1, 0], [0, -1], [0, -2], [-1, -2], [-1, -3]] + ), + ), + ( + numpy.array( + [ + [-2, -2.5], + [2, 2.5], + ] + ), + numpy.array( + [[-2, -3], [-2, -2], [-1, -2], [-1, -1], [0, 0], [0, 1], [1, 1], [1, 2]] + ), + ), + ( + numpy.array( + [ + [-2, -2], + [2, 2], + ] + ), + numpy.array([[-2, -2], [-1, -1], [0, 0], [1, 1]]), + ), + ( + numpy.array( + [ + [-3, 3], + [3, -3], + ] + ), + numpy.array([[-3, 2], [-2, 1], [-1, 0], [0, -1], [1, -2], [2, -3]]), + ), + ( + numpy.array( + [ + [-3, -2], + [3, 2], + ] + ), + numpy.array( + [[-3, -2], [-2, -2], [-2, -1], [-1, -1], [0, 0], [1, 0], [1, 1], [2, 1]] + ), + ), + ( + numpy.array( + [ + [-3, -1], + [3, 1], + ] + ), + numpy.array([[-3, -1], [-2, -1], [-1, -1], [0, 0], [1, 0], [2, 0]]), + ), + ( + numpy.array( + [ + [2, 1], + [-2, -1], + ] + ), + numpy.array([[1, 0], [0, 0], [-1, -1], [-2, -1]]), + ), + ), +) +def test_cells_intersecting_line(line, expected_ids): + grid = RectGrid(dx=1, dy=1) + + line = numpy.array(line) + ids = grid.cells_intersecting_line(line) + numpy.testing.assert_allclose(ids, expected_ids) + + # Make sure the returned indices don't depend on the order of the two points + ids_rev = grid.cells_intersecting_line(line[::-1]) + numpy.testing.assert_allclose(ids, ids_rev[::-1])