From c5e5b4174cda747142d227369812fc799753f085 Mon Sep 17 00:00:00 2001 From: Thomas Buckley-Houston Date: Wed, 22 Apr 2026 15:57:51 +0100 Subject: [PATCH 1/4] fix: Disallow non-centred tiles There's probably a way for us to support non-centred tiles in the application, but for now it's best to just disallow them. Non-centred tiles compute fine but their final geo referencing causes misalignment if displayed on a global map. --- README.md | 2 +- benchmarks/run.sh | 1 - crates/total-viewsheds/src/config.rs | 9 ---- crates/total-viewsheds/src/tile.rs | 62 ++++++++++------------------ 4 files changed, 23 insertions(+), 51 deletions(-) diff --git a/README.md b/README.md index 23c261a..ffbc3d2 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`. diff --git a/benchmarks/run.sh b/benchmarks/run.sh index e08b926..80cffad 100755 --- a/benchmarks/run.sh +++ b/benchmarks/run.sh @@ -33,7 +33,6 @@ rm $sqlite_db_path || true 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/crates/total-viewsheds/src/config.rs b/crates/total-viewsheds/src/config.rs index e2c66e5..c530874 100644 --- a/crates/total-viewsheds/src/config.rs +++ b/crates/total-viewsheds/src/config.rs @@ -116,15 +116,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, diff --git a/crates/total-viewsheds/src/tile.rs b/crates/total-viewsheds/src/tile.rs index 34d912d..0da2886 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_centre_by_raster(dataset: &gdal::Dataset) -> Result { let (width, height) = dataset.raster_size(); #[expect( @@ -126,24 +118,24 @@ 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 { + let geometric_centre = Self::get_centre_by_raster(dataset)?; + if geometric_centre.x != 0.0f64 || geometric_centre.y != 0.0f64 { + color_eyre::eyre::bail!( + "Tile centre ({},{}) must be at 0,0.", + geometric_centre.x, + geometric_centre.y + ); + } + let projection = &dataset.spatial_ref()?.to_proj4()?; let lat_0: f64 = projection @@ -174,7 +166,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 +179,9 @@ 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_centre_by_raster(&dataset).unwrap(); + assert_eq!(centre, geo::Coord { x: 0.0, y: 0.0 }); } } From 4a8980863e501f1f8305717a432230a5968e8269 Mon Sep 17 00:00:00 2001 From: Thomas Buckley-Houston Date: Wed, 22 Apr 2026 16:00:48 +0100 Subject: [PATCH 2/4] docs: Update Usage in README --- README.md | 52 ++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 44 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index ffbc3d2..5904183 100644 --- a/README.md +++ b/README.md @@ -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 From e3b60bc3683b23187b2b07e43436e285750058b5 Mon Sep 17 00:00:00 2001 From: Thomas Buckley-Houston Date: Fri, 24 Apr 2026 14:44:37 +0100 Subject: [PATCH 3/4] feat: Allow <1mm off-centre tiles When creating tiles that are not at exactly 100m resolution their centres are often off by a few fractions of a milimetre. Which should be allowed. --- benchmarks/samples/aeqd_bad_centre.tiff | Bin 0 -> 1634 bytes crates/total-viewsheds/src/tile.rs | 22 +++++++++++++++++----- 2 files changed, 17 insertions(+), 5 deletions(-) create mode 100644 benchmarks/samples/aeqd_bad_centre.tiff diff --git a/benchmarks/samples/aeqd_bad_centre.tiff b/benchmarks/samples/aeqd_bad_centre.tiff new file mode 100644 index 0000000000000000000000000000000000000000..f9e26115f212558a10913b7cfcaeaf2341577864 GIT binary patch literal 1634 zcmeH{%TE(g6vlruQ_7LPW$z?I3EcRz!I*1Q$oQc{yB;WZCzEXluyURFev>Ovo@yiQQ^6)Nrmw1!i(ceq`*BaX8z9F9WQU_zbOYf8AEoHPZyL; zWHW=MTB9PDl`MW8KW?k6Sj?(?KqEX#-Fe!&Fsu4G&ql@v#-fQ(fA746s+%Qk%jDn6 zrCp_*cL1AeYFH&I)KSlB*3dwZMw(blGwW!fmGx{O#6~vJMmw9?LI+#vq>F8AX9qjk zMK`S31Xb&6md>7#2JP;%Q+H^aGp`d z7-xbDOmdM+T;>W_nc^DLE6goEQO9#ji*V(ZxSaJjv59hva>x_0m2ykB@b}7MG72i! zyDCfMs%-jlS@ljq?d6P8yF&s~tNp6!SJnNl>Sa`G^~D#q>b;yXU9>dX_%wF5RoTid z76FYuYBieFX)Id#&G(o0zUE`j=Z}5O`N{li^+~1Ucb9p&yXl+hTXX68=wKo=G?hw* z!Xa;Icy=b0OvyrYHknM%%}n2zSFTiJ!qM;(3(+5+PA0E&PR8+SR@NQTESHDpiXP#f kE9%RNu6QoobHyj&YW#P_zqw-JcDU+8B$Id1(T6FIA4>(d0RR91 literal 0 HcmV?d00001 diff --git a/crates/total-viewsheds/src/tile.rs b/crates/total-viewsheds/src/tile.rs index 0da2886..964214b 100644 --- a/crates/total-viewsheds/src/tile.rs +++ b/crates/total-viewsheds/src/tile.rs @@ -98,7 +98,7 @@ impl Tile { } /// Get the metric centre of the tile simply by dividing the extent in half. - fn get_centre_by_raster(dataset: &gdal::Dataset) -> Result { + fn get_metric_centre(dataset: &gdal::Dataset) -> Result { let (width, height) = dataset.raster_size(); #[expect( @@ -127,10 +127,13 @@ impl Tile { /// Get the lat/lon centre of the tile by querying the projection's definition. fn get_centre_by_projection(dataset: &gdal::Dataset) -> Result { - let geometric_centre = Self::get_centre_by_raster(dataset)?; - if geometric_centre.x != 0.0f64 || geometric_centre.y != 0.0f64 { + 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 at 0,0.", + "Tile centre ({},{}) must be within {ALLOWED_DEVIATION} of 0,0.", geometric_centre.x, geometric_centre.y ); @@ -181,7 +184,16 @@ mod test { #[test] fn get_metric_centre() { let dataset = gdal::Dataset::open("../../benchmarks/samples/aeqd_10x10.tiff").unwrap(); - let centre = Tile::get_centre_by_raster(&dataset).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}'"); + } } From 7967ab8fb6ec34bcee497b48876346ea46175e66 Mon Sep 17 00:00:00 2001 From: Thomas Buckley-Houston Date: Fri, 24 Apr 2026 14:46:05 +0100 Subject: [PATCH 4/4] feat: Remove `--scale` CLI argument We support GeoTiffs now so the scale always comes from the tile. --- benchmarks/run.sh | 1 - crates/total-viewsheds/src/config.rs | 5 +-- .../total-viewsheds/src/cpu/unrolled_los.rs | 3 +- crates/total-viewsheds/src/dem.rs | 27 +++------------- crates/total-viewsheds/src/main.rs | 32 +++++++++---------- crates/total-viewsheds/src/run/compute.rs | 2 +- 6 files changed, 24 insertions(+), 46 deletions(-) diff --git a/benchmarks/run.sh b/benchmarks/run.sh index 80cffad..eba1776 100755 --- a/benchmarks/run.sh +++ b/benchmarks/run.sh @@ -32,7 +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 \ --rings-per-km 3 \ --backend "$backend" \ --process all \ diff --git a/crates/total-viewsheds/src/config.rs b/crates/total-viewsheds/src/config.rs index c530874..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, @@ -133,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, })