diff --git a/backend/Cargo.lock b/backend/Cargo.lock index c5f152a..60e53ca 100644 --- a/backend/Cargo.lock +++ b/backend/Cargo.lock @@ -73,6 +73,7 @@ dependencies = [ "bincode", "console_error_panic_hook", "console_log", + "contour", "geo", "geojson", "log", @@ -146,6 +147,18 @@ dependencies = [ "web-sys", ] +[[package]] +name = "contour" +version = "0.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6bcf05496233f4d666d941c9536e0dec58f9f04d7c280874dfecb4df88a72d9b" +dependencies = [ + "geo-types", + "lazy_static", + "rustc-hash", + "slab", +] + [[package]] name = "crc32fast" version = "1.3.2" @@ -677,6 +690,12 @@ dependencies = [ "smallvec", ] +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + [[package]] name = "rustc_version" version = "0.4.0" @@ -759,6 +778,15 @@ dependencies = [ "serde", ] +[[package]] +name = "slab" +version = "0.4.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f92a496fb766b417c996b9c5e57daf2f7ad3b0bebe1ccfca4856390e3d3bb67" +dependencies = [ + "autocfg", +] + [[package]] name = "smallvec" version = "1.11.2" diff --git a/backend/Cargo.toml b/backend/Cargo.toml index f11b6e7..5177588 100644 --- a/backend/Cargo.toml +++ b/backend/Cargo.toml @@ -23,3 +23,4 @@ wasm-bindgen = "0.2.87" web-sys = { version = "0.3.64", features = ["console"] } bincode = "1.3.3" petgraph = "0.6.4" +contour = "0.12.0" diff --git a/backend/src/lib.rs b/backend/src/lib.rs index 7fdfc66..668b469 100644 --- a/backend/src/lib.rs +++ b/backend/src/lib.rs @@ -12,11 +12,13 @@ use wasm_bindgen::prelude::*; use self::cells::Cell; use self::neighbourhood::Neighbourhood; +use self::render_cells::RenderCells; use self::shortcuts::Shortcuts; mod cells; mod mercator; mod neighbourhood; +mod render_cells; mod scrape; mod shortcuts; mod tags; diff --git a/backend/src/neighbourhood.rs b/backend/src/neighbourhood.rs index 59fc94c..5f0ce52 100644 --- a/backend/src/neighbourhood.rs +++ b/backend/src/neighbourhood.rs @@ -6,12 +6,14 @@ use geo::{ }; use geojson::{Feature, GeoJson, Geometry}; -use crate::{Cell, IntersectionID, MapModel, RoadID, Shortcuts}; +use crate::render_cells::Color; +use crate::{Cell, IntersectionID, MapModel, RenderCells, RoadID, Shortcuts}; pub struct Neighbourhood { pub interior_roads: HashSet, crosses: HashMap, pub border_intersections: HashSet, + pub boundary_polygon: Polygon, } impl Neighbourhood { @@ -60,6 +62,7 @@ impl Neighbourhood { interior_roads, crosses, border_intersections, + boundary_polygon: boundary, }) } @@ -69,16 +72,9 @@ impl Neighbourhood { // TODO Decide how/where state lives let modal_filters = BTreeMap::new(); let cells = Cell::find_all(map, self, &modal_filters); + let render_cells = RenderCells::new(map, self, &cells); let shortcuts = Shortcuts::new(map, self); - // TODO Temporary rendering - let mut road_to_cell = HashMap::new(); - for (idx, cell) in cells.iter().enumerate() { - for (r, _) in &cell.roads { - road_to_cell.insert(*r, idx); - } - } - for r in &self.interior_roads { let mut f = map.get_r(*r).to_gj(&map.mercator); f.set_property("kind", "interior_road"); @@ -86,11 +82,6 @@ impl Neighbourhood { "shortcuts", shortcuts.count_per_road.get(r).cloned().unwrap_or(0), ); - if let Some(cell) = road_to_cell.get(r) { - f.set_property("cell", *cell); - } else { - error!("A road has no cell"); - } features.push(f); } for (r, pct) in &self.crosses { @@ -105,6 +96,19 @@ impl Neighbourhood { features.push(f); } + for (polygons, color) in render_cells + .polygons_per_cell + .into_iter() + .zip(render_cells.colors) + { + let mut f = Feature::from(Geometry::from(&map.mercator.to_wgs84(&polygons))); + match color { + Color::Disconnected => f.set_property("color", "disconnected"), + Color::Cell(idx) => f.set_property("color", idx), + } + features.push(f); + } + GeoJson::from(features) } } diff --git a/backend/src/render_cells.rs b/backend/src/render_cells.rs new file mode 100644 index 0000000..c0e4bb3 --- /dev/null +++ b/backend/src/render_cells.rs @@ -0,0 +1,323 @@ +use std::collections::{HashSet, VecDeque}; + +use geo::{BoundingRect, Densify, LineString, MultiPolygon, Polygon, Rect, Scale, Translate}; + +use crate::{Cell, MapModel, Neighbourhood}; + +const NUM_COLORS: usize = 10; +const RESOLUTION_M: f64 = 10.0; + +pub enum Color { + Disconnected, + Cell(usize), +} + +pub struct RenderCells { + /// Rarely, this might be empty if the area is very small + pub polygons_per_cell: Vec, + /// Colors per cell, such that adjacent cells are colored differently + pub colors: Vec, +} + +impl RenderCells { + /// Partition a neighbourhood's boundary polygon based on the cells. This discretizes space into + /// a grid, and then extracts a polygon from the raster. The results don't look perfect, but + /// it's fast. + pub fn new(map: &MapModel, neighbourhood: &Neighbourhood, cells: &Vec) -> RenderCells { + let boundary_polygon = neighbourhood.boundary_polygon.clone(); + // Make a 2D grid covering the polygon. Each tile in the grid contains a cell index, which + // will become a color by the end. None means no cell is assigned yet. + //let bounds: Rect = boundary_polygon.bounding_rect().into().unwrap(); + let bounds = >::from(boundary_polygon.bounding_rect()).unwrap(); + let mut grid: Grid> = Grid::new( + (bounds.width() / RESOLUTION_M).ceil() as usize, + (bounds.height() / RESOLUTION_M).ceil() as usize, + None, + ); + + // Initially fill out the grid based on the roads in each cell + let mut warn_leak = true; + for (cell_idx, cell) in cells.iter().enumerate() { + for (r, interval) in &cell.roads { + let road = map.get_r(*r); + // Some roads with a filter are _very_ short, and this fails. The connecting roads + // on either side should contribute a grid cell and wind up fine. + if let Some(slice) = + slice_linestring(&road.linestring, interval.start, interval.end) + { + // Walk along the center line. We could look at the road's thickness and fill + // out points based on that, but the diffusion should take care of it. + // TODO Clean up old comments like this + for pt in slice.densify(RESOLUTION_M / 2.0).0 { + let grid_idx = grid.idx( + ((pt.x - bounds.min().x) / RESOLUTION_M) as usize, + ((pt.y - bounds.min().y) / RESOLUTION_M) as usize, + ); + // Due to tunnels/bridges, sometimes a road belongs to a neighbourhood, but + // leaks outside the neighbourhood's boundary. Avoid crashing. The real fix + // is to better define boundaries in the face of z-order changes. + // + // Example is https://www.openstreetmap.org/way/87298633 + if grid_idx >= grid.data.len() { + if warn_leak { + warn!( + "{} leaks outside its neighbourhood's boundary polygon, near {:?}", + road.id, pt + ); + // In some neighbourhoods, there are so many warnings that logging + // causes noticeable slowdown! + warn_leak = false; + } + continue; + } + + // If roads from two different cells are close enough to clobber + // originally, oh well? + grid.data[grid_idx] = Some(cell_idx); + } + } + } + } + // Also mark the boundary polygon, so we can prevent the diffusion from "leaking" outside + // the area. The grid covers the rectangular bounds of the polygon. Rather than make an + // enum with 3 cases, just assign a new index to mean "boundary." + let boundary_marker = cells.len(); + for pt in boundary_polygon.exterior().densify(RESOLUTION_M / 2.0).0 { + // TODO Refactor helpers to transform between map-space and the grid tiles. Possibly + // Grid should know about this. + let grid_idx = grid.idx( + ((pt.x - bounds.min().x) / RESOLUTION_M) as usize, + ((pt.y - bounds.min().y) / RESOLUTION_M) as usize, + ); + grid.data[grid_idx] = Some(boundary_marker); + } + + let adjacencies = diffusion(&mut grid, boundary_marker); + let mut cell_colors = color_cells(cells.len(), adjacencies); + + // Color some special cells + for (idx, cell) in cells.iter().enumerate() { + if cell.is_disconnected() { + cell_colors[idx] = Color::Disconnected; + } + } + + finalize(grid, cell_colors, bounds, boundary_polygon) + } +} + +fn finalize( + grid: Grid>, + cell_colors: Vec, + bounds: Rect, + boundary_polygon: Polygon, +) -> RenderCells { + let mut result = RenderCells { + polygons_per_cell: Vec::new(), + colors: Vec::new(), + }; + + for (idx, color) in cell_colors.into_iter().enumerate() { + // contour will find where the grid is >= a threshold value. The main grid has one + // number per cell, so we can't directly use it -- the area >= some cell index is + // meaningless. Per cell, make a new grid that just has that cell. + let grid: Grid = Grid { + width: grid.width, + height: grid.height, + data: grid + .data + .iter() + .map( + |maybe_cell| { + if maybe_cell == &Some(idx) { + 1.0 + } else { + 0.0 + } + }, + ) + .collect(), + }; + + let smooth = false; + let contour_builder = + contour::ContourBuilder::new(grid.width as u32, grid.height as u32, smooth); + let thresholds = vec![1.0]; + + let mut cell_polygons = Vec::new(); + for contour in contour_builder.contours(&grid.data, &thresholds).unwrap() { + // TODO Check this API again + let (polygons, _) = contour.into_inner(); + for p in polygons { + if let Ok(poly) = Polygon::try_from(p) { + cell_polygons.push( + poly.scale(RESOLUTION_M) + .translate(bounds.min().x, bounds.min().y), + ); + } + } + } + + // Sometimes one cell "leaks" out of the neighbourhood boundary. Not sure why. But we + // can just clip the result. + let mut clipped = Vec::new(); + for p in cell_polygons { + // If clipping fails, just use the original polygon. + /*if let Ok(list) = p.intersection(&result.boundary_polygon) { + clipped.extend(list); + } else {*/ + clipped.push(p); + //} + } + + result.polygons_per_cell.push(MultiPolygon::new(clipped)); + result.colors.push(color); + } + + result +} + +/// Returns a set of adjacent indices. The pairs are symmetric -- (x, y) and (y, x) will both be +/// populated. Adjacency with boundary_marker doesn't count. +fn diffusion(grid: &mut Grid>, boundary_marker: usize) -> HashSet<(usize, usize)> { + // Grid indices to propagate + let mut queue: VecDeque = VecDeque::new(); + + // Initially seed the queue with all colored tiles + for (idx, value) in grid.data.iter().enumerate() { + if let Some(x) = value { + // Don't expand the boundary tiles + if *x != boundary_marker { + queue.push_back(idx); + } + } + } + + let mut adjacencies = HashSet::new(); + + while !queue.is_empty() { + let current_idx = queue.pop_front().unwrap(); + let current_color = grid.data[current_idx].unwrap(); + let (current_x, current_y) = grid.xy(current_idx); + // Don't flood to diagonal neighbors. That would usually result in "leaking" out past the + // boundary tiles when the boundary polygon isn't axis-aligned. + // TODO But this still does "leak" out sometimes -- the cell covering 22nd/Lynn, for + // example. + for (next_x, next_y) in grid.orthogonal_neighbors(current_x, current_y) { + let next_idx = grid.idx(next_x, next_y); + if let Some(prev_color) = grid.data[next_idx] { + // If the color doesn't match our current_color, we've found the border between two + // cells. + if current_color != prev_color + && current_color != boundary_marker + && prev_color != boundary_marker + { + adjacencies.insert((current_color, prev_color)); + adjacencies.insert((prev_color, current_color)); + } + // If a color has been assigned, don't flood any further. + } else { + grid.data[next_idx] = Some(current_color); + queue.push_back(next_idx); + } + } + } + + adjacencies +} + +fn color_cells(num_cells: usize, adjacencies: HashSet<(usize, usize)>) -> Vec { + // This is the same greedy logic as Perimeter::calculate_coloring + let mut assigned_colors = Vec::new(); + for this_idx in 0..num_cells { + let mut available_colors: Vec = std::iter::repeat(true).take(NUM_COLORS).collect(); + // Find all neighbors + for other_idx in 0..num_cells { + if adjacencies.contains(&(this_idx, other_idx)) { + // We assign colors in order, so any neighbor index smaller than us has been + // chosen + if other_idx < this_idx { + available_colors[assigned_colors[other_idx]] = false; + } + } + } + + // If there are multiple colors available, prefer one that hasn't been used anywhere yet. + // Cells far apart shouldn't seem related to the user. + let mut choice = None; + let mut backup = None; + for (idx, available) in available_colors.into_iter().enumerate() { + if !available { + continue; + } + if assigned_colors.iter().any(|x| *x == idx) { + if backup.is_none() { + backup = Some(idx); + } + } else { + choice = Some(idx); + break; + } + } + assigned_colors.push( + choice + .or(backup) + .unwrap_or_else(|| assigned_colors.len() % NUM_COLORS), + ); + } + assigned_colors + .into_iter() + .map(|idx| Color::Cell(idx)) + .collect() +} + +/// A 2D grid containing some arbitrary data. +pub struct Grid { + /// Logically represents a 2D vector. Row-major ordering. + pub data: Vec, + pub width: usize, + pub height: usize, +} + +impl Grid { + pub fn new(width: usize, height: usize, default: T) -> Grid { + Grid { + data: std::iter::repeat(default).take(width * height).collect(), + width, + height, + } + } + + /// Calculate the index from a given (x, y). Doesn't do any bounds checking. + pub fn idx(&self, x: usize, y: usize) -> usize { + y * self.width + x + } + + /// The inverse of `idx`. No bounds checking. + pub fn xy(&self, idx: usize) -> (usize, usize) { + let y = idx / self.width; + let x = idx % self.width; + (x, y) + } + + /// From one tile, calculate the 4 orthogonal neighbors. Includes bounds checking. + pub fn orthogonal_neighbors(&self, center_x: usize, center_y: usize) -> Vec<(usize, usize)> { + let center_x = center_x as isize; + let center_y = center_y as isize; + let mut results = Vec::new(); + for (dx, dy) in [(-1, 0), (0, -1), (0, 1), (1, 0)] { + let x = center_x + dx; + let y = center_y + dy; + if x < 0 || (x as usize) >= self.width || y < 0 || (y as usize) >= self.height { + continue; + } + results.push((x as usize, y as usize)); + } + results + } +} + +// TODO Bring in that geo PR +fn slice_linestring(linestring: &LineString, start: f64, end: f64) -> Option { + Some(linestring.clone()) +} diff --git a/backend/src/shortcuts.rs b/backend/src/shortcuts.rs index 1c96b16..dcab7d2 100644 --- a/backend/src/shortcuts.rs +++ b/backend/src/shortcuts.rs @@ -2,7 +2,7 @@ use std::collections::HashMap; use petgraph::graphmap::DiGraphMap; -use crate::{IntersectionID, MapModel, Neighbourhood, RoadID}; +use crate::{MapModel, Neighbourhood, RoadID}; pub struct Shortcuts { pub count_per_road: HashMap, diff --git a/web/src/NeighbourhoodLayer.svelte b/web/src/NeighbourhoodLayer.svelte index dababa9..9fa8375 100644 --- a/web/src/NeighbourhoodLayer.svelte +++ b/web/src/NeighbourhoodLayer.svelte @@ -5,6 +5,7 @@ CircleLayer, FillLayer, GeoJSON, + hoverStateFilter, LineLayer, Popup, } from "svelte-maplibre"; @@ -12,26 +13,49 @@ constructMatchExpression, isLine, isPoint, + isPolygon, PropertiesTable, } from "./common"; export let model: MapModel; export let boundary: Feature; + // A qualitative palette from colorbrewer2.org, skipping the red hue (used + // for levels of shortcutting) and grey (too close to the basemap) + let cell_colors = [ + "#8dd3c7", + "#ffffb3", + "#bebada", + "#80b1d3", + "#fdb462", + "#b3de69", + "#fccde5", + "#bc80bd", + "#ccebc5", + "#ffed6f", + ]; + let details = JSON.parse(model.analyzeNeighbourhood(boundary)); let maxShortcuts = Math.max( ...details.features.map((f) => f.properties.shortcuts ?? 0) ); + for (let f of details.features) { + if (f.properties.color == "disconnected") { + f.properties.color = "red"; + } else if (f.properties.color) { + f.properties.color = cell_colors[f.properties.color % cell_colors.length]; + } + } - + +