diff --git a/Cargo.lock b/Cargo.lock index 00a4b16d6..6cb773b7b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,22 @@ # It is not intended for manual editing. version = 4 +[[package]] +name = "ab_glyph" +version = "0.2.32" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "01c0457472c38ea5bd1c3b5ada5e368271cb550be7a4ca4a0b4634e9913f6cc2" +dependencies = [ + "ab_glyph_rasterizer", + "owned_ttf_parser", +] + +[[package]] +name = "ab_glyph_rasterizer" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "366ffbaa4442f4684d91e2cd7c5ea7c4ed8add41959a31447066e279e432b618" + [[package]] name = "actix-codec" version = "0.5.2" @@ -2487,7 +2503,7 @@ dependencies = [ "memmap2", "slotmap", "tinyvec", - "ttf-parser", + "ttf-parser 0.24.1", ] [[package]] @@ -3307,6 +3323,24 @@ dependencies = [ "quick-error", ] +[[package]] +name = "imageproc" +version = "0.25.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2393fb7808960751a52e8a154f67e7dd3f8a2ef9bd80d1553078a7b4e8ed3f0d" +dependencies = [ + "ab_glyph", + "approx", + "getrandom 0.2.16", + "image", + "itertools 0.12.1", + "nalgebra", + "num 0.4.3", + "rand 0.8.5", + "rand_distr", + "rayon", +] + [[package]] name = "imagesize" version = "0.13.0" @@ -3406,6 +3440,15 @@ version = "1.70.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a6cb138bb79a146c1bd460005623e142ef0181e3d0219cb493e02f7d08a35695" +[[package]] +name = "itertools" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba291022dbbd398a455acf126c1e341954079855bc60dfdda641363bd6922569" +dependencies = [ + "either", +] + [[package]] name = "itertools" version = "0.13.0" @@ -3818,6 +3861,7 @@ dependencies = [ "env_logger", "futures", "image", + "imageproc", "insta", "itertools 0.14.0", "maplibre_native", @@ -3877,6 +3921,16 @@ version = "0.8.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "47e1ffaa40ddd1f3ed91f717a33c8c0ee23fff369e3aa8772b9605cc1d22f4c3" +[[package]] +name = "matrixmultiply" +version = "0.3.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08" +dependencies = [ + "autocfg", + "rawpointer", +] + [[package]] name = "maybe-async" version = "0.2.10" @@ -4041,6 +4095,21 @@ dependencies = [ "serde", ] +[[package]] +name = "nalgebra" +version = "0.32.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7b5c17de023a86f59ed79891b2e5d5a94c705dbe904a5b5c9c952ea6221b03e4" +dependencies = [ + "approx", + "matrixmultiply", + "num-complex 0.4.6", + "num-rational 0.4.2", + "num-traits", + "simba", + "typenum", +] + [[package]] name = "native-tls" version = "0.2.14" @@ -4389,6 +4458,15 @@ version = "0.5.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1a80800c0488c3a21695ea981a54918fbb37abf04f4d0720c453632255e2ff0e" +[[package]] +name = "owned_ttf_parser" +version = "0.25.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "36820e9051aca1014ddc75770aab4d68bc1e9e632f0f5627c4086bc216fb583b" +dependencies = [ + "ttf-parser 0.25.1", +] + [[package]] name = "oxipng" version = "9.1.5" @@ -5264,6 +5342,16 @@ dependencies = [ "getrandom 0.3.4", ] +[[package]] +name = "rand_distr" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31" +dependencies = [ + "num-traits", + "rand 0.8.5", +] + [[package]] name = "rav1e" version = "0.8.1" @@ -5314,6 +5402,12 @@ dependencies = [ "rgb", ] +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + [[package]] name = "rayon" version = "1.11.0" @@ -5791,7 +5885,7 @@ dependencies = [ "core_maths", "log", "smallvec", - "ttf-parser", + "ttf-parser 0.24.1", "unicode-bidi-mirroring", "unicode-ccc", "unicode-properties", @@ -5804,6 +5898,15 @@ version = "1.0.22" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a50f4cf475b65d88e057964e0e9bb1f0aa9bbb2036dc65c64596b42932536984" +[[package]] +name = "safe_arch" +version = "0.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "96b02de82ddbe1b636e6170c21be622223aea188ef2e139be0a5b219ec215323" +dependencies = [ + "bytemuck", +] + [[package]] name = "same-file" version = "1.0.6" @@ -6132,6 +6235,19 @@ dependencies = [ "rand_core 0.6.4", ] +[[package]] +name = "simba" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "061507c94fc6ab4ba1c9a0305018408e312e17c041eb63bef8aa726fa33aceae" +dependencies = [ + "approx", + "num-complex 0.4.6", + "num-traits", + "paste", + "wide", +] + [[package]] name = "simd-adler32" version = "0.3.8" @@ -7323,6 +7439,12 @@ dependencies = [ "core_maths", ] +[[package]] +name = "ttf-parser" +version = "0.25.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2df906b07856748fa3f6e0ad0cbaa047052d4a7dd609e231c4f72cee8c36f31" + [[package]] name = "twox-hash" version = "2.1.2" @@ -7728,6 +7850,16 @@ dependencies = [ "web-sys", ] +[[package]] +name = "wide" +version = "0.7.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03" +dependencies = [ + "bytemuck", + "safe_arch", +] + [[package]] name = "winapi" version = "0.3.9" diff --git a/Cargo.toml b/Cargo.toml index d01bcf798..c8cdb19a1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -56,6 +56,7 @@ flate2 = "1" flume = "0.12" futures = "0.3" image = "0.25.8" +imageproc = "0.25.0" indoc = "2" insta = "1.44.3" itertools = "0.14" diff --git a/martin-core/Cargo.toml b/martin-core/Cargo.toml index 16b5629eb..2e08fe4af 100644 --- a/martin-core/Cargo.toml +++ b/martin-core/Cargo.toml @@ -33,7 +33,7 @@ postgres = [ "dep:serde_json", "_tiles", ] -unstable-cog = ["dep:png", "dep:tiff", "_tiles"] +unstable-cog = ["dep:png", "dep:tiff", "_tiles", "dep:image", "dep:imageproc"] unstable-rendering = ["styles", "dep:maplibre_native", "dep:image"] fonts = [ "dep:bit-set", @@ -67,6 +67,8 @@ bit-set = { workspace = true, optional = true } dashmap = { workspace = true, optional = true } deadpool-postgres = { workspace = true, optional = true } futures = { workspace = true, optional = true } +image = { workspace = true, optional = true } +imageproc = { workspace = true, optional = true } itertools = { workspace = true, optional = true } martin-tile-utils.workspace = true mbtiles = { workspace = true, optional = true } diff --git a/martin-core/src/tiles/cog/errors.rs b/martin-core/src/tiles/cog/errors.rs index e1e327258..d5423c114 100644 --- a/martin-core/src/tiles/cog/errors.rs +++ b/martin-core/src/tiles/cog/errors.rs @@ -83,6 +83,10 @@ pub enum CogError { )] GetFullResolutionFailed(PathBuf), + /// Failed to create image buffer + #[error("Failed to create image buffer for {0}: {1}")] + ImageBufferCreationFailed(PathBuf, String), + /// IO error. #[error("IO error {0}: {1}")] IoError(#[source] std::io::Error, PathBuf), diff --git a/martin-core/src/tiles/cog/image.rs b/martin-core/src/tiles/cog/image.rs index 7e4fef885..68559bec3 100644 --- a/martin-core/src/tiles/cog/image.rs +++ b/martin-core/src/tiles/cog/image.rs @@ -2,11 +2,13 @@ use std::fs::File; use std::io::BufWriter; use std::path::Path; +use image::{ImageBuffer, Rgba}; use martin_tile_utils::{TileCoord, TileData}; +use tiff::ColorType; use tiff::decoder::{Decoder, DecodingResult}; +use super::CogError; use crate::tiles::MartinCoreResult; -use crate::tiles::cog::CogError; /// Image represents a single image in a COG file. A tiff file may contain many images. /// This struct contains information and methods for taking tiles from the image. @@ -14,22 +16,157 @@ use crate::tiles::cog::CogError; pub struct Image { /// The Image File Directory index represents IDF entry with the image pointers to the actual image data. ifd_index: usize, + /// The extent of the image in model units, represented as [`min_x`, `min_y`, `max_x`, `max_y`]. + extent: [f64; 4], + /// The origin of the image in model units. + origin: [f64; 3], /// Number of tiles in a row of this image tiles_across: u32, - /// Number of tiles in a column of this image + /// Number of tiles in a column of this image tiles_down: u32, + /// Tile size in pixels + tile_size: (u32, u32), + /// Resolution of the image in model units per pixel + resolution: (f64, f64), } impl Image { - /// Creates a new image with the specified IFD index and tile dimensions. - pub fn new(ifd_index: usize, tiles_across: u32, tiles_down: u32) -> Self { + pub fn new( + ifd_index: usize, + extent: [f64; 4], + origin: [f64; 3], + tiles_across: u32, + tiles_down: u32, + tile_size: (u32, u32), + resolution: (f64, f64), + ) -> Self { Self { ifd_index, + extent, + origin, tiles_across, tiles_down, + tile_size, + resolution, } } + #[allow(clippy::cast_possible_truncation)] + pub fn get_tile_webmercator( + &self, + decoder: &mut Decoder, + xyz: TileCoord, + nodata: Option, + path: &Path, + ) -> Result { + let bbox = martin_tile_utils::xyz_to_bbox_webmercator(xyz.z, xyz.x, xyz.y, xyz.x, xyz.y); + #[allow(clippy::cast_sign_loss)] + let nodata_u8 = nodata.map(|v| v as u8); + let bytes = self.clip(decoder, bbox, 256, nodata_u8, path)?; + Ok(bytes) + } + + /// Clips the image to the specified bounding box and returns the PNG data. + #[allow(clippy::cast_sign_loss)] + #[allow(clippy::cast_possible_truncation)] + fn clip( + &self, + decoder: &mut Decoder, + bbox: [f64; 4], + output_size: u32, + nodata: Option, + path: &Path, + ) -> Result { + decoder + .seek_to_image(self.ifd_index()) + .map_err(|e| CogError::IfdSeekFailed(e, self.ifd_index(), path.to_path_buf()))?; + + let target_w = ((bbox[2] - bbox[0]) / self.resolution.0).round() as u32; + let target_h = ((bbox[3] - bbox[1]) / self.resolution.1.abs()).round() as u32; + let mut target = vec![0; (target_w * target_h * 4) as usize]; + + // draw each tile on the target + let intersected_tiles = tiles_intersected( + self.tile_size, + self.resolution, + self.extent, + (self.tiles_across, self.tiles_down), + bbox, + ); + for (col, row) in intersected_tiles { + let Some(idx) = self.get_tile_index(TileCoord { + z: 0, // actually this z is not used, so we can use 0 here + x: col, + y: row, + }) else { + continue; + }; + + let origin_x = self.origin[0]; + let origin_y = self.origin[1]; + let tile_min_x = origin_x + f64::from(col * self.tile_size.0) * self.resolution.0; + let tile_max_y = origin_y - f64::from(row * self.tile_size.1) * self.resolution.1.abs(); + + let geo_offset_x = tile_min_x - bbox[0]; + let geo_offset_y = bbox[3] - tile_max_y; // Use window's max Y + + let offset_x = (geo_offset_x / self.resolution.0).round() as i64; + let offset_y = (geo_offset_y / self.resolution.1.abs()).round() as i64; + + let color_type = decoder + .colortype() + .map_err(|e| CogError::InvalidTiffFile(e, path.to_path_buf()))?; //FIXME: maybe make color_type as prop of Image Struct? + let components_count = match color_type { + ColorType::RGB(_) => 3, + ColorType::RGBA(_) => 4, + ct => { + return Err(CogError::NotSupportedColorTypeAndBitDepth( + ct, + path.to_path_buf(), + ))?; + } + }; + + let (tile_w, tile_h) = decoder.chunk_data_dimensions(idx); + let tile_data = decoder.read_chunk(idx).map_err(|e| { + CogError::ReadChunkFailed(e, idx, self.ifd_index(), path.to_path_buf()) + })?; + match (tile_data, color_type) { + (DecodingResult::U8(vec), ColorType::RGB(_) | ColorType::RGBA(_)) => draw_tile( + &vec, + components_count, + nodata, + (tile_w, tile_h), + (target_w, target_h), + (offset_x, offset_y), + &mut target, + ), + (_, _) => { + return Err(CogError::NotSupportedColorTypeAndBitDepth( + color_type, + path.to_path_buf(), + )); + } //todo do others in next PRs, a lot of discussion would be needed + } + } + + let result_image: ImageBuffer, Vec> = + ImageBuffer::from_raw(target_w, target_h, target).ok_or_else(|| { + CogError::ImageBufferCreationFailed( + path.to_path_buf(), + format!("Failed to create image buffer with dimensions {target_w}x{target_h}"), + ) + })?; + let resized = image::imageops::resize( + &result_image, + output_size, + output_size, + image::imageops::FilterType::Nearest, //FIXME should make this configurable + ); + let png = encode_rgba_as_png(output_size, output_size, resized.as_raw(), path)?; + Ok(png) + } + /// Retrieves a tile from the image, decodes it, and converts it to PNG format. #[expect(clippy::cast_sign_loss, clippy::cast_possible_truncation)] pub fn get_tile( @@ -61,7 +198,7 @@ impl Image { // FIXME: do more research on the not u8 case, is this the right way to do it? let png_file_bytes = match (decode_result, color_type) { - (DecodingResult::U8(vec), tiff::ColorType::RGB(_)) => rgb_to_png( + (DecodingResult::U8(vec), ColorType::RGB(_)) => rgb_to_png( vec, (tile_width, tile_height), (data_width, data_height), @@ -69,7 +206,7 @@ impl Image { nodata.map(|v| v as u8), path, ), - (DecodingResult::U8(vec), tiff::ColorType::RGBA(_)) => rgb_to_png( + (DecodingResult::U8(vec), ColorType::RGBA(_)) => rgb_to_png( vec, (tile_width, tile_height), (data_width, data_height), @@ -91,6 +228,14 @@ impl Image { self.ifd_index } + pub fn resolution(&self) -> (f64, f64) { + self.resolution + } + + pub fn tile_size(&self) -> (u32, u32) { + self.tile_size + } + fn get_tile_index(&self, xyz: TileCoord) -> Option { if xyz.y >= self.tiles_down || xyz.x >= self.tiles_across { return None; @@ -100,6 +245,84 @@ impl Image { } } +/// Calculates the tiles that intersect with the given window. +#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] +fn tiles_intersected( + tile_size: (u32, u32), + resolution: (f64, f64), + extent: [f64; 4], + (tiles_across, tiles_down): (u32, u32), + window: [f64; 4], +) -> Vec<(u32, u32)> { + let epsilon = 1e-6; + + let tile_span_x = f64::from(tile_size.0) * resolution.0; + // resolution[1] is typically negative, use its absolute value for span calculation + let tile_span_y = f64::from(tile_size.1) * resolution.1.abs(); + + let tile_matrix_min_x = extent[0]; + // Use max Y from extent as the top edge for row calculation + let tile_matrix_max_y = extent[3]; + + let matrix_width = tiles_across; + let matrix_height = tiles_down; + + // Calculate tile index ranges based on the provided formula + let tile_min_col_f = ((window[0] - tile_matrix_min_x) / tile_span_x + epsilon).floor(); + let tile_max_col_f = ((window[2] - tile_matrix_min_x) / tile_span_x - epsilon).floor(); + let tile_min_row_f = ((tile_matrix_max_y - window[3]) / tile_span_y + epsilon).floor(); + let tile_max_row_f = ((tile_matrix_max_y - window[1]) / tile_span_y - epsilon).floor(); + + // Convert to integer type for clamping and iteration + let mut tile_min_col = tile_min_col_f as i64; + let mut tile_max_col = tile_max_col_f as i64; + let mut tile_min_row = tile_min_row_f as i64; + let mut tile_max_row = tile_max_row_f as i64; + + // Clamp minimum values to 0 + if tile_min_col < 0 { + tile_min_col = 0; + } + if tile_min_row < 0 { + tile_min_row = 0; + } + + // Clamp maximum values to matrix dimensions - 1 + let matrix_width_i64 = i64::from(matrix_width); + let matrix_height_i64 = i64::from(matrix_height); + + if tile_max_col >= matrix_width_i64 { + tile_max_col = matrix_width_i64 - 1; + } + if tile_max_row >= matrix_height_i64 { + tile_max_row = matrix_height_i64 - 1; + } + + // If the calculated range is invalid (max < min), return empty vector + if tile_max_col < tile_min_col || tile_max_row < tile_min_row { + return Vec::new(); + } + + // Convert to u32 for the final result type + let tile_min_col = tile_min_col as u32; + let tile_max_col = tile_max_col as u32; + let tile_min_row = tile_min_row as u32; + let tile_max_row = tile_max_row as u32; + + let mut covered_tiles = Vec::new(); + // Iterate through the valid tile range and collect the indexes + for row in tile_min_row..=tile_max_row { + for col in tile_min_col..=tile_max_col { + // Double check bounds (should be guaranteed by clamping, but safe) + if col < matrix_width && row < matrix_height { + covered_tiles.push((col, row)); + } + } + } + + covered_tiles +} + /// Converts RGB/RGBA tile data to PNG format. fn rgb_to_png( data: Vec, @@ -128,8 +351,7 @@ fn ensure_pixels_valid( nodata: Option, ) -> Vec { let is_padded = data_width != tile_width || data_height != tile_height; - // FIXME: why not `== 3`? - let add_alpha = components_count != 4; + let add_alpha = components_count == 3; // 1. Check if the tile is padded. If so, we need to add padding part back. // The decoded value might be smaller than the tile size. // TIFF crate always cut off the padding part, so we would need to add the padding part back. @@ -138,36 +360,15 @@ fn ensure_pixels_valid( // See https://gdal.org/en/stable/drivers/raster/gtiff.html#nodata-value if nodata.is_some() || add_alpha || is_padded { let mut result_vec = vec![0; (tile_width * tile_height * 4) as usize]; - for row in 0..data_height { - 'outer: for col in 0..data_width { - let idx_chunk = row * data_width * components_count + col * components_count; - let idx_result = row * tile_width * 4 + col * 4; - - // Copy component values one by one - for component_idx in 0..components_count { - // Before copying, check if this component == nodata. If so, skip because it's transparent. - // FIXME: Should we copy the RGB values anyway and just set alpha to 0? - // The visual result is the same (transparent), but the component values would differ. - // But it might be a little slower as we don't skip the copy. - // Source pixel: [4, 1, 2, 3] nodata: Some(1) - // Skip: - // result pixel: [4, 0, 0, 0] - // Do not skip: - // result pixel: [4, 1, 2, 0] - // So the visual result is the same, but the component values are different. - - let value = data[(idx_chunk + component_idx) as usize]; - if nodata == Some(value) { - continue 'outer; - } - // Copy this component to the result vector - result_vec[(idx_result + component_idx) as usize] = value; - } - if add_alpha { - result_vec[(idx_result + 3) as usize] = 255; // opaque - } - } - } + draw_tile( + &data, + components_count, + nodata, + (data_width, data_height), + (tile_width, tile_height), + (0, 0), + &mut result_vec, + ); result_vec } else { data @@ -200,11 +401,56 @@ fn encode_rgba_as_png( Ok(result_file_buffer) } +/// Cover a tile(rgb/rgba) on a rgba buffer +#[allow(clippy::cast_possible_truncation)] +#[allow(clippy::cast_sign_loss)] +fn draw_tile( + data: &[u8], + components_count: u32, + nodata: Option, + (data_width, data_height): (u32, u32), + (target_width, target_height): (u32, u32), + (offset_x, offset_y): (i64, i64), + target: &mut [u8], +) { + let add_alpha = components_count != 4; + for row in 0..data_height { + 'outer: for col in 0..data_width { + let idx_chunk = row * data_width * components_count + col * components_count; + let target_row = i64::from(row) + offset_y; + let target_col = i64::from(col) + offset_x; + if target_row < 0 + || target_col < 0 + || target_row >= i64::from(target_height) + || target_col >= i64::from(target_width) + { + continue 'outer; + } + let target_row = target_row as usize; + let target_col = target_col as usize; + let idx_result = target_row * target_width as usize * 4 + target_col * 4; + for component_idx in 0..components_count { + let value = data[(idx_chunk + component_idx) as usize]; + if let Some(v) = nodata + && value == v + { + continue 'outer; + } + // Copy this component to the result vector + target[idx_result + component_idx as usize] = value; + } + if add_alpha { + target[idx_result + 3_usize] = 255; // opaque + } + } + } +} + #[cfg(test)] mod tests { use std::path::Path; - use martin_tile_utils::TileCoord; + use martin_tile_utils::{TileCoord, xyz_to_bbox_webmercator}; use rstest::rstest; use crate::tiles::cog::image::Image; @@ -213,8 +459,12 @@ mod tests { fn can_calc_tile_idx() { let image = Image { ifd_index: 0, + origin: [0.0, 0.0, 0.0], + extent: [0.0, 0.0, 0.0, 0.0], tiles_across: 3, tiles_down: 3, + resolution: (1.0, 1.0), + tile_size: (256, 256), }; assert_eq!( Some(0), @@ -288,4 +538,106 @@ mod tests { insta::assert_binary_snapshot!(expected_file_path, png_bytes); } + + // test bbox which aligned with tile boundary + // these are edge cases need to be ensure + #[rstest] + #[case(0, 0, 0, 0, 0)] + #[case(1, 0, 0, 0, 0)] + #[case(1, 0, 0, 1, 1)] + #[case(2, 0, 0, 0, 0)] + #[case(2, 3, 3, 3, 3)] + #[case(2, 1, 1, 2, 2)] + #[case(2, 0, 0, 3, 3)] + #[case(3, 4, 5, 6, 7)] + #[case(4, 0, 0, 7, 0)] + #[case(4, 7, 7, 7, 12)] + #[case(1, 1, 1, 1, 1)] + fn can_get_intersected_tiles_index( + #[case] zoom: u8, + #[case] min_x: u32, + #[case] min_y: u32, + #[case] max_x: u32, + #[case] max_y: u32, + ) { + let tile_size = (256, 256); + let extent = [ + -20_037_508.342_789_2, + -20_037_508.342_789_2, + 20_037_508.342_789_2, + 20_037_508.342_789_2, + ]; + let across = 2u32.pow(u32::from(zoom)); + let down = 2u32.pow(u32::from(zoom)); + + let resolution = ( + (20_037_508.342_789_2 * 2.0) / (f64::from(across) * f64::from(tile_size.0)), + -(20_037_508.342_789_2 * 2.0) / (f64::from(down) * f64::from(tile_size.1)), + ); + + let bbox = xyz_to_bbox_webmercator(zoom, min_x, min_y, max_x, max_y); + + let actual = super::tiles_intersected(tile_size, resolution, extent, (across, down), bbox); + assert_eq!( + actual.len() as usize, + (max_x - min_x + 1) as usize * (max_y - min_y + 1) as usize + ); + for row in min_y..=max_y { + for col in min_x..=max_x { + assert!( + actual.contains(&(col, row)), + "Tile ({col}, {row}) not found in the result" + ); + } + } + } + + // test bbox which not aligned with tile boundary + #[rstest] + #[case(0, [-20_037_508.342_789_2 - 1000.0, + -20_037_508.342_789_2 - 1000.0, + 20_037_508.342_789_2 + 1000.0, + 20_037_508.342_789_2 + 1000.0], (0,0,0,0))] // bigger than extent at aoom 0, should be [0,0,0,0] + #[case(0, [-20_037_508.342_789_2 + 1000.0, + -20_037_508.342_789_2 + 1000.0, + 20_037_508.342_789_2 - 1000.0, + 20_037_508.342_789_2 - 1000.0], (0,0,0,0))] // smaller than extent at aoom 0, should be [0,0,0,0] + #[case(1, [-2000.0,1000.0,-1000.0,2000.0] ,(0,0,0,0))] + #[case(1, [1000.0,-2000.0,2000.0,-1000.0] ,(1,1,1,1))] + #[case(1, [-1000.0, + -1000.0,1000.0,1000.0], (0,0,1,1))] + fn tiles_intersected_with_bbox( + #[case] zoom: u8, + #[case] bbox: [f64; 4], + #[case] expected: (u32, u32, u32, u32), + ) { + let tile_size = (256, 256); + let extent = [ + -20_037_508.342_789_2, + -20_037_508.342_789_2, + 20_037_508.342_789_2, + 20_037_508.342_789_2, + ]; + let across = 2u32.pow(u32::from(zoom)); + let down = 2u32.pow(u32::from(zoom)); + let resolution = ( + (20_037_508.342_789_2 * 2.0) / (f64::from(across) * f64::from(tile_size.0)), + -(20_037_508.342_789_2 * 2.0) / (f64::from(down) * f64::from(tile_size.1)), + ); + + let actual = super::tiles_intersected(tile_size, resolution, extent, (across, down), bbox); + + let (min_x, min_y, max_x, max_y) = expected; + let expected_count = (max_x - min_x + 1) as usize * (max_y - min_y + 1) as usize; + assert_eq!(actual.len(), expected_count, "Unexpected tile count"); + + for row in min_y..=max_y { + for col in min_x..=max_x { + assert!( + actual.contains(&(col, row)), + "Tile ({col}, {row}) not found" + ); + } + } + } } diff --git a/martin-core/src/tiles/cog/model.rs b/martin-core/src/tiles/cog/model.rs index 33d9357d2..a535eb706 100644 --- a/martin-core/src/tiles/cog/model.rs +++ b/martin-core/src/tiles/cog/model.rs @@ -4,7 +4,7 @@ use std::path::Path; use tiff::decoder::Decoder; use tiff::tags::Tag; -use crate::tiles::cog::CogError; +use super::CogError; /// These tags define the relationship between raster space and model space. /// See [ogc doc](https://docs.ogc.org/is/19-008r4/19-008r4.html#_coordinate_transformations) for details. diff --git a/martin-core/src/tiles/cog/source.rs b/martin-core/src/tiles/cog/source.rs index 7b8add6c7..7905a0d40 100644 --- a/martin-core/src/tiles/cog/source.rs +++ b/martin-core/src/tiles/cog/source.rs @@ -5,17 +5,19 @@ use std::path::{Path, PathBuf}; use std::vec; use async_trait::async_trait; -use martin_tile_utils::{Format, TileCoord, TileData, TileInfo}; +use martin_tile_utils::{EARTH_CIRCUMFERENCE, Format, TileCoord, TileData, TileInfo}; use tiff::decoder::{ChunkType, Decoder}; use tiff::tags::Tag::{self, GdalNodata}; use tilejson::{TileJSON, tilejson}; use tracing::warn; -use crate::tiles::cog::CogError; -use crate::tiles::cog::image::Image; -use crate::tiles::cog::model::ModelInfo; +use super::CogError; +use super::image::Image; +use super::model::ModelInfo; use crate::tiles::{MartinCoreResult, Source, UrlQuery}; +const MAX_ZOOM: u8 = 30; + /// Tile source that reads from `Cloud Optimized GeoTIFF` files. #[derive(Clone, Debug)] pub struct CogSource { @@ -32,12 +34,13 @@ pub struct CogSource { nodata: Option, tilejson: TileJSON, tileinfo: TileInfo, + auto_webmercator: bool, } impl CogSource { #[expect(clippy::cast_possible_truncation)] /// Creates a new COG tile source from a file path. - pub fn new(id: String, path: PathBuf) -> Result { + pub fn new(id: String, path: PathBuf, auto_webmercator: bool) -> Result { let tileinfo = TileInfo::new(Format::Png, martin_tile_utils::Encoding::Uncompressed); let tif_file = File::open(&path).map_err(|e: std::io::Error| CogError::IoError(e, path.clone()))?; @@ -80,7 +83,14 @@ impl CogSource { .map_or_else(|_| true, |v| v & 4 != 4); // see https://www.verypdf.com/document/tiff6/pg_0036.htm if is_image { // TODO: We should not ignore mask in the next PRs - images.push(get_image(&mut decoder, &path, ifd_index)?); + images.push(get_image( + &mut decoder, + &path, + ifd_index, + origin, + extent, + (full_width, full_length), + )?); } else { warn!( "A subfile of {} is ignored in the tiff file as Martin currently does not support mask subfile in tiff. IFD={ifd_index}", @@ -96,20 +106,33 @@ impl CogSource { break; } } - let min_zoom = 0; - let max_zoom = (images.len() - 1) as u8; + + let images_len = images.len() as u8; let images: HashMap = images .into_iter() .enumerate() .map(|(idx, image)| { - let zoom = max_zoom.saturating_sub(idx as u8); + let zoom = if auto_webmercator { + nearest_web_mercator_zoom(image.resolution(), image.tile_size()) + } else { + (images_len - 1).saturating_sub(idx as u8) + }; (zoom, image) }) .collect(); + + let min_zoom = *images + .keys() + .min() + .ok_or_else(|| CogError::NoImagesFound(path.clone()))?; + let max_zoom = *images + .keys() + .max() + .ok_or_else(|| CogError::NoImagesFound(path.clone()))?; let tilejson = tilejson! { tiles: vec![], minzoom: min_zoom, - maxzoom: max_zoom + maxzoom: max_zoom, }; Ok(CogSource { id, @@ -124,10 +147,29 @@ impl CogSource { nodata, tilejson, tileinfo, + auto_webmercator, }) } } +/// find a zoom level of [WebMercatorQuad](https://docs.ogc.org/is/17-083r2/17-083r2.html#72) that is closest to the given resolution +fn nearest_web_mercator_zoom(resolution: (f64, f64), tile_size: (u32, u32)) -> u8 { + let tile_width_in_model = resolution.0 * f64::from(tile_size.0); + let mut nearest_zoom = 0u8; + let mut min_diff = f64::INFINITY; + + for zoom in 0..MAX_ZOOM { + let tile_length = EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom); + let current_diff = (tile_width_in_model - tile_length).abs(); + + if current_diff < min_diff { + min_diff = current_diff; + nearest_zoom = zoom; + } + } + nearest_zoom +} + #[async_trait] impl Source for CogSource { fn get_id(&self) -> &str { @@ -173,8 +215,14 @@ impl Source for CogSource { CogError::ZoomOutOfRange(xyz.z, self.path.clone(), self.min_zoom, self.max_zoom) })?; - let bytes = image.get_tile(&mut decoder, xyz, self.nodata, &self.path)?; - Ok(bytes) + if self.auto_webmercator { + // just clip the image to get the tile in web mercator + let bytes = image.get_tile_webmercator(&mut decoder, xyz, self.nodata, &self.path)?; + Ok(bytes) + } else { + let bytes = image.get_tile(&mut decoder, xyz, self.nodata, &self.path)?; + Ok(bytes) + } } } @@ -261,13 +309,28 @@ fn get_image( decoder: &mut Decoder, path: &Path, ifd_index: usize, + origin: [f64; 3], + extent: [f64; 4], + (width_in_model, length_in_model): (f64, f64), ) -> Result { let (tile_width, tile_height) = (decoder.chunk_dimensions().0, decoder.chunk_dimensions().1); let (image_width, image_length) = dimensions_in_pixel(decoder, path, ifd_index)?; let tiles_across = image_width.div_ceil(tile_width); let tiles_down = image_length.div_ceil(tile_height); - - Ok(Image::new(ifd_index, tiles_across, tiles_down)) + let resolution = ( + // all the images share a same extent, so to get resolution of current image, we can use the full width and lenght to divide the image width and length + width_in_model / f64::from(image_width), + length_in_model / f64::from(image_length), + ); + Ok(Image::new( + ifd_index, + extent, + origin, + tiles_across, + tiles_down, + (tile_width, tile_height), + resolution, + )) } /// Gets image pixel dimensions from TIFF decoder @@ -585,4 +648,17 @@ mod tests { assert_abs_diff_eq!(full_resolution[0], expected[0], epsilon = 0.00001); assert_abs_diff_eq!(full_resolution[1], expected[1], epsilon = 0.00001); } + + #[rstest] + #[case((9.554_628_535_644_346,-9.554_628_535_646_77),(256,256), 14)] + fn test_nearest_web_mercator_zoom( + #[case] resolution: (f64, f64), + #[case] tile_size: (u32, u32), + #[case] expected_zoom: u8, + ) { + use crate::tiles::cog::source::nearest_web_mercator_zoom; + + let result = nearest_web_mercator_zoom(resolution, tile_size); + assert_eq!(result, expected_zoom); + } } diff --git a/martin-tile-utils/src/lib.rs b/martin-tile-utils/src/lib.rs index 3303275f4..42752c537 100644 --- a/martin-tile-utils/src/lib.rs +++ b/martin-tile-utils/src/lib.rs @@ -279,6 +279,30 @@ pub fn tile_index(lng: f64, lat: f64, zoom: u8) -> (u32, u32) { (col, row) } +/// Convert min/max XYZ tile coordinates to a bounding box values in Web Mercator. +#[must_use] +pub fn xyz_to_bbox_webmercator( + zoom: u8, + min_x: u32, + min_y: u32, + max_x: u32, + max_y: u32, +) -> [f64; 4] { + assert!(zoom <= MAX_ZOOM, "zoom {zoom} must be <= {MAX_ZOOM}"); + + let tile_length = EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom); + + let left_down_bbox = tile_bbox(min_x, max_y, tile_length); + let right_top_bbox = tile_bbox(max_x, min_y, tile_length); + + [ + left_down_bbox[0], + left_down_bbox[1], + right_top_bbox[2], + right_top_bbox[3], + ] +} + /// Convert min/max XYZ tile coordinates to a bounding box values. /// /// The result is `[min_lng, min_lat, max_lng, max_lat]` diff --git a/martin/src/config/file/tiles/cog.rs b/martin/src/config/file/tiles/cog.rs index 92d90e38f..e7272d8d3 100644 --- a/martin/src/config/file/tiles/cog.rs +++ b/martin/src/config/file/tiles/cog.rs @@ -13,6 +13,17 @@ use crate::config::file::{ #[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)] pub struct CogConfig { + /// Default true + /// + /// |option |Reprojecting to WebMercatorQuad|Pro And Con| + /// |---|---| + /// |true|Server Side|1. A little bit slow as Martin would do some cliping and merging
2. No any extra configuration needed for map viewers as WebMercatorQuad is the most default| + /// |false|Client Side|1. Most efficient 2. Need extra configuration for map viewers| + /// + /// + /// As martin currently has no support for CRS not 3857, we strongly recommend to enable this option + pub auto_webmercator: Option, + #[serde(flatten, skip_serializing)] pub unrecognized: UnrecognizedValues, }