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/gridkit/base_grid.py b/gridkit/base_grid.py index 9fe6e0595..38b0b9a92 100644 --- a/gridkit/base_grid.py +++ b/gridkit/base_grid.py @@ -535,6 +535,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 136c9fac0..b246f3cf7 100644 --- a/gridkit/rect_grid.py +++ b/gridkit/rect_grid.py @@ -620,6 +620,15 @@ 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/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 6796b0305..5ba3fb5a9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -261,6 +261,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..4fe7a178c 100644 --- a/src/rect_grid.rs +++ b/src/rect_grid.rs @@ -1,5 +1,10 @@ +use geo::intersects; 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 +177,255 @@ impl RectGrid { } nearby_cells } + + 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, + // 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. + // + // 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()); + + // 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 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]; + } + } + } 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]; + } + } + 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 + } + } // 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()); + + // 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 + continue; + } + 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; + 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: + // 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 { + 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)]}]), + ]; + + side_intersection_id = -1; + for i in 0..sides.len() { + if skip_sides[i] { // Discount the intersection towards previous cell + continue; + } + let intersects = sides[i].intersects(&line); + if intersects { + if sides[i].intersects(&point2) { + reached_point2 = true; + } + side_intersection_id = i as i64; + 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: + // 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; + skip_sides = array![false, false, true, false]; + } + 1 => { // Right side + cell_id[Ix2(0,0)] += 1; + skip_sides = array![false, false, false, true]; + } + 2 => { // Top side + cell_id[Ix2(0,1)] += 1; + skip_sides = array![true, false, false, false]; + } + 3 => { // Left side + cell_id[Ix2(0,0)] -= 1; + skip_sides = array![false, true, false, false]; + } + _ => { // No intersection + // Reached the end point of the line, break infinite loop and return from function + break; + } + } + let _ = ids.push(Axis(0), cell_id.slice(s![0, ..]).view()); + } + return ids + } + } diff --git a/tests/test_gridkit/test_rect_grid.py b/tests/test_gridkit/test_rect_grid.py index ef5d438bf..72123c8ae 100644 --- a/tests/test_gridkit/test_rect_grid.py +++ b/tests/test_gridkit/test_rect_grid.py @@ -506,6 +506,123 @@ def test_dx_dy_setter(): 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]) + + def test_size_setter(): grid = RectGrid(dx=1.23, dy=4.56, rotation=10) assert grid.size is None