diff --git a/README.md b/README.md index 23c261a..5904183 100644 --- a/README.md +++ b/README.md @@ -31,7 +31,7 @@ The primary use case for this project is to find the longest line of sight on th Any file format supported by `gdal` should be supported. It must follow this requirements: * It must be perfectly square. * The width must be divisible by 48. -* It should be in an AEQD or similar metric projection whose anchor is the centre of the tile. +* It must be in an AEQD or similar metric projection whose anchor is the centre of the tile. * All points must be the same metric distance apart. The best source of elevation data I have found is here: https://www.viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm It's mostly in the `.hgt` format, but can easily be converted to `.bt` using `gdal_translate`. @@ -57,9 +57,10 @@ Once that's run, you will have access to what's called "ring data", stored by de Usage: tvs Commands: - compute Run main computations - viewshed Reconstruct a viewshed - help Print this message or the help of the given subcommand(s) + compute Run main computations + viewshed Reconstruct a viewshed + post-process Post-processing tasks + help Print this message or the help of the given subcommand(s) Options: -h, --help @@ -98,7 +99,7 @@ Once that's run, you will have access to what's called "ring data", stored by de - cpu: Optimised cache-efficient CPU kernel - cuda: TBC - [default: vulkan] + [default: cpu] --output-dir Directory to save results in @@ -134,11 +135,25 @@ Once that's run, you will have access to what's called "ring data", stored by de [default: 0.13] - --thread-count + --thread-count Thread count used for CPU parallelism [default: 8] + --disable-image-render + Controls line of sight and total viewshed image generation + + --viewsheds-db-path + Where to store the viewshed data. Requires build with `--features=ring_data` + + [default: ./viewsheds.db] + + --aoi-point + Lon/lat coordinates for a polygon that represents an Area of Interest within the DEM. Points outside this polygon will be ignored + + --database-per-thread + Creates a database per thread, managed by an additional worker thread + -h, --help Print help (see a summary with '-h') @@ -146,11 +161,14 @@ Once that's run, you will have access to what's called "ring data", stored by de * Viewshed Reconstruct a viewshed - Usage: tvs viewshed [COORDINATES]... + Usage: tvs viewshed [COORDINATES]... Arguments: - - Directory where compute output was saved + + Path or directory where viewshed DB was saved + + + Directory to save viewshed in [COORDINATES]... Coordinates to reconstruct viewsheds for @@ -158,6 +176,24 @@ Once that's run, you will have access to what's called "ring data", stored by de Options: -h, --help Print help + + +* Post-process + Post-processing tasks + + Usage: tvs post-process [OPTIONS] + + Arguments: + + Directory where viewshed DBs were saved + + Options: + --threads + Number of threads to use. Defaults to letting Rayon choose the thread count + + -h, --help + Print help + ``` ## Building Vulkan shader diff --git a/benchmarks/run.sh b/benchmarks/run.sh index e08b926..eba1776 100755 --- a/benchmarks/run.sh +++ b/benchmarks/run.sh @@ -32,8 +32,6 @@ rm $sqlite_db_path || true # Do the calculations time cargo run --features ring_data --release -- \ compute "$PROJECT_ROOT/benchmarks/cardiff.tiff" \ - --scale 100 \ - --centre-from-projection \ --rings-per-km 3 \ --backend "$backend" \ --process all \ diff --git a/benchmarks/samples/aeqd_bad_centre.tiff b/benchmarks/samples/aeqd_bad_centre.tiff new file mode 100644 index 0000000..f9e2611 Binary files /dev/null and b/benchmarks/samples/aeqd_bad_centre.tiff differ diff --git a/crates/total-viewsheds/src/config.rs b/crates/total-viewsheds/src/config.rs index e2c66e5..699e838 100644 --- a/crates/total-viewsheds/src/config.rs +++ b/crates/total-viewsheds/src/config.rs @@ -72,10 +72,6 @@ pub struct Compute { )] pub output_dir: std::path::PathBuf, - /// Override the calculated DEM points scale from the DEM file. Units in meters. - #[arg(long, value_name = "DEM scale (meters)")] - pub scale: Option, - /// What to compute. #[arg( long, @@ -116,15 +112,6 @@ pub struct Compute { #[arg(long, value_name = "Render image", default_value = "false")] pub disable_image_render: bool, - /// Derive the tile's centre from the tile's anchored projection. This can be more accurate for - /// large metric-projected tiles. - #[arg( - long, - value_name = "Get centre from projection", - default_value = "false" - )] - pub centre_from_projection: bool, - /// Where to store the viewshed data. Requires build with `--features=ring_data`. #[arg( long, @@ -142,6 +129,7 @@ pub struct Compute { value_name = "Lon/lat coords of region") ] pub aoi_point: Vec<(f32, f32)>, + /// Creates a database per thread, managed by an additional worker thread #[arg( long, diff --git a/crates/total-viewsheds/src/cpu/unrolled_los.rs b/crates/total-viewsheds/src/cpu/unrolled_los.rs index 7ec5b4d..48f6bcf 100644 --- a/crates/total-viewsheds/src/cpu/unrolled_los.rs +++ b/crates/total-viewsheds/src/cpu/unrolled_los.rs @@ -87,7 +87,8 @@ impl UnrolledVectorLos Result { let size: u64 = u64::from(width) * u64::from(width); - #[expect( - clippy::cast_possible_truncation, - clippy::as_conversions, - clippy::cast_sign_loss, - clippy::cast_precision_loss, - reason = "This shouldn't be a problem in most sane cases" - )] - let max_los_as_points = (max_line_of_sight as f32 / scale) as u32; let max_possible_los_as_points = width.div_euclid(2); - if max_los_as_points > max_possible_los_as_points { + if max_line_of_sight_as_points > max_possible_los_as_points { color_eyre::eyre::bail!( - "The maximum line of sight ({max_line_of_sight}m) is longer than the maximum \ + "The maximum line of sight ({max_line_of_sight_as_points}) is longer than the maximum \ distance that any point can completely see ({}m).", f64::from(max_possible_los_as_points) * f64::from(scale) ); @@ -78,23 +70,12 @@ impl DEM { size, scale, centre: centre_latlon, - max_los_as_points, + max_los_as_points: max_line_of_sight_as_points, computable_points_count, }; Ok(dem) } - #[expect( - clippy::as_conversions, - clippy::cast_possible_truncation, - clippy::cast_sign_loss, - reason = "The scale is positive and never be beyond the limts of f32" - )] - /// Cast the scale into an integer - pub const fn scale_u32(&self) -> u32 { - self.scale as u32 - } - /// Is the DEM ID within the DEM such that a valid viewshed can be made for it? pub fn is_point_computable(&self, dem_id: u32) -> bool { let coord = self.convert_dem_id_to_coord(dem_id).0; diff --git a/crates/total-viewsheds/src/main.rs b/crates/total-viewsheds/src/main.rs index 86aedce..c743cb0 100644 --- a/crates/total-viewsheds/src/main.rs +++ b/crates/total-viewsheds/src/main.rs @@ -140,7 +140,19 @@ fn compute(config: &config::Compute) -> Result<()> { } let mut tile = tile::Tile::load(config)?; - let scale = config.scale.unwrap_or(tile.scale); + + let max_line_of_sight_as_points = tile.width.div_euclid(3); + + let mut dem = crate::dem::DEM::new( + tile.centre, + tile.width, + tile.scale, + max_line_of_sight_as_points, + )?; + + dem.elevations = mem::take(&mut tile.data); + + tracing::debug!("Created DEM: {dem:?}"); #[expect( clippy::as_conversions, @@ -149,30 +161,18 @@ fn compute(config: &config::Compute) -> Result<()> { clippy::cast_precision_loss, reason = "Sign loss and truncation aren't relevant" )] - let max_line_of_sight_metres = (tile.width.div_euclid(3) as f32 * scale) as u32; - - let mut dem = crate::dem::DEM::new(tile.centre, tile.width, scale, max_line_of_sight_metres)?; - - dem.elevations = mem::take(&mut tile.data); + let max_line_of_sight_as_metres = (max_line_of_sight_as_points as f32 * tile.scale) as u32; // Free up RAM drop(tile); - tracing::debug!("Created DEM: {dem:?}"); - - let max_line_of_sight = if matches!(config.backend, config::Backend::CPU) { - dem.max_los_as_points - } else { - max_line_of_sight_metres - }; - let dem_metadata = crate::cpu::storage::metadata::MetaData { width: dem.width, scale: dem.scale, - max_line_of_sight, + max_line_of_sight: max_line_of_sight_as_points, reserved_ring_size: run::compute::Compute::ring_count_per_band( config.rings_per_km, - dem.max_los_as_points * dem.scale_u32(), + max_line_of_sight_as_metres, ), centre: dem.centre, }; diff --git a/crates/total-viewsheds/src/run/compute.rs b/crates/total-viewsheds/src/run/compute.rs index b4157fe..acd167f 100644 --- a/crates/total-viewsheds/src/run/compute.rs +++ b/crates/total-viewsheds/src/run/compute.rs @@ -194,7 +194,7 @@ impl<'compute> Compute<'compute> { Ok(crate::cpu::storage::metadata::MetaData { width: self.dem.width, scale: self.dem.scale, - max_line_of_sight: self.dem.max_los_as_points * self.dem.scale_u32(), + max_line_of_sight: self.dem.max_los_as_points, reserved_ring_size: usize::try_from(self.constants.reserved_rings_per_band)?, centre: self.dem.centre, }) diff --git a/crates/total-viewsheds/src/tile.rs b/crates/total-viewsheds/src/tile.rs index 34d912d..964214b 100644 --- a/crates/total-viewsheds/src/tile.rs +++ b/crates/total-viewsheds/src/tile.rs @@ -6,8 +6,6 @@ use color_eyre::eyre::{ContextCompat as _, Result}; -use crate::projection; - /// The basic raw data needed to compute total viewsheds. pub struct Tile { /// The width of the tile. @@ -43,11 +41,7 @@ impl Tile { ); } - let centre = if config.centre_from_projection { - Self::get_centre_by_projection(&dataset)? - } else { - Self::get_centre_by_raster(&dataset)? - }; + let centre = Self::get_centre_by_projection(&dataset)?; let lon = centre.0.x; let lat = centre.0.y; @@ -103,10 +97,8 @@ impl Tile { } } - /// Get the lat/lon of the centre of the tile simply by dividing the extent in half. This isn't - /// ideal as there's no guarantee that it's the same centre that the creator used to generate - /// the tile. - fn get_centre_by_raster(dataset: &gdal::Dataset) -> Result { + /// Get the metric centre of the tile simply by dividing the extent in half. + fn get_metric_centre(dataset: &gdal::Dataset) -> Result { let (width, height) = dataset.raster_size(); #[expect( @@ -126,24 +118,27 @@ impl Tile { x_top_left + (centre_x * pixel_width), y_top_left - (centre_y * pixel_height), ); - let crs = &dataset.spatial_ref()?.to_proj4()?; - - let mut converted = (x_world, y_world, 0.0f64); - proj4rs::transform::transform( - &proj4rs::Proj::from_proj_string(crs)?, - &projection::Converter::degrees_projection()?, - &mut converted, - )?; - - Ok(crate::projection::LonLatCoord(geo::coord! { - x: converted.0.to_degrees(), - y: converted.1.to_degrees() - })) + + Ok(geo::coord! { + x: x_world, + y: y_world + }) } - /// Get the lat/lon centre of the tile by querying the projection's definition. This is more - /// likely to guarantee that the tile's centre matches the creator's intended centre. + /// Get the lat/lon centre of the tile by querying the projection's definition. fn get_centre_by_projection(dataset: &gdal::Dataset) -> Result { + const ALLOWED_DEVIATION: f64 = 0.001; // 1 milimetre + let geometric_centre = Self::get_metric_centre(dataset)?; + if geometric_centre.x.abs() > ALLOWED_DEVIATION + || geometric_centre.y.abs() > ALLOWED_DEVIATION + { + color_eyre::eyre::bail!( + "Tile centre ({},{}) must be within {ALLOWED_DEVIATION} of 0,0.", + geometric_centre.x, + geometric_centre.y + ); + } + let projection = &dataset.spatial_ref()?.to_proj4()?; let lat_0: f64 = projection @@ -174,7 +169,6 @@ mod test { fn get_centre_by_projection() { let config = crate::config::Compute { input: "../../benchmarks/samples/aeqd_10x10.tiff".into(), - centre_from_projection: true, ..Default::default() }; let tile = Tile::load(&config).unwrap(); @@ -188,18 +182,18 @@ mod test { } #[test] - fn get_centre_by_raster() { - let config = crate::config::Compute { - input: "../../benchmarks/samples/aeqd_10x10.tiff".into(), - ..Default::default() - }; - let tile = Tile::load(&config).unwrap(); - assert_eq!( - tile.centre, - crate::projection::LonLatCoord(geo::Coord { - x: -0.1278, - y: 51.5074 - }) - ); + fn get_metric_centre() { + let dataset = gdal::Dataset::open("../../benchmarks/samples/aeqd_10x10.tiff").unwrap(); + let centre = Tile::get_metric_centre(&dataset).unwrap(); + assert_eq!(centre, geo::Coord { x: 0.0, y: 0.0 }); + } + + #[test] + fn bad_centre() { + let dataset = gdal::Dataset::open("../../benchmarks/samples/aeqd_bad_centre.tiff").unwrap(); + let error = Tile::get_centre_by_projection(&dataset) + .unwrap_err() + .to_string(); + assert!(error.contains("Tile centre"), "Wrong error: '{error}'"); } }