diff --git a/martin/src/cog/errors.rs b/martin/src/cog/errors.rs index d4b5a49c1..d40643bc1 100644 --- a/martin/src/cog/errors.rs +++ b/martin/src/cog/errors.rs @@ -58,4 +58,9 @@ pub enum CogError { "Calculating the tile origin failed for {0}: the length of ModelTiepointTag should be >= 6, or the length of ModelTransformationTag should be >= 12" )] GetOriginFailed(PathBuf), + + #[error( + "Get full resolution failed for {0}: either a valid ModelPixelScaleTag or ModelPixelScaleTag is required" + )] + GetFullResolutionFailed(PathBuf), } diff --git a/martin/src/cog/source.rs b/martin/src/cog/source.rs index 3bc2c4316..62cc64e02 100644 --- a/martin/src/cog/source.rs +++ b/martin/src/cog/source.rs @@ -25,6 +25,8 @@ struct Meta { model: ModelInfo, // The geo coords of pixel(0, 0, 0) ordering in [x, y, z] origin: [f64; 3], + // [minx, miny, maxx, maxy] in its model space coordinate system + extent: [f64; 4], zoom_and_ifd: HashMap, zoom_and_tile_across_down: HashMap, nodata: Option, @@ -325,6 +327,27 @@ fn get_meta(path: &PathBuf) -> Result { model.transformation.as_deref(), path, )?; + let (full_width_pixel, full_length_pixel) = decoder.dimensions().map_err(|e| { + CogError::TagsNotFound( + e, + vec![Tag::ImageWidth.to_u16(), Tag::ImageLength.to_u16()], + 0, // we are at ifd 0, the first image, haven't seek to others + path.clone(), + ) + })?; + let full_resolution = get_full_resolution( + model.pixel_scale.as_deref(), + model.transformation.as_deref(), + path, + )?; + let full_width = full_resolution[0] * f64::from(full_width_pixel); + let full_length = full_resolution[1] * f64::from(full_length_pixel); + let extent = get_extent( + &origin, + model.transformation.as_deref(), + (full_width_pixel, full_length_pixel), + (full_width, full_length), + ); verify_requirements(&mut decoder, &model, path)?; let mut zoom_and_ifd: HashMap = HashMap::new(); let mut zoom_and_tile_across_down: HashMap = HashMap::new(); @@ -360,6 +383,7 @@ fn get_meta(path: &PathBuf) -> Result { max_zoom: images_ifd.len() as u8 - 1, model, origin, + extent, zoom_and_ifd, zoom_and_tile_across_down, nodata, @@ -456,6 +480,87 @@ fn get_origin( } } +fn get_full_resolution( + pixel_scale: Option<&[f64]>, + transformation: Option<&[f64]>, + path: &Path, +) -> Result<[f64; 2], CogError> { + match (pixel_scale, transformation) { + // ModelPixelScaleTag = (ScaleX, ScaleY, ScaleZ) + (Some(scale), _) => Ok([scale[0], scale[1]]), + // here we adopted the 2-d matrix form based on the geotiff spec, the z-axis is dropped intentionally, see https://docs.ogc.org/is/19-008r4/19-008r4.html#_geotiff_tags_for_coordinate_transformations + // It looks like this: + /* + |- -| |- -| |- -| + | X | | a b 0 d | | I | + | | | | | | | + | Y | | e f 0 h | | J | + | | = | | | | + | Z | | 0 0 0 0 | | K | + | | | | | | | + | 1 | | 0 0 0 1 | | 1 | + |- -| |- -| |- -| + */ + (_, Some(matrix)) => { + let mut x_res = (matrix[0] * matrix[0] + matrix[4] * matrix[4]).sqrt(); + x_res = x_res.copysign(matrix[0]); + let mut y_res = (matrix[1] * matrix[1] + matrix[5] * matrix[5]).sqrt(); + // A positive y_res indicates that model space Y cordinates decrease as raster space J indices increase. This is the standard vertical relationship between raster space and model space + y_res = y_res.copysign(-matrix[5]); + Ok([x_res, y_res]) // drop the z scale directly as we don't use it + } + (None, None) => Err(CogError::GetFullResolutionFailed(path.to_path_buf())), + } +} + +fn raster2model(i: u32, j: u32, matrix: &[f64]) -> (f64, f64) { + let i = f64::from(i); + let j = f64::from(j); + let x = matrix[3] + (matrix[0] * i) + (matrix[1] * j); + let y = matrix[7] + (matrix[4] * i) + (matrix[5] * j); + (x, y) +} + +/// Computes the bounding box (`[min_x, min_y, max_x, max_y]`) based on the transformation matrix, origin, width and hieght. +/// +/// Applies a transformation matrix to corner pixels if provided; +/// otherwise, computes extent from origin and raster size in model units. +fn get_extent( + origin: &[f64; 3], + transformation: Option<&[f64]>, + (full_width_pixel, full_height_pixel): (u32, u32), + (full_width, full_height): (f64, f64), +) -> [f64; 4] { + if let Some(matrix) = transformation { + let corner_pixels = [ + (0, 0), // Top-left + (0, full_height_pixel), // Bottom-left + (full_width_pixel, 0), // Top-right + (full_width_pixel, full_height_pixel), // Bottom-right + ]; + + // Transform the first corner to initialize min/max values + let (mut min_x, mut min_y) = raster2model(corner_pixels[0].0, corner_pixels[0].1, matrix); + let mut max_x = min_x; + let mut max_y = min_y; + + // Iterate over the rest of the corners + for &(i, j) in corner_pixels.iter().skip(1) { + let (x, y) = raster2model(i, j, matrix); + min_x = min_x.min(x); + min_y = min_y.min(y); + max_x = max_x.max(x); + max_y = max_y.max(y); + } + return [min_x, min_y, max_x, max_y]; + } + let [x1, y1, _] = origin; + let x2 = x1 + full_width; + let y2 = y1 - full_height; + + [x1.min(x2), y1.min(y2), x1.max(x2), y1.max(y2)] +} + #[cfg(test)] mod tests { use std::fs::File; @@ -604,4 +709,99 @@ mod tests { } } } + + #[rstest] + #[case( + None,Some(vec![10.0, 10.0,0.0]),Some(vec![0.0, 0.0, 0.0, 1_620_750.250_8, 4_277_012.715_3, 0.0]),(512,512), + [1_620_750.250_8, 4_271_892.715_3, 1_625_870.250_8, 4_277_012.715_3] + )] + #[case( + Some(vec![ + 10.0,0.0,0.0,1_620_750.250_8, + 0.0,-10.0,0.0,4_277_012.715_3, + 0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,1.0 + ]),None,None,(512,512), + [1_620_750.250_8, 4_271_892.715_3, 1_625_870.250_8, 4_277_012.715_3] + )] + #[case( + Some(vec![ + 0.010_005_529_647_693, 0.0, 0.0, -7.583_906_932_854_38, + 0.0, -0.009_986_188_755_447_6, 0.0, 38.750_354_738_325_9, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.0 + ]), None, None, (598, 279), + [-7.583_906_9, 35.964_208_1, -1.600_600_2, 38.750_354_7] + )] + fn can_get_extent( + #[case] matrix: Option>, + #[case] pixel_scale: Option>, + #[case] tie_point: Option>, + #[case] (full_width_pixel, full_length_pixel): (u32, u32), + #[case] expected_extent: [f64; 4], + ) { + use approx::assert_abs_diff_eq; + + use crate::cog::source::{get_extent, get_full_resolution, get_origin}; + + let origin = get_origin( + tie_point.as_deref(), + matrix.as_deref(), + &PathBuf::from("not_exist.tif"), + ) + .unwrap(); + let full_resolution = get_full_resolution( + pixel_scale.as_deref(), + matrix.as_deref(), + &PathBuf::from("not_exist.tif"), + ) + .unwrap(); + + let full_width = full_resolution[0] * f64::from(full_width_pixel); + let full_length = full_resolution[1] * f64::from(full_length_pixel); + + let extent = get_extent( + &origin, + matrix.as_deref(), + (full_width_pixel, full_length_pixel), + (full_width, full_length), + ); + + assert_abs_diff_eq!(extent[0], expected_extent[0], epsilon = 0.00001); + assert_abs_diff_eq!(extent[1], expected_extent[1], epsilon = 0.00001); + assert_abs_diff_eq!(extent[2], expected_extent[2], epsilon = 0.00001); + assert_abs_diff_eq!(extent[3], expected_extent[3], epsilon = 0.00001); + } + + #[rstest] + #[case( + None,Some(vec![118.450_587_6, 118.450_587_6, 0.0]), [118.450_587_6, 118.450_587_6] + )] + #[case( + None,Some(vec![100.00, -100.0]), [100.0, -100.0] + )] + #[ + case( + Some(vec![ + 0.010_005_529_647_693_3, 0.0, 0.0, -7.583_906_932_854_38, 0.0, -0.009_986_188_755_447_63, 0.0, 38.750_354_738_325_9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]), + None, [0.010_005_529_647_693, 0.009_986_188_755_448]) + ] + fn can_get_full_resolution( + #[case] matrix: Option>, + #[case] pixel_scale: Option>, + #[case] expected: [f64; 2], + ) { + use approx::assert_abs_diff_eq; + + use crate::cog::source::get_full_resolution; + + let full_resolution = get_full_resolution( + pixel_scale.as_deref(), + matrix.as_deref(), + &PathBuf::from("not_exist.tif"), + ) + .unwrap(); + assert_abs_diff_eq!(full_resolution[0], expected[0], epsilon = 0.00001); + assert_abs_diff_eq!(full_resolution[1], expected[1], epsilon = 0.00001); + } }