From f1f4797ebf3e6bb30a03a0736d1a170c5864d9db Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Mon, 7 Oct 2024 20:49:16 +0200 Subject: [PATCH 01/15] feat: Determine correct coordinate transform from TIFF tags --- src/coordinate_transform.rs | 90 +++++++++++++++++++++++++++++++++++++ src/decoder_ext.rs | 31 +++++++++++++ src/lib.rs | 5 +++ 3 files changed, 126 insertions(+) create mode 100644 src/coordinate_transform.rs diff --git a/src/coordinate_transform.rs b/src/coordinate_transform.rs new file mode 100644 index 0000000..fdcf10a --- /dev/null +++ b/src/coordinate_transform.rs @@ -0,0 +1,90 @@ +use tiff::{TiffError, TiffFormatError, TiffResult}; + +const MODEL_TIE_POINT_TAG: &str = "ModelTiePointTag"; +const MODEL_PIXEL_SCALE_TAG: &str = "ModelPixelScaleTag"; +const MODEL_TRANSFORMATION_TAG: &str = "ModelTransformationTag"; + +#[derive(Debug)] +pub(super) enum CoordinateTransform { + AffineTransform, + TiePointAndPixelScale, + TiePoints, +} + +impl CoordinateTransform { + pub(super) fn from_tag_data( + pixel_scale_data: Option>, + model_tie_points_data: Option>, + model_transformation_data: Option>, + ) -> TiffResult { + let pixel_scale = pixel_scale_data + .map(|data| { + <[f64; 3]>::try_from(data).map_err(|_| { + TiffError::FormatError(TiffFormatError::Format(format!( + "Number values in {MODEL_PIXEL_SCALE_TAG} must be equal to 3" + ))) + }) + }) + .transpose()?; + let tie_points = model_tie_points_data + .map(|data| { + let len = data.len(); + if len == 0 { + return Err(TiffError::FormatError(TiffFormatError::Format(format!( + "Number of values in {MODEL_TIE_POINT_TAG} must be greater than 0" + )))); + } + + if len % 6 != 0 { + return Err(TiffError::FormatError(TiffFormatError::Format(format!( + "Number of values in {MODEL_TIE_POINT_TAG} must be divisible by 6" + )))); + } + + Ok(data) + }) + .transpose()?; + let transformation_matrix = model_transformation_data + .map(|data| { + <[f64; 16]>::try_from(data).map_err(|_| { + TiffError::FormatError(TiffFormatError::Format(format!( + "Number of values in {MODEL_TRANSFORMATION_TAG} must be equal to 16" + ))) + }) + }) + .transpose()?; + + if let Some(transformation_matrix) = transformation_matrix { + if pixel_scale.is_some() { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_PIXEL_SCALE_TAG} must not be specified when {MODEL_TRANSFORMATION_TAG} is present"), + ))); + } + if tie_points.is_some() { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_TIE_POINT_TAG} must not be specified when {MODEL_TRANSFORMATION_TAG} is present"), + ))); + } + + Ok(CoordinateTransform::AffineTransform) + } else { + let Some(tie_points) = tie_points else { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_TIE_POINT_TAG} must be present when {MODEL_TRANSFORMATION_TAG} is missing"), + ))); + }; + + if tie_points.len() == 6 { + let Some(pixel_scale) = pixel_scale else { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_PIXEL_SCALE_TAG} must be specified when {MODEL_TIE_POINT_TAG} contains 6 values"), + ))); + }; + + Ok(Self::TiePointAndPixelScale) + } else { + Ok(Self::TiePoints) + } + } + } +} diff --git a/src/decoder_ext.rs b/src/decoder_ext.rs index 5ea1a83..adfd560 100644 --- a/src/decoder_ext.rs +++ b/src/decoder_ext.rs @@ -4,13 +4,44 @@ use tiff::decoder::Decoder; use tiff::tags::Tag; use tiff::TiffResult; +use crate::coordinate_transform::CoordinateTransform; use crate::geo_key_directory::GeoKeyDirectory; pub(super) trait DecoderExt { + fn coordinate_transform(&mut self) -> TiffResult>; + fn geo_key_directory(&mut self) -> TiffResult; } impl DecoderExt for Decoder { + fn coordinate_transform(&mut self) -> TiffResult> { + let pixel_scale_data = self + .find_tag(Tag::ModelPixelScaleTag)? + .map(|value| value.into_f64_vec()) + .transpose()?; + let tie_points_data = self + .find_tag(Tag::ModelTiepointTag)? + .map(|value| value.into_f64_vec()) + .transpose()?; + let model_transformation_data = self + .find_tag(Tag::ModelTransformationTag)? + .map(|value| value.into_f64_vec()) + .transpose()?; + + if pixel_scale_data.is_none() + && tie_points_data.is_none() + && model_transformation_data.is_none() + { + return Ok(None); + } + + Ok(Some(CoordinateTransform::from_tag_data( + pixel_scale_data, + tie_points_data, + model_transformation_data, + )?)) + } + fn geo_key_directory(&mut self) -> TiffResult { let Some(directory_data) = self .find_tag(Tag::GeoKeyDirectoryTag)? diff --git a/src/lib.rs b/src/lib.rs index 6ac7137..a51f758 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -9,9 +9,11 @@ use tiff::TiffResult; pub use crate::geo_key_directory::*; +use crate::coordinate_transform::*; use crate::decoder_ext::*; use crate::raster_data::*; +mod coordinate_transform; mod decoder_ext; mod geo_key_directory; mod raster_data; @@ -39,6 +41,7 @@ pub struct GeoTiff { pub raster_width: usize, pub raster_height: usize, pub num_samples: usize, + coordinate_transform: Option, raster_data: RasterData, } @@ -47,6 +50,7 @@ impl GeoTiff { let mut decoder = Decoder::new(reader)?; let geo_key_directory = decoder.geo_key_directory()?; + let coordinate_transform = decoder.coordinate_transform()?; let (raster_width, raster_height) = decoder .dimensions() @@ -74,6 +78,7 @@ impl GeoTiff { raster_width, raster_height, num_samples, + coordinate_transform, raster_data, }) } From 3e386625a8ef66ed079bfe7c1ea36256a5d62059 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Mon, 7 Oct 2024 20:54:17 +0200 Subject: [PATCH 02/15] refactor: Make raster type an enum --- src/geo_key_directory.rs | 18 +++++++++++++++--- tests/integration.rs | 4 ++-- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/geo_key_directory.rs b/src/geo_key_directory.rs index c74b7df..f23b6a8 100644 --- a/src/geo_key_directory.rs +++ b/src/geo_key_directory.rs @@ -12,7 +12,7 @@ pub struct GeoKeyDirectory { pub key_revision: u16, pub minor_revision: u16, pub model_type: Option, - pub raster_type: Option, + pub raster_type: Option, pub citation: Option, pub geographic_type: Option, pub geog_citation: Option, @@ -99,8 +99,13 @@ impl GeoKeyDirectory { Self::get_short(key_tag, location_tag, *count, *value_or_offset)?.into() } GeoKeyDirectoryTag::RasterType => { - directory.raster_type = - Self::get_short(key_tag, location_tag, *count, *value_or_offset)?.into() + let raster_type = + Self::get_short(key_tag, location_tag, *count, *value_or_offset)?; + directory.raster_type = Some(RasterType::try_from(raster_type).map_err(|_| { + TiffError::FormatError(TiffFormatError::Format(format!( + "Unknown raster type: {raster_type}" + ))) + })?) } GeoKeyDirectoryTag::Citation => { directory.citation = Self::get_string( @@ -657,3 +662,10 @@ enum GeoKeyDirectoryTag { VerticalDatum = 4098, VerticalUnits = 4099, } + +#[derive(Debug, Clone, Copy, PartialEq, TryFromPrimitive, IntoPrimitive)] +#[repr(u16)] +pub enum RasterType { + RasterPixelIsArea = 1, + RasterPixelIsPoint = 2, +} diff --git a/tests/integration.rs b/tests/integration.rs index 2a2181f..a766306 100644 --- a/tests/integration.rs +++ b/tests/integration.rs @@ -1,7 +1,7 @@ use std::fs::File; use std::path::Path; -use geotiff::{GeoKeyDirectory, GeoTiff}; +use geotiff::{GeoKeyDirectory, GeoTiff, RasterType}; fn read_geotiff>(path: P) -> GeoTiff { GeoTiff::read(File::open(path).expect("File I/O error")).expect("File I/O error") @@ -56,7 +56,7 @@ fn test_load_merc() { key_revision: 1, minor_revision: 2, model_type: Some(1), - raster_type: Some(1), + raster_type: Some(RasterType::RasterPixelIsArea), geog_geodetic_datum: Some(6267), geog_ellipsoid: Some(7008), projected_type: Some(32767), From 60a8a999707f5d6c9b02ece255a72c0e67d39be2 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Mon, 7 Oct 2024 21:05:31 +0200 Subject: [PATCH 03/15] feat: Implement transformation by tie point and pixel scale --- Cargo.lock | 488 ++++++++++++++++++ Cargo.toml | 4 + ...ie_point_and_pixel_scale_pixel_is_area.tif | Bin 0 -> 2374 bytes ...e_point_and_pixel_scale_pixel_is_point.tif | Bin 0 -> 2374 bytes src/coordinate_transform.rs | 43 +- .../tie_point_and_pixel_scale.rs | 50 ++ src/lib.rs | 83 ++- tests/integration.rs | 260 +++++++++- 8 files changed, 906 insertions(+), 22 deletions(-) create mode 100644 resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif create mode 100644 resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif create mode 100644 src/coordinate_transform/tie_point_and_pixel_scale.rs diff --git a/Cargo.lock b/Cargo.lock index 8822144..4ef4f85 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -8,18 +8,103 @@ version = "2.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "512761e0bb2578dd7380c6baaa0f4ce03e84f95e960231d1dec8bf4d7d6e2627" +[[package]] +name = "aho-corasick" +version = "1.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916" +dependencies = [ + "memchr", +] + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + [[package]] name = "autocfg" version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" +[[package]] +name = "bindgen" +version = "0.66.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2b84e06fc203107bfbad243f4aba2af864eb7db3b1cf46ea0a023b0b433d2a7" +dependencies = [ + "bitflags", + "cexpr", + "clang-sys", + "lazy_static", + "lazycell", + "log", + "peeking_take_while", + "prettyplease", + "proc-macro2", + "quote", + "regex", + "rustc-hash", + "shlex", + "syn", + "which", +] + +[[package]] +name = "bitflags" +version = "2.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b048fb63fd8b5923fc5aa7b340d8e156aec7ec02f0c78fa8a6ddc2613f6f71de" + +[[package]] +name = "cc" +version = "1.1.28" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2e80e3b6a3ab07840e1cae9b0666a63970dc28e8ed5ffbcdacbfc760c281bfc1" +dependencies = [ + "shlex", +] + +[[package]] +name = "cexpr" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6fac387a98bb7c37292057cffc56d62ecb629900026402633ae9160df93a8766" +dependencies = [ + "nom", +] + [[package]] name = "cfg-if" version = "1.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" +[[package]] +name = "clang-sys" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b023947811758c97c59bf9d1c188fd619ad4718dcaa767947df1cadb14f39f4" +dependencies = [ + "glob", + "libc", + "libloading", +] + +[[package]] +name = "cmake" +version = "0.1.51" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb1e43aa7fd152b1f968787f7dbcdeb306d1867ff373c69955211876c053f91a" +dependencies = [ + "cc", +] + [[package]] name = "crc32fast" version = "1.4.2" @@ -29,12 +114,40 @@ dependencies = [ "cfg-if", ] +[[package]] +name = "either" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" + [[package]] name = "equivalent" version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5443807d6dff69373d433ab9ef5378ad8df50ca6298caf15de6e52e24aaf54d5" +[[package]] +name = "errno" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "534c5cf6194dfab3db3242765c03bbe257cf92f22b38f6bc0c58d59108a820ba" +dependencies = [ + "libc", + "windows-sys 0.52.0", +] + +[[package]] +name = "filetime" +version = "0.2.25" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "35c0522e981e68cbfa8c3f978441a5f34b30b96e146b33cd3359176b50fe8586" +dependencies = [ + "cfg-if", + "libc", + "libredox", + "windows-sys 0.59.0", +] + [[package]] name = "flate2" version = "1.0.33" @@ -45,21 +158,49 @@ dependencies = [ "miniz_oxide", ] +[[package]] +name = "geo-types" +version = "0.7.13" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ff16065e5720f376fbced200a5ae0f47ace85fd70b7e54269790281353b6d61" +dependencies = [ + "approx", + "num-traits", + "serde", +] + [[package]] name = "geotiff" version = "0.0.2" dependencies = [ + "geo-types", "num-traits", "num_enum", + "proj", "tiff", ] +[[package]] +name = "glob" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2fabcfbdc87f4758337ca535fb41a6d701b65693ce38287d856d1674551ec9b" + [[package]] name = "hashbrown" version = "0.14.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" +[[package]] +name = "home" +version = "0.5.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3d1354bf6b7235cb4a0576c2619fd4ed18183f689b12b006a0ee7329eeff9a5" +dependencies = [ + "windows-sys 0.52.0", +] + [[package]] name = "indexmap" version = "2.5.0" @@ -76,12 +217,75 @@ version = "0.3.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f5d4a7da358eff58addd2877a45865158f0d78c911d43a5784ceb7bbf52833b0" +[[package]] +name = "lazy_static" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" + +[[package]] +name = "lazycell" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "830d08ce1d1d941e6b30645f1a0eb5643013d835ce3779a5fc208261dbe10f55" + +[[package]] +name = "libc" +version = "0.2.159" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "561d97a539a36e26a9a5fad1ea11a3039a67714694aaa379433e580854bc3dc5" + +[[package]] +name = "libloading" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4979f22fdb869068da03c9f7528f8297c6fd2606bc3a4affe42e6a823fdb8da4" +dependencies = [ + "cfg-if", + "windows-targets", +] + +[[package]] +name = "libm" +version = "0.2.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" + +[[package]] +name = "libredox" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c0ff37bd590ca25063e35af745c343cb7a0271906fb7b37e4813e8f79f00268d" +dependencies = [ + "bitflags", + "libc", + "redox_syscall", +] + +[[package]] +name = "linux-raw-sys" +version = "0.4.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78b3ae25bc7c8c38cec158d1f2757ee79e9b3740fbc7ccf0e59e4b08d793fa89" + +[[package]] +name = "log" +version = "0.4.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a7a70ba024b9dc04c27ea2f0c0548feb474ec5c54bba33a7f72f873a39d07b24" + [[package]] name = "memchr" version = "2.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" +[[package]] +name = "minimal-lexical" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "68354c5c6bd36d73ff3feceb05efa59b6acb7626617f4962be322a825e61f79a" + [[package]] name = "miniz_oxide" version = "0.8.0" @@ -91,6 +295,16 @@ dependencies = [ "adler2", ] +[[package]] +name = "nom" +version = "7.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d273983c5a657a70a3e8f2a01329822f3b8c8172b73826411a55751e404a0a4a" +dependencies = [ + "memchr", + "minimal-lexical", +] + [[package]] name = "num-traits" version = "0.2.19" @@ -98,6 +312,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", + "libm", ] [[package]] @@ -121,6 +336,34 @@ dependencies = [ "syn", ] +[[package]] +name = "once_cell" +version = "1.20.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1261fe7e33c73b354eab43b1273a57c8f967d0391e80353e51f764ac02cf6775" + +[[package]] +name = "peeking_take_while" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19b17cddbe7ec3f8bc800887bab5e717348c95ea2ca0b1bf0837fb964dc67099" + +[[package]] +name = "pkg-config" +version = "0.3.31" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "953ec861398dccce10c670dfeaf3ec4911ca479e9c02154b3a215178c5f566f2" + +[[package]] +name = "prettyplease" +version = "0.2.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "479cf940fbbb3426c32c5d5176f62ad57549a0bb84773423ba8be9d089f5faba" +dependencies = [ + "proc-macro2", + "syn", +] + [[package]] name = "proc-macro-crate" version = "3.2.0" @@ -139,6 +382,32 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "proj" +version = "0.27.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ad1830ad8966eba22c76e78440458f07bd812bef5c3efdf335dec55cd1085ab" +dependencies = [ + "geo-types", + "libc", + "num-traits", + "proj-sys", + "thiserror", +] + +[[package]] +name = "proj-sys" +version = "0.23.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "601bf4fa1e17fde1a56d303f7bed5c65969cf1822c6baf5d6c2c12c593638fec" +dependencies = [ + "bindgen", + "cmake", + "flate2", + "pkg-config", + "tar", +] + [[package]] name = "quote" version = "1.0.37" @@ -148,6 +417,89 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "redox_syscall" +version = "0.5.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b6dfecf2c74bce2466cabf93f6664d6998a69eb21e39f4207930065b27b771f" +dependencies = [ + "bitflags", +] + +[[package]] +name = "regex" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "38200e5ee88914975b69f657f0801b6f6dccafd44fd9326302a4aaeecfacb1d8" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "368758f23274712b504848e9d5a6f010445cc8b87a7cdb4d7cbee666c1288da3" +dependencies = [ + "aho-corasick", + "memchr", + "regex-syntax", +] + +[[package]] +name = "regex-syntax" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "rustix" +version = "0.38.37" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8acb788b847c24f28525660c4d7758620a7210875711f79e7f663cc152726811" +dependencies = [ + "bitflags", + "errno", + "libc", + "linux-raw-sys", + "windows-sys 0.52.0", +] + +[[package]] +name = "serde" +version = "1.0.210" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c8e3592472072e6e22e0a54d5904d9febf8508f65fb8552499a1abc7d1078c3a" +dependencies = [ + "serde_derive", +] + +[[package]] +name = "serde_derive" +version = "1.0.210" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "243902eda00fad750862fc144cea25caca5e20d615af0a81bee94ca738f1df1f" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "shlex" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" + [[package]] name = "syn" version = "2.0.77" @@ -159,6 +511,37 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "tar" +version = "0.4.42" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ff6c40d3aedb5e06b57c6f669ad17ab063dd1e63d977c6a88e7f4dfa4f04020" +dependencies = [ + "filetime", + "libc", + "xattr", +] + +[[package]] +name = "thiserror" +version = "1.0.64" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d50af8abc119fb8bb6dbabcfa89656f46f84aa0ac7688088608076ad2b459a84" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.64" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08904e7672f5eb876eaaf87e0ce17857500934f4981c4a0ab2b4aa98baac7fc3" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + [[package]] name = "tiff" version = "0.9.1" @@ -199,6 +582,100 @@ version = "0.1.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "53a85b86a771b1c87058196170769dd264f66c0782acf1ae6cc51bfd64b39082" +[[package]] +name = "which" +version = "4.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "87ba24419a2078cd2b0f2ede2691b6c66d8e47836da3b6db8265ebad47afbfc7" +dependencies = [ + "either", + "home", + "once_cell", + "rustix", +] + +[[package]] +name = "windows-sys" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "282be5f36a8ce781fad8c8ae18fa3f9beff57ec1b52cb3de0789201425d9a33d" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-sys" +version = "0.59.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-targets" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_gnullvm", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469" + +[[package]] +name = "windows_i686_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b" + +[[package]] +name = "windows_i686_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66" + +[[package]] +name = "windows_i686_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" + [[package]] name = "winnow" version = "0.6.18" @@ -207,3 +684,14 @@ checksum = "68a9bda4691f099d435ad181000724da8e5899daa10713c2d432552b9ccd3a6f" dependencies = [ "memchr", ] + +[[package]] +name = "xattr" +version = "1.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8da84f1a25939b27f6820d92aed108f83ff920fdf11a7b19366c27c4cda81d4f" +dependencies = [ + "libc", + "linux-raw-sys", + "rustix", +] diff --git a/Cargo.toml b/Cargo.toml index 70fa012..b4eba9c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,6 +7,10 @@ authors = ["Dominik Bucher "] repository = "https://github.com/georust/geotiff" [dependencies] +geo-types = "0.7" num_enum = "0.7" num-traits = "0.2" tiff = "0.9" + +[dev-dependencies] +proj = "0.27" diff --git a/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif b/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif new file mode 100644 index 0000000000000000000000000000000000000000..b18b0b0fbe1d258542e7724c939b39b1f18b684c GIT binary patch literal 2374 zcmebEWzb?^VBla7WMp7q2C^6#e}f1Jn_(Z2%>-q00NKnCcEcthn++-s(j>^h0;G6= z_!grmR2*m^qZpJe0A!29)G%@}N-?m4v;mP7qcoiD#HfL!*8#)^JISk=4{SD&!|v76 z1ZJ}V*-7o|!3uz;vrk&m$iM_Nlf^*K07!uVJJ6Nvw}AL3(1S8SY!Ae7KwJyN3xN1A z5IdcaL|c3?l;oIzN&I literal 0 HcmV?d00001 diff --git a/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif b/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif new file mode 100644 index 0000000000000000000000000000000000000000..2bb7fa3b214540ee9ce104067914037276448743 GIT binary patch literal 2374 zcmebEWzb?^VBla7WMp7q2C^6#e}f1Jn_(Z2%>-q00NKnCcEcthn++-s(j>^h0;G6= z_!grmR2*m^qZpJe0A!29)G%@}N-?m4v;mP7qcoiD#HfL!*8#)^JISk=4{SD&!|v76 z1ZJ}V*-7o|!3uz;vrk&m$iM_Nlf^*K07!uVJJ6Nvw}AL3(1S8SY!Ae7KwJyN3xN1A z5IdcaL|c3?l;oJKT~7 literal 0 HcmV?d00001 diff --git a/src/coordinate_transform.rs b/src/coordinate_transform.rs index fdcf10a..13a7a8f 100644 --- a/src/coordinate_transform.rs +++ b/src/coordinate_transform.rs @@ -1,5 +1,8 @@ +use geo_types::Coord; use tiff::{TiffError, TiffFormatError, TiffResult}; +mod tie_point_and_pixel_scale; + const MODEL_TIE_POINT_TAG: &str = "ModelTiePointTag"; const MODEL_PIXEL_SCALE_TAG: &str = "ModelPixelScaleTag"; const MODEL_TRANSFORMATION_TAG: &str = "ModelTransformationTag"; @@ -7,7 +10,11 @@ const MODEL_TRANSFORMATION_TAG: &str = "ModelTransformationTag"; #[derive(Debug)] pub(super) enum CoordinateTransform { AffineTransform, - TiePointAndPixelScale, + TiePointAndPixelScale { + raster_point: Coord, + model_point: Coord, + pixel_scale: Coord, + }, TiePoints, } @@ -81,10 +88,42 @@ impl CoordinateTransform { ))); }; - Ok(Self::TiePointAndPixelScale) + Self::from_tie_point_and_pixel_scale(&tie_points, &pixel_scale) } else { Ok(Self::TiePoints) } } } + + pub fn transform_to_model(&self, coord: &Coord) -> Coord { + match self { + CoordinateTransform::TiePointAndPixelScale { + raster_point, + model_point, + pixel_scale, + } => Self::transform_to_model_by_tie_point_and_pixel_scale( + raster_point, + model_point, + pixel_scale, + coord, + ), + _ => unimplemented!() + } + } + + pub(super) fn transform_to_raster(&self, coord: &Coord) -> Coord { + match self { + CoordinateTransform::TiePointAndPixelScale { + raster_point, + model_point, + pixel_scale, + } => Self::transform_to_raster_by_tie_point_and_pixel_scale( + raster_point, + model_point, + pixel_scale, + coord, + ), + _ => unimplemented!() + } + } } diff --git a/src/coordinate_transform/tie_point_and_pixel_scale.rs b/src/coordinate_transform/tie_point_and_pixel_scale.rs new file mode 100644 index 0000000..6075752 --- /dev/null +++ b/src/coordinate_transform/tie_point_and_pixel_scale.rs @@ -0,0 +1,50 @@ +use geo_types::Coord; +use tiff::TiffResult; + +use crate::coordinate_transform::CoordinateTransform; + +impl CoordinateTransform { + pub(super) fn from_tie_point_and_pixel_scale( + tie_points: &[f64], + pixel_scale: &[f64], + ) -> TiffResult { + Ok(CoordinateTransform::TiePointAndPixelScale { + raster_point: Coord { + x: tie_points[0], + y: tie_points[1], + }, + model_point: Coord { + x: tie_points[3], + y: tie_points[4], + }, + pixel_scale: Coord { + x: pixel_scale[0], + y: pixel_scale[1], + }, + }) + } + + pub(super) fn transform_to_model_by_tie_point_and_pixel_scale( + raster_point: &Coord, + model_point: &Coord, + pixel_scale: &Coord, + coord: &Coord, + ) -> Coord { + Coord { + x: (coord.x - raster_point.x) * pixel_scale.x + model_point.x, + y: (coord.y - raster_point.y) * -pixel_scale.y + model_point.y, + } + } + + pub(super) fn transform_to_raster_by_tie_point_and_pixel_scale( + raster_point: &Coord, + model_point: &Coord, + pixel_scale: &Coord, + coord: &Coord, + ) -> Coord { + Coord { + x: (coord.x - model_point.x) / pixel_scale.x + raster_point.x, + y: (coord.y - model_point.y) / -pixel_scale.y + raster_point.y, + } + } +} diff --git a/src/lib.rs b/src/lib.rs index a51f758..438cb22 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,7 +1,7 @@ //! A [GeoTIFF](https://www.ogc.org/standard/geotiff) library for Rust use std::any::type_name; use std::io::{Read, Seek}; - +use geo_types::{Coord, Rect}; use num_traits::FromPrimitive; use tiff::decoder::{Decoder, DecodingResult}; use tiff::tags::Tag; @@ -83,23 +83,35 @@ impl GeoTiff { }) } - pub fn get_value_at(&self, x: usize, y: usize, sample: usize) -> T { - let GeoTiff { - raster_width, - num_samples, - raster_data, - .. - } = self; + pub fn model_extent(&self) -> Rect { + let offset = self.raster_offset(); + let lower = Coord { + x: offset, + y: offset, + }; + let upper = Coord { + x: self.raster_width as f64 + offset, + y: self.raster_height as f64 + offset, + }; - if &sample >= num_samples { - panic!( - "sample out of bounds: the number of samples is {} but the sample is {}", - num_samples, sample + if let Some(coordinate_transform) = &self.coordinate_transform { + Rect::new( + coordinate_transform.transform_to_model(&lower), + coordinate_transform.transform_to_model(&upper), ) + } else { + Rect::new(lower, upper) } + } + + pub fn get_value_at( + &self, + coord: &Coord, + sample: usize, + ) -> Option { + let index = self.compute_index(coord, sample)?; - let index = (y * raster_width + x) * num_samples + sample; - match raster_data { + Some(match &self.raster_data { RasterData::U8(data) => unwrap_primitive_type!(T::from_u8(data[index]), u8, T), RasterData::U16(data) => unwrap_primitive_type!(T::from_u16(data[index]), u16, T), RasterData::U32(data) => unwrap_primitive_type!(T::from_u32(data[index]), u32, T), @@ -110,6 +122,49 @@ impl GeoTiff { RasterData::I16(data) => unwrap_primitive_type!(T::from_i16(data[index]), i16, T), RasterData::I32(data) => unwrap_primitive_type!(T::from_i32(data[index]), i32, T), RasterData::I64(data) => unwrap_primitive_type!(T::from_i64(data[index]), i64, T), + }) + } + + fn compute_index(&self, coord: &Coord, sample: usize) -> Option { + let GeoTiff { + raster_width, + raster_height, + num_samples, + coordinate_transform, + .. + } = self; + + if &sample >= num_samples { + panic!( + "sample out of bounds: the number of samples is {} but the sample is {}", + num_samples, sample + ) + } + + let mut coord = match coordinate_transform { + None => *coord, + Some(transform) => transform.transform_to_raster(coord), + }; + + let raster_offset = self.raster_offset(); + coord.x -= raster_offset; + coord.y -= raster_offset; + + if coord.x < 0.0 + || coord.x >= *raster_width as f64 + || coord.y < 0.0 + || coord.y >= *raster_height as f64 + { + return None; + } + + Some((coord.y as usize * raster_width + coord.x as usize) * num_samples + sample) + } + + fn raster_offset(&self) -> f64 { + match self.geo_key_directory.raster_type { + Some(RasterType::RasterPixelIsPoint) => -0.5, + _ => 0.0, } } } diff --git a/tests/integration.rs b/tests/integration.rs index a766306..efdb2de 100644 --- a/tests/integration.rs +++ b/tests/integration.rs @@ -1,7 +1,68 @@ use std::fs::File; use std::path::Path; +use geo_types::{Coord, Rect}; use geotiff::{GeoKeyDirectory, GeoTiff, RasterType}; +use proj::Proj; + +trait RectExt { + fn rounded(self) -> Self; +} + +impl RectExt for Rect { + fn rounded(mut self) -> Self { + let factor = 10f64.powi(8); + let mut min = self.min(); + let mut max = self.max(); + min.x = (min.x * factor).round() / factor; + min.y = (min.y * factor).round() / factor; + max.x = (max.x * factor).round() / factor; + max.y = (max.y * factor).round() / factor; + self.set_min(min); + self.set_max(max); + self + } +} + +const BREGENZ: Coord = Coord { + x: 9.74926, + y: 47.50315, +}; +const EISENSTADT: Coord = Coord { + x: 15.43301, + y: 47.06298, +}; +const GRAZ: Coord = Coord { + x: 15.43301, + y: 47.06298, +}; +const INNSBRUCK: Coord = Coord { + x: 11.39960, + y: 47.26239, +}; +const KLAGENFURT: Coord = Coord { + x: 14.31528, + y: 46.62366, +}; +const LINZ: Coord = Coord { + x: 14.30571, + y: 48.27532, +}; +const SALZBURG: Coord = Coord { + x: 13.05345, + y: 47.80763, +}; +const SANKT_POELTEN: Coord = Coord { + x: 15.62291, + y: 48.20440, +}; +const VIENNA: Coord = Coord { + x: 16.37499, + y: 48.22158, +}; + +const WHITE: u8 = 255; +const BLACK: u8 = 0; fn read_geotiff>(path: P) -> GeoTiff { GeoTiff::read(File::open(path).expect("File I/O error")).expect("File I/O error") @@ -15,9 +76,28 @@ fn test_load_marbles() { assert_eq!(geotiff.raster_width, 1419); assert_eq!(geotiff.raster_height, 1001); assert_eq!(geotiff.num_samples, 3); - assert_eq!(geotiff.get_value_at::(761, 599, 0), 147); - assert_eq!(geotiff.get_value_at::(761, 599, 1), 128); - assert_eq!(geotiff.get_value_at::(761, 599, 2), 165); + assert_eq!( + geotiff.model_extent(), + Rect::new( + Coord { x: 0.0, y: 0.0 }, + Coord { + x: 1419.0, + y: 1001.0 + } + ) + ); + assert_eq!( + geotiff.get_value_at::(&Coord { x: 761.0, y: 599.0 }, 0), + Some(147) + ); + assert_eq!( + geotiff.get_value_at::(&Coord { x: 761.0, y: 599.0 }, 1), + Some(128) + ); + assert_eq!( + geotiff.get_value_at::(&Coord { x: 761.0, y: 599.0 }, 2), + Some(165) + ); } #[test] @@ -28,9 +108,49 @@ fn test_load_zh_dem_25() { assert_eq!(geotiff.raster_width, 399); assert_eq!(geotiff.raster_height, 366); assert_eq!(geotiff.num_samples, 1); - assert_eq!(geotiff.get_value_at::(0, 0, 0), 551); - assert_eq!(geotiff.get_value_at::(67, 45, 0), 530); - assert_eq!(geotiff.get_value_at::(325, 142, 0), 587); + assert_eq!( + geotiff.model_extent(), + Rect::new( + Coord { + x: 677562.5, + y: 243862.5 + }, + Coord { + x: 687537.5, + y: 253012.5 + } + ) + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: 677575.0, + y: 253000.0 + }, + 0 + ), + Some(551) + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: 679250.0, + y: 251875.0 + }, + 0 + ), + Some(530) + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: 685700.0, + y: 249450.0 + }, + 0 + ), + Some(587) + ); assert_eq!( geotiff.geo_key_directory, @@ -73,4 +193,132 @@ fn test_load_merc() { ..Default::default() } ); + + assert_eq!( + geotiff.model_extent(), + Rect::new( + Coord { + x: 1871032.9538880002, + y: 662408.6726400064 + }, + Coord { + x: 1901982.949391994, + y: 693358.6681440001 + } + ) + ); +} + +#[test] +fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + +fn test_transform>(path: P, raster_type: RasterType) { + let geotiff = read_geotiff(path); + + let proj = Proj::new_known_crs( + "EPSG:4326", + &format!("EPSG:{}", geotiff.geo_key_directory.projected_type.unwrap()), + None, + ) + .unwrap(); + + let mut bregenz = proj.project(BREGENZ, false).unwrap(); + let mut eisenstadt = proj.project(EISENSTADT, false).unwrap(); + let mut graz = proj.project(GRAZ, false).unwrap(); + let mut innsbruck = proj.project(INNSBRUCK, false).unwrap(); + let mut klagenfurt = proj.project(KLAGENFURT, false).unwrap(); + let mut linz = proj.project(LINZ, false).unwrap(); + let mut salzburg = proj.project(SALZBURG, false).unwrap(); + let mut sankt_poelten = proj.project(SANKT_POELTEN, false).unwrap(); + let mut vienna = proj.project(VIENNA, false).unwrap(); + let mut expected_model_extent = Rect::new( + Coord { + x: 4302000.0, + y: 2621000.0, + }, + Coord { + x: 4809000.0, + y: 2811000.0, + }, + ); + + // If raster type is RasterPixelIsPoint, shift coordinates accordingly + if raster_type == RasterType::RasterPixelIsPoint { + let mut min = expected_model_extent.min(); + let mut max = expected_model_extent.max(); + + let mut coords = [ + &mut bregenz, + &mut eisenstadt, + &mut graz, + &mut innsbruck, + &mut klagenfurt, + &mut linz, + &mut salzburg, + &mut sankt_poelten, + &mut vienna, + &mut min, + &mut max, + ]; + + for coord in coords.iter_mut() { + coord.x -= 500.0; + coord.y += 500.0; + } + + expected_model_extent.set_min(min); + expected_model_extent.set_max(max); + } + + assert_eq!(geotiff.model_extent().rounded(), expected_model_extent); + + // Non-capital location should return background color + assert_eq!( + geotiff.get_value_at::(&expected_model_extent.center(), 0), + Some(WHITE) + ); + + // A location outside the model extent should return None + + let min = expected_model_extent.min(); + let max = expected_model_extent.max(); + + assert_eq!( + geotiff.get_value_at::(&Coord { x: min.x, y: min.y }, 0), + None + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: max.x + 1.0, + y: max.y + 1.0, + }, + 0 + ), + None + ); + + // Each capital location should return Some(BLACK) + assert_eq!(geotiff.get_value_at::(&bregenz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&eisenstadt, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&graz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&innsbruck, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&klagenfurt, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&linz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&salzburg, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&sankt_poelten, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&vienna, 0), Some(BLACK)); } From 721fdfed147edcea6462c3beca52b72a25a03cf1 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Mon, 7 Oct 2024 21:09:12 +0200 Subject: [PATCH 04/15] feat: Implement transformation by affine transform --- ...als_model_transformation_pixel_is_area.tif | Bin 0 -> 2418 bytes ...ls_model_transformation_pixel_is_point.tif | Bin 0 -> 2418 bytes src/coordinate_transform.rs | 16 ++++-- src/coordinate_transform/affine_transform.rs | 46 ++++++++++++++++++ tests/integration.rs | 16 ++++++ 5 files changed, 75 insertions(+), 3 deletions(-) create mode 100644 resources/austrian_capitals_model_transformation_pixel_is_area.tif create mode 100644 resources/austrian_capitals_model_transformation_pixel_is_point.tif create mode 100644 src/coordinate_transform/affine_transform.rs diff --git a/resources/austrian_capitals_model_transformation_pixel_is_area.tif b/resources/austrian_capitals_model_transformation_pixel_is_area.tif new file mode 100644 index 0000000000000000000000000000000000000000..9a976dc915643652f230ed9eb03b906818d6bb0f GIT binary patch literal 2418 zcmebEWzb?^VBla7U}Rum2C^6#e}f1Jn_(Z2%>-q00NKnCcEcthn++-s(j>^h0;G6= z_#C4sR2*m^qZpJe0A!29)%;_WVqgVn10p>}X*k=AQ3FY@1BeZFl3nW!u-QOU*zMZa zg9SnCeM=e{n1JT680Z-QDG*@42lNaZ5UT>Q2M}iiaTgG;1mZJ5{GNe-q00NKnCcEcthn++-s(j>^h0;G6= z_#C4sR2*m^qZpJe0A!29)%;_WVqgVn10p>}X*k=AQ3FY@1BeZFl3nW!u-QOU*zMZa zg9SnCeM=e{n1JT680Z-QDG*@42lNaZ5UT>Q2M}iiaTgG;1mZJ5{GNe Coord { match self { + CoordinateTransform::AffineTransform { transform, .. } => { + Self::transform_by_affine_transform(transform, coord) + } CoordinateTransform::TiePointAndPixelScale { raster_point, model_point, @@ -113,6 +120,9 @@ impl CoordinateTransform { pub(super) fn transform_to_raster(&self, coord: &Coord) -> Coord { match self { + CoordinateTransform::AffineTransform { + inverse_transform, .. + } => Self::transform_by_affine_transform(inverse_transform, coord), CoordinateTransform::TiePointAndPixelScale { raster_point, model_point, diff --git a/src/coordinate_transform/affine_transform.rs b/src/coordinate_transform/affine_transform.rs new file mode 100644 index 0000000..d21bd46 --- /dev/null +++ b/src/coordinate_transform/affine_transform.rs @@ -0,0 +1,46 @@ +use geo_types::Coord; +use tiff::{TiffError, TiffFormatError, TiffResult}; + +use crate::coordinate_transform::CoordinateTransform; + +impl CoordinateTransform { + pub fn from_transformation_matrix(transformation_matrix: [f64; 16]) -> TiffResult { + let mut transform = [0f64; 6]; + transform[0] = transformation_matrix[0]; + transform[1] = transformation_matrix[1]; + transform[2] = transformation_matrix[3]; + transform[3] = transformation_matrix[4]; + transform[4] = transformation_matrix[5]; + transform[5] = transformation_matrix[7]; + + let det = transform[0] * transform[4] - transform[1] * transform[3]; + if det.abs() < 0.000000000000001 { + return Err(TiffError::FormatError(TiffFormatError::Format( + "Provided transformation matrix is not invertible".into(), + ))); + } + + let mut inverse_transform = [0f64; 6]; + inverse_transform[0] = transform[4] / det; + inverse_transform[1] = -transform[1] / det; + inverse_transform[2] = (transform[1] * transform[5] - transform[2] * transform[4]) / det; + inverse_transform[3] = -transform[3] / det; + inverse_transform[4] = transform[0] / det; + inverse_transform[5] = (-transform[0] * transform[5] + transform[2] * transform[3]) / det; + + Ok(CoordinateTransform::AffineTransform { + transform, + inverse_transform, + }) + } + + pub(super) fn transform_by_affine_transform( + transform: &[f64; 6], + coord: &Coord, + ) -> Coord { + Coord { + x: coord.x * transform[0] + coord.y * transform[1] + transform[2], + y: coord.x * transform[3] + coord.y * transform[4] + transform[5], + } + } +} diff --git a/tests/integration.rs b/tests/integration.rs index efdb2de..27ad2e6 100644 --- a/tests/integration.rs +++ b/tests/integration.rs @@ -225,6 +225,22 @@ fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_point() { ) } +#[test] +fn test_transform_by_model_transformation_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_transformation_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_transformation_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_transformation_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + fn test_transform>(path: P, raster_type: RasterType) { let geotiff = read_geotiff(path); From 2d7e565c99433b9d885125f9ea65777fc4275549 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Mon, 7 Oct 2024 21:23:11 +0200 Subject: [PATCH 05/15] feat: Implement transformation by multiple tie points --- Cargo.lock | 66 +++ Cargo.toml | 4 +- ...apitals_model_tie_points_pixel_is_area.tif | Bin 0 -> 2482 bytes ...pitals_model_tie_points_pixel_is_point.tif | Bin 0 -> 2482 bytes src/coordinate_transform.rs | 26 +- src/coordinate_transform/tie_points.rs | 384 ++++++++++++++++++ tests/integration.rs | 16 + 7 files changed, 491 insertions(+), 5 deletions(-) create mode 100644 resources/austrian_capitals_model_tie_points_pixel_is_area.tif create mode 100644 resources/austrian_capitals_model_tie_points_pixel_is_point.tif create mode 100644 src/coordinate_transform/tie_points.rs diff --git a/Cargo.lock b/Cargo.lock index 4ef4f85..2e7e073 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -61,6 +61,12 @@ version = "2.6.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b048fb63fd8b5923fc5aa7b340d8e156aec7ec02f0c78fa8a6ddc2613f6f71de" +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + [[package]] name = "cc" version = "1.1.28" @@ -114,6 +120,15 @@ dependencies = [ "cfg-if", ] +[[package]] +name = "delaunator" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ab46e386c7a38300a0d93b0f3e484bc2ee0aded66c47b14762ec9ab383934fa" +dependencies = [ + "robust", +] + [[package]] name = "either" version = "1.13.0" @@ -166,6 +181,7 @@ checksum = "9ff16065e5720f376fbced200a5ae0f47ace85fd70b7e54269790281353b6d61" dependencies = [ "approx", "num-traits", + "rstar", "serde", ] @@ -173,10 +189,12 @@ dependencies = [ name = "geotiff" version = "0.0.2" dependencies = [ + "delaunator", "geo-types", "num-traits", "num_enum", "proj", + "rstar", "tiff", ] @@ -186,12 +204,31 @@ version = "0.3.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d2fabcfbdc87f4758337ca535fb41a6d701b65693ce38287d856d1674551ec9b" +[[package]] +name = "hash32" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "47d60b12902ba28e2730cd37e95b8c9223af2808df9e902d4df49588d1470606" +dependencies = [ + "byteorder", +] + [[package]] name = "hashbrown" version = "0.14.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" +[[package]] +name = "heapless" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0bfb9eb618601c89945a70e254898da93b13be0388091d42117462b265bb3fad" +dependencies = [ + "hash32", + "stable_deref_trait", +] + [[package]] name = "home" version = "0.5.9" @@ -455,6 +492,23 @@ version = "0.8.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" +[[package]] +name = "robust" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5864e7ef1a6b7bcf1d6ca3f655e65e724ed3b52546a0d0a663c991522f552ea" + +[[package]] +name = "rstar" +version = "0.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "133315eb94c7b1e8d0cb097e5a710d850263372fd028fff18969de708afc7008" +dependencies = [ + "heapless", + "num-traits", + "smallvec", +] + [[package]] name = "rustc-hash" version = "1.1.0" @@ -500,6 +554,18 @@ version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" +[[package]] +name = "smallvec" +version = "1.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" + +[[package]] +name = "stable_deref_trait" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3" + [[package]] name = "syn" version = "2.0.77" diff --git a/Cargo.toml b/Cargo.toml index b4eba9c..54441bf 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,9 +7,11 @@ authors = ["Dominik Bucher "] repository = "https://github.com/georust/geotiff" [dependencies] -geo-types = "0.7" +delaunator = "1.0" +geo-types = { version = "0.7", features = ["approx", "rstar_0_12"] } num_enum = "0.7" num-traits = "0.2" +rstar = "0.12" tiff = "0.9" [dev-dependencies] diff --git a/resources/austrian_capitals_model_tie_points_pixel_is_area.tif b/resources/austrian_capitals_model_tie_points_pixel_is_area.tif new file mode 100644 index 0000000000000000000000000000000000000000..3875c2cdad879fd54b46b8bc597fe4bb4fed880c GIT binary patch literal 2482 zcmebEWzb?^VBla7U}Rum2C^6#e}f1Jn_(Z2%>-q00NKnCcEcthn++-s(j>^h0;G6= z_#C4sR2*m^qZpJe0A!29)%;_WVqgVn10p>}X*k=AQ3FY@1BeZFl3hy^*lY;~26nsl z^H;@AfJ7W!Gb)f}97->MxenqG=_s5^k*Ah8^SJePpDj;=*ENW!3+j1gFaGJ`5_5Ss;3yD@?)at223+ZDIoUUTFH zR_siMe`d2udtdX~!LInQ+lztWfe&LDs0z18wl-n<2UK^^5Ty6eg0KF$-1-ZEOb}?; z&M3cufe|Qufe~zs!td$VWd18Kae`zEoJ1KI4yZGKFJfc>DFKOLPQG#nnZ&%4LFQ-+ZN0FQ*6jsO4v literal 0 HcmV?d00001 diff --git a/resources/austrian_capitals_model_tie_points_pixel_is_point.tif b/resources/austrian_capitals_model_tie_points_pixel_is_point.tif new file mode 100644 index 0000000000000000000000000000000000000000..b8b7dc2ef2c189617e0fc5945b7ed6874e8d44ee GIT binary patch literal 2482 zcmebEWzb?^VBla7U}Rum2C^6#e}f1Jn_(Z2%>-q00NKnCcEcthn++-s(j>^h0;G6= z_#C4sR2*m^qZpJe0A!29)%;_WVqgVn10p>}X*k=AQ3FY@1BeZFl3hy^*lY;~26nsl z^H;@AfJ7W!Gb)f}97->MxenqG=_s5^k*Ah8^SJePpDj;=*ENW!3+j1gFaGJ`5_5Ss;3yD@<(at223+ZDIoUUTFH zR_siMe`d2udtdX~!LInQ+lztWfe&LDs0z18wl-n<2UK^^5Ty6eg0KF$-1-ZEOb}?; z&M3cufe|Qufe~zs!td$VWd18Kae`zEoJ1KI4yZGKFJfc>DFKOLPQG#nnZ&%4LFQ-+ZN0Ffb_j{pDw literal 0 HcmV?d00001 diff --git a/src/coordinate_transform.rs b/src/coordinate_transform.rs index e56f557..33280a1 100644 --- a/src/coordinate_transform.rs +++ b/src/coordinate_transform.rs @@ -1,8 +1,12 @@ +use std::rc::Rc; use geo_types::Coord; +use rstar::RTree; use tiff::{TiffError, TiffFormatError, TiffResult}; +use crate::coordinate_transform::tie_points::{Face, TreeItem}; mod affine_transform; mod tie_point_and_pixel_scale; +mod tie_points; const MODEL_TIE_POINT_TAG: &str = "ModelTiePointTag"; const MODEL_PIXEL_SCALE_TAG: &str = "ModelPixelScaleTag"; @@ -19,7 +23,12 @@ pub enum CoordinateTransform { model_point: Coord, pixel_scale: Coord, }, - TiePoints, + TiePoints { + raster_mesh: Rc>, + raster_index: RTree, + model_mesh: Rc>, + model_index: RTree, + }, } impl CoordinateTransform { @@ -94,7 +103,7 @@ impl CoordinateTransform { Self::from_tie_point_and_pixel_scale(&tie_points, &pixel_scale) } else { - Ok(Self::TiePoints) + Self::from_tie_points(&tie_points) } } } @@ -114,7 +123,12 @@ impl CoordinateTransform { pixel_scale, coord, ), - _ => unimplemented!() + CoordinateTransform::TiePoints { + model_mesh, + raster_index, + .. + } => Self::transform_by_tie_points(raster_index, model_mesh, coord), + } } @@ -133,7 +147,11 @@ impl CoordinateTransform { pixel_scale, coord, ), - _ => unimplemented!() + CoordinateTransform::TiePoints { + model_index, + raster_mesh, + .. + } => Self::transform_by_tie_points(model_index, raster_mesh, coord), } } } diff --git a/src/coordinate_transform/tie_points.rs b/src/coordinate_transform/tie_points.rs new file mode 100644 index 0000000..2179254 --- /dev/null +++ b/src/coordinate_transform/tie_points.rs @@ -0,0 +1,384 @@ +use std::array; +use std::rc::Rc; + +use delaunator::{Point, Triangulation}; +use geo_types::Coord; +use rstar::{RTree, RTreeObject, AABB}; +use tiff::TiffResult; + +use crate::coordinate_transform::CoordinateTransform; + +impl CoordinateTransform { + pub(super) fn from_tie_points(tie_points: &[f64]) -> TiffResult { + let capacity = tie_points.iter().len() / 6; + let mut raster_points = Vec::with_capacity(capacity); + let mut model_points = Vec::with_capacity(capacity); + + for chunk in tie_points.chunks(6) { + raster_points.push(Point { + x: chunk[0], + y: chunk[1], + }); + model_points.push(Point { + x: chunk[3], + y: chunk[4], + }); + } + + let triangulation = delaunator::triangulate(&raster_points); + let raster_mesh = Rc::new(Self::build_faces(raster_points, &triangulation)); + let model_mesh = Rc::new(Self::build_faces(model_points, &triangulation)); + + Ok(Self::TiePoints { + raster_mesh: raster_mesh.clone(), + raster_index: RTree::bulk_load( + (0..raster_mesh.len()) + .map(|index| TreeItem::new(raster_mesh.clone(), index)) + .collect(), + ), + model_mesh: model_mesh.clone(), + model_index: RTree::bulk_load( + (0..model_mesh.len()) + .map(|index| TreeItem::new(model_mesh.clone(), index)) + .collect(), + ), + }) + } + + fn build_faces(points: Vec, triangulation: &Triangulation) -> Vec { + let Triangulation { + triangles, hull, .. + } = triangulation; + + let len = hull.len(); + let mut angle_bisectors = vec![None; points.len()]; + for i in 0..len { + let pi = hull[i]; + let ci = hull[(i + 1) % len]; + let ni = hull[(i + 2) % len]; + + let prev = &points[pi]; + let curr = &points[ci]; + let next = &points[ni]; + + let prev_curr = Coord { + x: curr.x - prev.x, + y: curr.y - prev.y, + } + .normalize(); + let next_curr = Coord { + x: curr.x - next.x, + y: curr.y - next.y, + } + .normalize(); + let direction = Coord { + x: prev_curr.x + next_curr.x, + y: prev_curr.y + next_curr.y, + } + .normalize(); + + angle_bisectors[ci] = Some(direction); + } + triangles + .chunks(3) + .map(|chunk| { + let i1 = chunk[0]; + let i2 = chunk[1]; + let i3 = chunk[2]; + + let b12 = hull.as_slice().contains_sequence(&chunk[0..2]); + let b23 = hull.as_slice().contains_sequence(&chunk[1..3]); + let b31 = hull.as_slice().contains_sequence(&[i3, i1]); + + let c1 = Coord { + x: points[i1].x, + y: points[i1].y, + }; + let c2 = Coord { + x: points[i2].x, + y: points[i2].y, + }; + let c3 = Coord { + x: points[i3].x, + y: points[i3].y, + }; + + let boundary = if b12 { + if b23 { + if b31 { + // Open + None + } else { + // Closed at edge 3-1 + Some(Boundary::Open { + coords: vec![c3, c1], + from_direction: angle_bisectors[i3].unwrap(), + to_direction: angle_bisectors[i1].unwrap(), + }) + } + } else if b31 { + // Closed at edge 2-3 + Some(Boundary::Open { + coords: vec![c2, c3], + from_direction: angle_bisectors[i2].unwrap(), + to_direction: angle_bisectors[i3].unwrap(), + }) + } else { + // Closed at edges 2-3 and 3-1 + Some(Boundary::Open { + coords: vec![c2, c3, c1], + from_direction: angle_bisectors[i2].unwrap(), + to_direction: angle_bisectors[i1].unwrap(), + }) + } + } else if b23 { + if b31 { + // Closed at edge 1-2 + Some(Boundary::Open { + coords: vec![c1, c2], + from_direction: angle_bisectors[i1].unwrap(), + to_direction: angle_bisectors[i2].unwrap(), + }) + } else { + // Closed at edges 1-2 and 3-1 + Some(Boundary::Open { + coords: vec![c3, c1, c2], + from_direction: angle_bisectors[i3].unwrap(), + to_direction: angle_bisectors[i2].unwrap(), + }) + } + } else if b31 { + // Closed at edges 1-2 and 2-3 + Some(Boundary::Open { + coords: vec![c1, c2, c3], + from_direction: angle_bisectors[i1].unwrap(), + to_direction: angle_bisectors[i3].unwrap(), + }) + } else { + // Closed + Some(Boundary::Closed { + coords: vec![c1, c2, c3, c1], + }) + }; + + Face { + boundary, + support_points: array::from_fn(|i| { + let point = &points[chunk[i]]; + Coord { + x: point.x, + y: point.y, + } + }), + } + }) + .collect() + } + + pub(super) fn transform_by_tie_points( + source_index: &RTree, + target_mesh: &Rc>, + coord: &Coord, + ) -> Coord { + let TreeItem { mesh, index, .. } = source_index + .locate_in_envelope_intersecting(&AABB::from_point(*coord)) + .find(|TreeItem { mesh, index, .. }| mesh[*index].contains(coord)) + .unwrap(); + let uv = mesh[*index].locate(coord); + target_mesh[*index].interpolate(uv) + } +} + +#[derive(Debug)] +pub struct Face { + boundary: Option, + support_points: [Coord; 3], +} + +impl Face { + fn contains(&self, coord: &Coord) -> bool { + let Some(boundary) = &self.boundary else { + return true; + }; + + let check = |c1: &Coord, c2: &Coord| -> bool { + ((c2.x - c1.x) * (coord.y - c1.y) - (c2.y - c1.y) * (coord.x - c1.x)).is_sign_positive() + }; + + match boundary { + Boundary::Open { + coords, + from_direction, + to_direction, + } => { + check(&(coords[0] + *from_direction), &coords[1]) + && check(&coords[1], &(coords[1] + *to_direction)) + && coords.windows(2).all(|w| check(&w[0], &w[1])) + } + Boundary::Closed { coords } => { + let len = self.support_points.len(); + (0..=len).all(|i| check(&coords[i], &coords[(i + 1) % len])) + } + } + } + + fn locate(&self, coord: &Coord) -> [f64; 2] { + let [a, b, c] = &self.support_points; + let d = c.x * (a.y - b.y) - b.x * (a.y - c.y) + a.x * (b.y - c.y); + [ + -(coord.x * (a.y - c.y) - c.x * (a.y - coord.y) + a.x * (c.y - coord.y)) / d, + (coord.x * (a.y - b.y) - b.x * (a.y - coord.y) + a.x * (b.y - coord.y)) / d, + ] + } + + fn interpolate(&self, params: [f64; 2]) -> Coord { + let [u, v] = params; + let [a, b, c] = &self.support_points; + Coord { + x: -u * a.x - v * a.x + a.x + u * b.x + v * c.x, + y: -u * a.y - v * a.y + a.y + u * b.y + v * c.y, + } + } +} + +#[derive(Debug)] +enum Boundary { + Open { + coords: Vec, + from_direction: Coord, + to_direction: Coord, + }, + Closed { + coords: Vec, + }, +} + +#[derive(Debug)] +pub struct TreeItem { + mesh: Rc>, + index: usize, + envelope: AABB, +} + +impl TreeItem { + fn new(mesh: Rc>, index: usize) -> Self { + let envelope = Self::compute_envelope(&mesh[index]); + Self { + mesh, + index, + envelope, + } + } + + fn compute_envelope(face: &Face) -> AABB { + let Some(boundary) = &face.boundary else { + return AABB::from_corners( + Coord { + x: f64::MIN, + y: f64::MIN, + }, + Coord { + x: f64::MAX, + y: f64::MAX, + }, + ); + }; + + match boundary { + Boundary::Open { + coords, + from_direction, + to_direction, + } => { + let mut lower = Coord { + x: f64::MAX, + y: f64::MAX, + }; + let mut upper = Coord { + x: f64::MIN, + y: f64::MIN, + }; + + for c in coords { + if c.x < lower.x { + lower.x = c.x; + } + if c.x > upper.x { + upper.x = c.x; + } + if c.y < lower.y { + lower.y = c.y; + } + if c.y > upper.y { + upper.y = c.y; + } + } + + for direction in [from_direction, to_direction] { + // Compute right-hand side normal vector + let nx = direction.y; + let ny = -direction.x; + + if nx.is_sign_positive() { + upper.x = f64::MAX; + } else { + lower.x = f64::MIN; + } + + if ny.is_sign_positive() { + upper.y = f64::MAX; + } else { + lower.y = f64::MIN; + } + } + + AABB::from_corners(lower, upper) + } + Boundary::Closed { coords } => AABB::from_points(coords), + } + } +} + +impl RTreeObject for TreeItem { + type Envelope = AABB; + + fn envelope(&self) -> Self::Envelope { + self.envelope + } +} + +pub(crate) trait CoordExt { + fn normalize(self) -> Self; +} + +impl CoordExt for Coord { + fn normalize(mut self) -> Self { + let len = (self.x.powi(2) + self.y.powi(2)).sqrt(); + self.x /= len; + self.y /= len; + self + } +} + +trait HullExt +where + T: PartialEq, +{ + fn contains_sequence(&self, seq: &[T]) -> bool; +} + +impl HullExt for &[T] +where + T: PartialEq, +{ + fn contains_sequence(&self, seq: &[T]) -> bool { + let len = self.len(); + for i in 0..len { + let (a, b) = seq.split_at(seq.len().min(len - i)); + if self[i..].starts_with(a) && self.starts_with(b) { + return true; + } + } + false + } +} diff --git a/tests/integration.rs b/tests/integration.rs index 27ad2e6..ef0856d 100644 --- a/tests/integration.rs +++ b/tests/integration.rs @@ -241,6 +241,22 @@ fn test_transform_by_model_transformation_pixel_is_point() { ) } +#[test] +fn test_transform_by_model_tie_points_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_tie_points_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_tie_points_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_tie_points_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + fn test_transform>(path: P, raster_type: RasterType) { let geotiff = read_geotiff(path); From 4bfe9741e47a3e3a9de5263b3ead99d2769da56c Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Mon, 7 Oct 2024 21:27:03 +0200 Subject: [PATCH 06/15] chore: Update list of authors --- Cargo.toml | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 54441bf..6595662 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,7 +3,10 @@ name = "geotiff" description = "A GeoTIFF library for Rust" version = "0.0.2" edition = "2021" -authors = ["Dominik Bucher "] +authors = [ + "Dominik Bucher ", + "Gunnar Schulze " +] repository = "https://github.com/georust/geotiff" [dependencies] From 16ddd6b21975fb0abeac2f0c3b60b1043891408e Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Tue, 8 Oct 2024 05:21:04 +0200 Subject: [PATCH 07/15] chore: Directly initialize transform and inverse transform --- src/coordinate_transform/affine_transform.rs | 35 ++++++++++---------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/coordinate_transform/affine_transform.rs b/src/coordinate_transform/affine_transform.rs index d21bd46..c97c3ba 100644 --- a/src/coordinate_transform/affine_transform.rs +++ b/src/coordinate_transform/affine_transform.rs @@ -5,13 +5,14 @@ use crate::coordinate_transform::CoordinateTransform; impl CoordinateTransform { pub fn from_transformation_matrix(transformation_matrix: [f64; 16]) -> TiffResult { - let mut transform = [0f64; 6]; - transform[0] = transformation_matrix[0]; - transform[1] = transformation_matrix[1]; - transform[2] = transformation_matrix[3]; - transform[3] = transformation_matrix[4]; - transform[4] = transformation_matrix[5]; - transform[5] = transformation_matrix[7]; + let transform = [ + transformation_matrix[0], + transformation_matrix[1], + transformation_matrix[3], + transformation_matrix[4], + transformation_matrix[5], + transformation_matrix[7], + ]; let det = transform[0] * transform[4] - transform[1] * transform[3]; if det.abs() < 0.000000000000001 { @@ -20,13 +21,14 @@ impl CoordinateTransform { ))); } - let mut inverse_transform = [0f64; 6]; - inverse_transform[0] = transform[4] / det; - inverse_transform[1] = -transform[1] / det; - inverse_transform[2] = (transform[1] * transform[5] - transform[2] * transform[4]) / det; - inverse_transform[3] = -transform[3] / det; - inverse_transform[4] = transform[0] / det; - inverse_transform[5] = (-transform[0] * transform[5] + transform[2] * transform[3]) / det; + let inverse_transform = [ + transform[4] / det, + -transform[1] / det, + (transform[1] * transform[5] - transform[2] * transform[4]) / det, + -transform[3] / det, + transform[0] / det, + (-transform[0] * transform[5] + transform[2] * transform[3]) / det, + ]; Ok(CoordinateTransform::AffineTransform { transform, @@ -34,10 +36,7 @@ impl CoordinateTransform { }) } - pub(super) fn transform_by_affine_transform( - transform: &[f64; 6], - coord: &Coord, - ) -> Coord { + pub(super) fn transform_by_affine_transform(transform: &[f64; 6], coord: &Coord) -> Coord { Coord { x: coord.x * transform[0] + coord.y * transform[1] + transform[2], y: coord.x * transform[3] + coord.y * transform[4] + transform[5], From 0e43e48d0ec6e997a1b33f279abf665777c79ffc Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Tue, 8 Oct 2024 08:49:45 +0200 Subject: [PATCH 08/15] refactor: Migrate to 'geo-index' crate --- Cargo.lock | 85 ++++++------- Cargo.toml | 4 +- src/coordinate_transform.rs | 19 +-- src/coordinate_transform/tie_points.rs | 157 ++++++++++--------------- 4 files changed, 115 insertions(+), 150 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 2e7e073..0bf02e8 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -62,10 +62,10 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b048fb63fd8b5923fc5aa7b340d8e156aec7ec02f0c78fa8a6ddc2613f6f71de" [[package]] -name = "byteorder" -version = "1.5.0" +name = "bytemuck" +version = "1.18.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" +checksum = "94bbb0ad554ad961ddc5da507a12a29b14e4ae5bda06b19f575a3e6079d2e2ae" [[package]] name = "cc" @@ -173,6 +173,25 @@ dependencies = [ "miniz_oxide", ] +[[package]] +name = "float_next_after" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8bf7cc16383c4b8d58b9905a8509f02926ce3058053c056376248d958c9df1e8" + +[[package]] +name = "geo-index" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "818dd7464a1edadbe7932a09b8e3672216597557100f23d1315f351e46c2c20e" +dependencies = [ + "bytemuck", + "float_next_after", + "num-traits", + "thiserror", + "tinyvec", +] + [[package]] name = "geo-types" version = "0.7.13" @@ -181,7 +200,6 @@ checksum = "9ff16065e5720f376fbced200a5ae0f47ace85fd70b7e54269790281353b6d61" dependencies = [ "approx", "num-traits", - "rstar", "serde", ] @@ -190,11 +208,11 @@ name = "geotiff" version = "0.0.2" dependencies = [ "delaunator", + "geo-index", "geo-types", "num-traits", "num_enum", "proj", - "rstar", "tiff", ] @@ -204,31 +222,12 @@ version = "0.3.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d2fabcfbdc87f4758337ca535fb41a6d701b65693ce38287d856d1674551ec9b" -[[package]] -name = "hash32" -version = "0.3.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "47d60b12902ba28e2730cd37e95b8c9223af2808df9e902d4df49588d1470606" -dependencies = [ - "byteorder", -] - [[package]] name = "hashbrown" version = "0.14.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" -[[package]] -name = "heapless" -version = "0.8.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0bfb9eb618601c89945a70e254898da93b13be0388091d42117462b265bb3fad" -dependencies = [ - "hash32", - "stable_deref_trait", -] - [[package]] name = "home" version = "0.5.9" @@ -498,17 +497,6 @@ version = "0.2.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e5864e7ef1a6b7bcf1d6ca3f655e65e724ed3b52546a0d0a663c991522f552ea" -[[package]] -name = "rstar" -version = "0.12.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "133315eb94c7b1e8d0cb097e5a710d850263372fd028fff18969de708afc7008" -dependencies = [ - "heapless", - "num-traits", - "smallvec", -] - [[package]] name = "rustc-hash" version = "1.1.0" @@ -554,18 +542,6 @@ version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" -[[package]] -name = "smallvec" -version = "1.13.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" - -[[package]] -name = "stable_deref_trait" -version = "1.2.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3" - [[package]] name = "syn" version = "2.0.77" @@ -619,6 +595,21 @@ dependencies = [ "weezl", ] +[[package]] +name = "tinyvec" +version = "1.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "445e881f4f6d382d5f27c034e25eb92edd7c784ceab92a0937db7f2e9471b938" +dependencies = [ + "tinyvec_macros", +] + +[[package]] +name = "tinyvec_macros" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20" + [[package]] name = "toml_datetime" version = "0.6.8" diff --git a/Cargo.toml b/Cargo.toml index 6595662..5bd7b79 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,10 +11,10 @@ repository = "https://github.com/georust/geotiff" [dependencies] delaunator = "1.0" -geo-types = { version = "0.7", features = ["approx", "rstar_0_12"] } +geo-index = "0.1" +geo-types = { version = "0.7" } num_enum = "0.7" num-traits = "0.2" -rstar = "0.12" tiff = "0.9" [dev-dependencies] diff --git a/src/coordinate_transform.rs b/src/coordinate_transform.rs index 33280a1..6046136 100644 --- a/src/coordinate_transform.rs +++ b/src/coordinate_transform.rs @@ -1,8 +1,10 @@ use std::rc::Rc; + +use geo_index::rtree::OwnedRTree; use geo_types::Coord; -use rstar::RTree; use tiff::{TiffError, TiffFormatError, TiffResult}; -use crate::coordinate_transform::tie_points::{Face, TreeItem}; + +use crate::coordinate_transform::tie_points::Face; mod affine_transform; mod tie_point_and_pixel_scale; @@ -25,9 +27,9 @@ pub enum CoordinateTransform { }, TiePoints { raster_mesh: Rc>, - raster_index: RTree, + raster_index: OwnedRTree, model_mesh: Rc>, - model_index: RTree, + model_index: OwnedRTree, }, } @@ -124,11 +126,11 @@ impl CoordinateTransform { coord, ), CoordinateTransform::TiePoints { - model_mesh, raster_index, + raster_mesh, + model_mesh, .. - } => Self::transform_by_tie_points(raster_index, model_mesh, coord), - + } => Self::transform_by_tie_points(raster_index, raster_mesh, model_mesh, coord), } } @@ -149,9 +151,10 @@ impl CoordinateTransform { ), CoordinateTransform::TiePoints { model_index, + model_mesh, raster_mesh, .. - } => Self::transform_by_tie_points(model_index, raster_mesh, coord), + } => Self::transform_by_tie_points(model_index, model_mesh, raster_mesh, coord), } } } diff --git a/src/coordinate_transform/tie_points.rs b/src/coordinate_transform/tie_points.rs index 2179254..4116a98 100644 --- a/src/coordinate_transform/tie_points.rs +++ b/src/coordinate_transform/tie_points.rs @@ -2,8 +2,9 @@ use std::array; use std::rc::Rc; use delaunator::{Point, Triangulation}; +use geo_index::rtree::sort::STRSort; +use geo_index::rtree::{OwnedRTree, RTreeBuilder, RTreeIndex}; use geo_types::Coord; -use rstar::{RTree, RTreeObject, AABB}; use tiff::TiffResult; use crate::coordinate_transform::CoordinateTransform; @@ -28,20 +29,14 @@ impl CoordinateTransform { let triangulation = delaunator::triangulate(&raster_points); let raster_mesh = Rc::new(Self::build_faces(raster_points, &triangulation)); let model_mesh = Rc::new(Self::build_faces(model_points, &triangulation)); + let raster_index = Self::build_index(&raster_mesh); + let model_index = Self::build_index(&model_mesh); Ok(Self::TiePoints { raster_mesh: raster_mesh.clone(), - raster_index: RTree::bulk_load( - (0..raster_mesh.len()) - .map(|index| TreeItem::new(raster_mesh.clone(), index)) - .collect(), - ), + raster_index, model_mesh: model_mesh.clone(), - model_index: RTree::bulk_load( - (0..model_mesh.len()) - .map(|index| TreeItem::new(model_mesh.clone(), index)) - .collect(), - ), + model_index, }) } @@ -175,17 +170,28 @@ impl CoordinateTransform { .collect() } + fn build_index(mesh: &[Face]) -> OwnedRTree { + let mut builder = RTreeBuilder::new(mesh.len()); + for face in mesh.iter() { + let (min_x, min_y, max_x, max_y) = face.compute_envelope(); + builder.add(min_x, min_y, max_x, max_y); + } + builder.finish::() + } + pub(super) fn transform_by_tie_points( - source_index: &RTree, + source_index: &OwnedRTree, + source_mesh: &Rc>, target_mesh: &Rc>, coord: &Coord, ) -> Coord { - let TreeItem { mesh, index, .. } = source_index - .locate_in_envelope_intersecting(&AABB::from_point(*coord)) - .find(|TreeItem { mesh, index, .. }| mesh[*index].contains(coord)) + let index = source_index + .search(coord.x, coord.y, coord.x, coord.y) + .into_iter() + .find(|face_index| source_mesh[*face_index].contains(coord)) .unwrap(); - let uv = mesh[*index].locate(coord); - target_mesh[*index].interpolate(uv) + let uv = source_mesh[index].locate(coord); + target_mesh[index].interpolate(uv) } } @@ -239,49 +245,32 @@ impl Face { y: -u * a.y - v * a.y + a.y + u * b.y + v * c.y, } } -} -#[derive(Debug)] -enum Boundary { - Open { - coords: Vec, - from_direction: Coord, - to_direction: Coord, - }, - Closed { - coords: Vec, - }, -} - -#[derive(Debug)] -pub struct TreeItem { - mesh: Rc>, - index: usize, - envelope: AABB, -} + fn compute_envelope(&self) -> (f64, f64, f64, f64) { + let Some(boundary) = &self.boundary else { + return (f64::MIN, f64::MIN, f64::MAX, f64::MAX); + }; -impl TreeItem { - fn new(mesh: Rc>, index: usize) -> Self { - let envelope = Self::compute_envelope(&mesh[index]); - Self { - mesh, - index, - envelope, - } - } + let mut min_x = f64::MAX; + let mut min_y = f64::MAX; + let mut max_x = f64::MIN; + let mut max_y = f64::MIN; - fn compute_envelope(face: &Face) -> AABB { - let Some(boundary) = &face.boundary else { - return AABB::from_corners( - Coord { - x: f64::MIN, - y: f64::MIN, - }, - Coord { - x: f64::MAX, - y: f64::MAX, - }, - ); + let mut update_envelope = |coords: &[Coord]| { + for c in coords { + if c.x < min_x { + min_x = c.x; + } + if c.x > min_x { + min_x = c.x; + } + if c.y < min_x { + min_x = c.y; + } + if c.y > min_x { + min_x = c.y; + } + } }; match boundary { @@ -290,29 +279,7 @@ impl TreeItem { from_direction, to_direction, } => { - let mut lower = Coord { - x: f64::MAX, - y: f64::MAX, - }; - let mut upper = Coord { - x: f64::MIN, - y: f64::MIN, - }; - - for c in coords { - if c.x < lower.x { - lower.x = c.x; - } - if c.x > upper.x { - upper.x = c.x; - } - if c.y < lower.y { - lower.y = c.y; - } - if c.y > upper.y { - upper.y = c.y; - } - } + update_envelope(coords); for direction in [from_direction, to_direction] { // Compute right-hand side normal vector @@ -320,31 +287,35 @@ impl TreeItem { let ny = -direction.x; if nx.is_sign_positive() { - upper.x = f64::MAX; + max_x = f64::MAX; } else { - lower.x = f64::MIN; + min_x = f64::MIN; } if ny.is_sign_positive() { - upper.y = f64::MAX; + max_y = f64::MAX; } else { - lower.y = f64::MIN; + min_y = f64::MIN; } } - - AABB::from_corners(lower, upper) } - Boundary::Closed { coords } => AABB::from_points(coords), + Boundary::Closed { coords } => update_envelope(coords), } + + (min_x, min_y, max_x, max_y) } } -impl RTreeObject for TreeItem { - type Envelope = AABB; - - fn envelope(&self) -> Self::Envelope { - self.envelope - } +#[derive(Debug)] +enum Boundary { + Open { + coords: Vec, + from_direction: Coord, + to_direction: Coord, + }, + Closed { + coords: Vec, + }, } pub(crate) trait CoordExt { From f1f8532e53b3b45d1f535631c7e5609c9c15993b Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Wed, 9 Oct 2024 09:58:49 +0200 Subject: [PATCH 09/15] refactor: Move coordinate transformation tests to separate file --- tests/common/mod.rs | 8 ++ tests/integration.rs | 216 +------------------------------------------ tests/transform.rs | 213 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 224 insertions(+), 213 deletions(-) create mode 100644 tests/common/mod.rs create mode 100644 tests/transform.rs diff --git a/tests/common/mod.rs b/tests/common/mod.rs new file mode 100644 index 0000000..f6774ee --- /dev/null +++ b/tests/common/mod.rs @@ -0,0 +1,8 @@ +use std::fs::File; +use std::path::Path; + +use geotiff::GeoTiff; + +pub fn read_geotiff>(path: P) -> GeoTiff { + GeoTiff::read(File::open(path).expect("File I/O error")).expect("File I/O error") +} diff --git a/tests/integration.rs b/tests/integration.rs index ef0856d..ba06fae 100644 --- a/tests/integration.rs +++ b/tests/integration.rs @@ -1,72 +1,8 @@ -use std::fs::File; -use std::path::Path; - +use common::read_geotiff; use geo_types::{Coord, Rect}; -use geotiff::{GeoKeyDirectory, GeoTiff, RasterType}; -use proj::Proj; - -trait RectExt { - fn rounded(self) -> Self; -} +use geotiff::{GeoKeyDirectory, RasterType}; -impl RectExt for Rect { - fn rounded(mut self) -> Self { - let factor = 10f64.powi(8); - let mut min = self.min(); - let mut max = self.max(); - min.x = (min.x * factor).round() / factor; - min.y = (min.y * factor).round() / factor; - max.x = (max.x * factor).round() / factor; - max.y = (max.y * factor).round() / factor; - self.set_min(min); - self.set_max(max); - self - } -} - -const BREGENZ: Coord = Coord { - x: 9.74926, - y: 47.50315, -}; -const EISENSTADT: Coord = Coord { - x: 15.43301, - y: 47.06298, -}; -const GRAZ: Coord = Coord { - x: 15.43301, - y: 47.06298, -}; -const INNSBRUCK: Coord = Coord { - x: 11.39960, - y: 47.26239, -}; -const KLAGENFURT: Coord = Coord { - x: 14.31528, - y: 46.62366, -}; -const LINZ: Coord = Coord { - x: 14.30571, - y: 48.27532, -}; -const SALZBURG: Coord = Coord { - x: 13.05345, - y: 47.80763, -}; -const SANKT_POELTEN: Coord = Coord { - x: 15.62291, - y: 48.20440, -}; -const VIENNA: Coord = Coord { - x: 16.37499, - y: 48.22158, -}; - -const WHITE: u8 = 255; -const BLACK: u8 = 0; - -fn read_geotiff>(path: P) -> GeoTiff { - GeoTiff::read(File::open(path).expect("File I/O error")).expect("File I/O error") -} +mod common; #[test] fn test_load_marbles() { @@ -208,149 +144,3 @@ fn test_load_merc() { ) ); } - -#[test] -fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_area() { - test_transform( - "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif", - RasterType::RasterPixelIsArea, - ) -} - -#[test] -fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_point() { - test_transform( - "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif", - RasterType::RasterPixelIsPoint, - ) -} - -#[test] -fn test_transform_by_model_transformation_pixel_is_area() { - test_transform( - "resources/austrian_capitals_model_transformation_pixel_is_area.tif", - RasterType::RasterPixelIsArea, - ) -} - -#[test] -fn test_transform_by_model_transformation_pixel_is_point() { - test_transform( - "resources/austrian_capitals_model_transformation_pixel_is_point.tif", - RasterType::RasterPixelIsPoint, - ) -} - -#[test] -fn test_transform_by_model_tie_points_pixel_is_area() { - test_transform( - "resources/austrian_capitals_model_tie_points_pixel_is_area.tif", - RasterType::RasterPixelIsArea, - ) -} - -#[test] -fn test_transform_by_model_tie_points_pixel_is_point() { - test_transform( - "resources/austrian_capitals_model_tie_points_pixel_is_point.tif", - RasterType::RasterPixelIsPoint, - ) -} - -fn test_transform>(path: P, raster_type: RasterType) { - let geotiff = read_geotiff(path); - - let proj = Proj::new_known_crs( - "EPSG:4326", - &format!("EPSG:{}", geotiff.geo_key_directory.projected_type.unwrap()), - None, - ) - .unwrap(); - - let mut bregenz = proj.project(BREGENZ, false).unwrap(); - let mut eisenstadt = proj.project(EISENSTADT, false).unwrap(); - let mut graz = proj.project(GRAZ, false).unwrap(); - let mut innsbruck = proj.project(INNSBRUCK, false).unwrap(); - let mut klagenfurt = proj.project(KLAGENFURT, false).unwrap(); - let mut linz = proj.project(LINZ, false).unwrap(); - let mut salzburg = proj.project(SALZBURG, false).unwrap(); - let mut sankt_poelten = proj.project(SANKT_POELTEN, false).unwrap(); - let mut vienna = proj.project(VIENNA, false).unwrap(); - let mut expected_model_extent = Rect::new( - Coord { - x: 4302000.0, - y: 2621000.0, - }, - Coord { - x: 4809000.0, - y: 2811000.0, - }, - ); - - // If raster type is RasterPixelIsPoint, shift coordinates accordingly - if raster_type == RasterType::RasterPixelIsPoint { - let mut min = expected_model_extent.min(); - let mut max = expected_model_extent.max(); - - let mut coords = [ - &mut bregenz, - &mut eisenstadt, - &mut graz, - &mut innsbruck, - &mut klagenfurt, - &mut linz, - &mut salzburg, - &mut sankt_poelten, - &mut vienna, - &mut min, - &mut max, - ]; - - for coord in coords.iter_mut() { - coord.x -= 500.0; - coord.y += 500.0; - } - - expected_model_extent.set_min(min); - expected_model_extent.set_max(max); - } - - assert_eq!(geotiff.model_extent().rounded(), expected_model_extent); - - // Non-capital location should return background color - assert_eq!( - geotiff.get_value_at::(&expected_model_extent.center(), 0), - Some(WHITE) - ); - - // A location outside the model extent should return None - - let min = expected_model_extent.min(); - let max = expected_model_extent.max(); - - assert_eq!( - geotiff.get_value_at::(&Coord { x: min.x, y: min.y }, 0), - None - ); - assert_eq!( - geotiff.get_value_at::( - &Coord { - x: max.x + 1.0, - y: max.y + 1.0, - }, - 0 - ), - None - ); - - // Each capital location should return Some(BLACK) - assert_eq!(geotiff.get_value_at::(&bregenz, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&eisenstadt, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&graz, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&innsbruck, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&klagenfurt, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&linz, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&salzburg, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&sankt_poelten, 0), Some(BLACK)); - assert_eq!(geotiff.get_value_at::(&vienna, 0), Some(BLACK)); -} diff --git a/tests/transform.rs b/tests/transform.rs new file mode 100644 index 0000000..735ae58 --- /dev/null +++ b/tests/transform.rs @@ -0,0 +1,213 @@ +use std::path::Path; + +use common::read_geotiff; +use geo_types::{Coord, Rect}; +use geotiff::RasterType; +use proj::Proj; + +mod common; + +trait RectExt { + fn rounded(self) -> Self; +} + +impl RectExt for Rect { + fn rounded(mut self) -> Self { + let factor = 10f64.powi(8); + let mut min = self.min(); + let mut max = self.max(); + min.x = (min.x * factor).round() / factor; + min.y = (min.y * factor).round() / factor; + max.x = (max.x * factor).round() / factor; + max.y = (max.y * factor).round() / factor; + self.set_min(min); + self.set_max(max); + self + } +} + +const BREGENZ: Coord = Coord { + x: 9.74926, + y: 47.50315, +}; +const EISENSTADT: Coord = Coord { + x: 15.43301, + y: 47.06298, +}; +const GRAZ: Coord = Coord { + x: 15.43301, + y: 47.06298, +}; +const INNSBRUCK: Coord = Coord { + x: 11.39960, + y: 47.26239, +}; +const KLAGENFURT: Coord = Coord { + x: 14.31528, + y: 46.62366, +}; +const LINZ: Coord = Coord { + x: 14.30571, + y: 48.27532, +}; +const SALZBURG: Coord = Coord { + x: 13.05345, + y: 47.80763, +}; +const SANKT_POELTEN: Coord = Coord { + x: 15.62291, + y: 48.20440, +}; +const VIENNA: Coord = Coord { + x: 16.37499, + y: 48.22158, +}; + +const WHITE: u8 = 255; +const BLACK: u8 = 0; + +#[test] +fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + +#[test] +fn test_transform_by_model_transformation_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_transformation_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_transformation_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_transformation_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + +#[test] +fn test_transform_by_model_tie_points_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_tie_points_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_tie_points_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_tie_points_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + +fn test_transform>(path: P, raster_type: RasterType) { + let geotiff = read_geotiff(path); + + let proj = Proj::new_known_crs( + "EPSG:4326", + &format!("EPSG:{}", geotiff.geo_key_directory.projected_type.unwrap()), + None, + ) + .unwrap(); + + let mut bregenz = proj.project(BREGENZ, false).unwrap(); + let mut eisenstadt = proj.project(EISENSTADT, false).unwrap(); + let mut graz = proj.project(GRAZ, false).unwrap(); + let mut innsbruck = proj.project(INNSBRUCK, false).unwrap(); + let mut klagenfurt = proj.project(KLAGENFURT, false).unwrap(); + let mut linz = proj.project(LINZ, false).unwrap(); + let mut salzburg = proj.project(SALZBURG, false).unwrap(); + let mut sankt_poelten = proj.project(SANKT_POELTEN, false).unwrap(); + let mut vienna = proj.project(VIENNA, false).unwrap(); + let mut expected_model_extent = Rect::new( + Coord { + x: 4302000.0, + y: 2621000.0, + }, + Coord { + x: 4809000.0, + y: 2811000.0, + }, + ); + + // If raster type is RasterPixelIsPoint, shift coordinates accordingly + if raster_type == RasterType::RasterPixelIsPoint { + let mut min = expected_model_extent.min(); + let mut max = expected_model_extent.max(); + + let mut coords = [ + &mut bregenz, + &mut eisenstadt, + &mut graz, + &mut innsbruck, + &mut klagenfurt, + &mut linz, + &mut salzburg, + &mut sankt_poelten, + &mut vienna, + &mut min, + &mut max, + ]; + + for coord in coords.iter_mut() { + coord.x -= 500.0; + coord.y += 500.0; + } + + expected_model_extent.set_min(min); + expected_model_extent.set_max(max); + } + + assert_eq!(geotiff.model_extent().rounded(), expected_model_extent); + + // Non-capital location should return background color + assert_eq!( + geotiff.get_value_at::(&expected_model_extent.center(), 0), + Some(WHITE) + ); + + // A location outside the model extent should return None + + let min = expected_model_extent.min(); + let max = expected_model_extent.max(); + + assert_eq!( + geotiff.get_value_at::(&Coord { x: min.x, y: min.y }, 0), + None + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: max.x + 1.0, + y: max.y + 1.0, + }, + 0 + ), + None + ); + + // Each capital location should return Some(BLACK) + assert_eq!(geotiff.get_value_at::(&bregenz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&eisenstadt, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&graz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&innsbruck, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&klagenfurt, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&linz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&salzburg, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&sankt_poelten, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&vienna, 0), Some(BLACK)); +} From 385163845f860e62898df6f39fa60c492510a39f Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Wed, 9 Oct 2024 13:31:33 +0200 Subject: [PATCH 10/15] build: Make tie points transformation configurable via feature flag --- Cargo.toml | 7 +++++-- src/coordinate_transform.rs | 18 +++++++++++++++++- tests/transform.rs | 2 ++ 3 files changed, 24 insertions(+), 3 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 5bd7b79..866af36 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -10,8 +10,8 @@ authors = [ repository = "https://github.com/georust/geotiff" [dependencies] -delaunator = "1.0" -geo-index = "0.1" +delaunator = { version = "1.0", optional = true } +geo-index = { version = "0.1", optional = true } geo-types = { version = "0.7" } num_enum = "0.7" num-traits = "0.2" @@ -19,3 +19,6 @@ tiff = "0.9" [dev-dependencies] proj = "0.27" + +[features] +tie-points = ["dep:delaunator", "dep:geo-index"] diff --git a/src/coordinate_transform.rs b/src/coordinate_transform.rs index 6046136..17dd25a 100644 --- a/src/coordinate_transform.rs +++ b/src/coordinate_transform.rs @@ -1,13 +1,17 @@ +#[cfg(feature = "tie-points")] use std::rc::Rc; +#[cfg(feature = "tie-points")] use geo_index::rtree::OwnedRTree; use geo_types::Coord; use tiff::{TiffError, TiffFormatError, TiffResult}; +#[cfg(feature = "tie-points")] use crate::coordinate_transform::tie_points::Face; mod affine_transform; mod tie_point_and_pixel_scale; +#[cfg(feature = "tie-points")] mod tie_points; const MODEL_TIE_POINT_TAG: &str = "ModelTiePointTag"; @@ -25,6 +29,7 @@ pub enum CoordinateTransform { model_point: Coord, pixel_scale: Coord, }, + #[cfg(feature = "tie-points")] TiePoints { raster_mesh: Rc>, raster_index: OwnedRTree, @@ -105,7 +110,16 @@ impl CoordinateTransform { Self::from_tie_point_and_pixel_scale(&tie_points, &pixel_scale) } else { - Self::from_tie_points(&tie_points) + #[cfg(feature = "tie-points")] + { + Self::from_tie_points(&tie_points) + } + #[cfg(not(feature = "tie-points"))] + { + Err(TiffError::FormatError(TiffFormatError::Format( + "Transformation by tie points is not supported".into(), + ))) + } } } } @@ -125,6 +139,7 @@ impl CoordinateTransform { pixel_scale, coord, ), + #[cfg(feature = "tie-points")] CoordinateTransform::TiePoints { raster_index, raster_mesh, @@ -149,6 +164,7 @@ impl CoordinateTransform { pixel_scale, coord, ), + #[cfg(feature = "tie-points")] CoordinateTransform::TiePoints { model_index, model_mesh, diff --git a/tests/transform.rs b/tests/transform.rs index 735ae58..2634933 100644 --- a/tests/transform.rs +++ b/tests/transform.rs @@ -98,6 +98,7 @@ fn test_transform_by_model_transformation_pixel_is_point() { ) } +#[cfg(feature = "tie-points")] #[test] fn test_transform_by_model_tie_points_pixel_is_area() { test_transform( @@ -106,6 +107,7 @@ fn test_transform_by_model_tie_points_pixel_is_area() { ) } +#[cfg(feature = "tie-points")] #[test] fn test_transform_by_model_tie_points_pixel_is_point() { test_transform( From 64a924545cf6d9ae20339e297b4ea16abdcdcc88 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Tue, 19 Nov 2024 11:34:46 +0100 Subject: [PATCH 11/15] feat: Add 'Undefined' and 'UserDefined' raster type enum variants --- src/geo_key_directory.rs | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/geo_key_directory.rs b/src/geo_key_directory.rs index f23b6a8..f0f890e 100644 --- a/src/geo_key_directory.rs +++ b/src/geo_key_directory.rs @@ -101,11 +101,12 @@ impl GeoKeyDirectory { GeoKeyDirectoryTag::RasterType => { let raster_type = Self::get_short(key_tag, location_tag, *count, *value_or_offset)?; - directory.raster_type = Some(RasterType::try_from(raster_type).map_err(|_| { - TiffError::FormatError(TiffFormatError::Format(format!( - "Unknown raster type: {raster_type}" - ))) - })?) + directory.raster_type = + Some(RasterType::try_from(raster_type).map_err(|_| { + TiffError::FormatError(TiffFormatError::Format(format!( + "Unknown raster type: {raster_type}" + ))) + })?) } GeoKeyDirectoryTag::Citation => { directory.citation = Self::get_string( @@ -663,9 +664,14 @@ enum GeoKeyDirectoryTag { VerticalUnits = 4099, } +/// The raster type establishes the raster space used. +/// +/// Ref: https://docs.ogc.org/is/19-008r4/19-008r4.html#_requirements_class_gtrastertypegeokey #[derive(Debug, Clone, Copy, PartialEq, TryFromPrimitive, IntoPrimitive)] #[repr(u16)] pub enum RasterType { + Undefined = 0, RasterPixelIsArea = 1, RasterPixelIsPoint = 2, + UserDefined = 32767, } From e732c38d16199ad4deaefa0aea58061ff1210014 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Tue, 19 Nov 2024 11:36:31 +0100 Subject: [PATCH 12/15] doc: Document 'CoordinateTransform' enum --- src/coordinate_transform.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/coordinate_transform.rs b/src/coordinate_transform.rs index 17dd25a..a6b0c20 100644 --- a/src/coordinate_transform.rs +++ b/src/coordinate_transform.rs @@ -18,6 +18,9 @@ const MODEL_TIE_POINT_TAG: &str = "ModelTiePointTag"; const MODEL_PIXEL_SCALE_TAG: &str = "ModelPixelScaleTag"; const MODEL_TRANSFORMATION_TAG: &str = "ModelTransformationTag"; +/// Defines the transformation between raster space and model space. +/// +/// Ref: https://docs.ogc.org/is/19-008r4/19-008r4.html#_raster_to_model_coordinate_transformation_requirements #[derive(Debug)] pub enum CoordinateTransform { AffineTransform { From 930fac458e58b3fbd5549bf917461165d2afa4d4 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Wed, 20 Nov 2024 11:58:56 +0100 Subject: [PATCH 13/15] doc: Document public methods in GeoTiff 'struct' --- src/lib.rs | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index 438cb22..c4f4566 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,6 +1,7 @@ //! A [GeoTIFF](https://www.ogc.org/standard/geotiff) library for Rust use std::any::type_name; use std::io::{Read, Seek}; + use geo_types::{Coord, Rect}; use num_traits::FromPrimitive; use tiff::decoder::{Decoder, DecodingResult}; @@ -46,6 +47,7 @@ pub struct GeoTiff { } impl GeoTiff { + /// Reads a GeoTIFF from the given source. pub fn read(reader: R) -> TiffResult { let mut decoder = Decoder::new(reader)?; @@ -83,6 +85,7 @@ impl GeoTiff { }) } + /// Returns the extent of the image in model space. pub fn model_extent(&self) -> Rect { let offset = self.raster_offset(); let lower = Coord { @@ -104,6 +107,8 @@ impl GeoTiff { } } + /// Returns the value at the given location for the specified sample. + /// The coordinates are in model space. pub fn get_value_at( &self, coord: &Coord, From a592a1978bcb828cb1f48080a3a546080638fe88 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Wed, 20 Nov 2024 12:13:44 +0100 Subject: [PATCH 14/15] doc: Add reference to standard for raster offset calculation --- src/lib.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/lib.rs b/src/lib.rs index c4f4566..f69276d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -151,6 +151,7 @@ impl GeoTiff { Some(transform) => transform.transform_to_raster(coord), }; + // See https://docs.ogc.org/is/19-008r4/19-008r4.html#_raster_space for reference let raster_offset = self.raster_offset(); coord.x -= raster_offset; coord.y -= raster_offset; From f13418ea8f2eb092b8b5375448c5f1277494b9e2 Mon Sep 17 00:00:00 2001 From: Gunnar Schulze Date: Wed, 20 Nov 2024 12:27:35 +0100 Subject: [PATCH 15/15] fix: Fix faulty envelope computation --- src/coordinate_transform/tie_points.rs | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/coordinate_transform/tie_points.rs b/src/coordinate_transform/tie_points.rs index 4116a98..4498170 100644 --- a/src/coordinate_transform/tie_points.rs +++ b/src/coordinate_transform/tie_points.rs @@ -258,18 +258,10 @@ impl Face { let mut update_envelope = |coords: &[Coord]| { for c in coords { - if c.x < min_x { - min_x = c.x; - } - if c.x > min_x { - min_x = c.x; - } - if c.y < min_x { - min_x = c.y; - } - if c.y > min_x { - min_x = c.y; - } + min_x = min_x.min(c.x); + min_y = min_y.min(c.y); + max_x = max_x.max(c.x); + max_y = max_y.max(c.y); } };