Add get origin func to COG

martin/Cargo.toml

@@ -119,6 +119,7 @@ static-files = { workspace = true, optional = true }
 
 [dev-dependencies]
 anyhow.workspace = true
+approx.workspace = true
 criterion.workspace = true
 ctor.workspace = true
 indoc.workspace = true

martin/src/cog/errors.rs

@@ -3,6 +3,7 @@ use std::path::PathBuf;
 use png::EncodingError;
 use tiff::TiffError;
 
+#[non_exhaustive]
 #[derive(thiserror::Error, Debug)]
 pub enum CogError {
     #[error("Couldn't decode {1} as tiff file: {0}")]
@@ -52,4 +53,9 @@ pub enum CogError {
         "The pixel size of the image {0} is not squared, the x_scale is {1}, the y_scale is {2}"
     )]
     NonSquaredImage(PathBuf, f64, f64),
+
+    #[error(
+        "Calculating the tile origin failed for {0}: the length of ModelTiepointTag should be >= 6, or the length of ModelTransformationTag should be >= 12"
+    )]
+    GetOriginFailed(PathBuf),
 }

martin/src/cog/source.rs

@@ -23,6 +23,8 @@ struct Meta {
     min_zoom: u8,
     max_zoom: u8,
     model: ModelInfo,
+    // The geo coords of pixel(0, 0, 0) ordering in [x, y, z]
+    origin: [f64; 3],
     zoom_and_ifd: HashMap<u8, usize>,
     zoom_and_tile_across_down: HashMap<u8, (u32, u32)>,
     nodata: Option<f64>,
@@ -318,6 +320,11 @@ fn get_meta(path: &PathBuf) -> Result<Meta, FileError> {
         .map_err(|e| CogError::InvalidTiffFile(e, path.clone()))?
         .with_limits(tiff::decoder::Limits::unlimited());
     let model = ModelInfo::decode(&mut decoder, path);
+    let origin = get_origin(
+        model.tie_points.as_deref(),
+        model.transformation.as_deref(),
+        path,
+    )?;
     verify_requirements(&mut decoder, &model, path)?;
     let mut zoom_and_ifd: HashMap<u8, usize> = HashMap::new();
     let mut zoom_and_tile_across_down: HashMap<u8, (u32, u32)> = HashMap::new();
@@ -352,6 +359,7 @@ fn get_meta(path: &PathBuf) -> Result<Meta, FileError> {
         min_zoom: 0,
         max_zoom: images_ifd.len() as u8 - 1,
         model,
+        origin,
         zoom_and_ifd,
         zoom_and_tile_across_down,
         nodata,
@@ -416,6 +424,38 @@ fn get_images_ifd(decoder: &mut Decoder<File>, path: &Path) -> Vec<usize> {
     res
 }
 
+fn get_origin(
+    tie_points: Option<&[f64]>,
+    transformation: Option<&[f64]>,
+    path: &Path,
+) -> Result<[f64; 3], CogError> {
+    // From geotiff spec: "This matrix tag should not be used if the ModelTiepointTag and the ModelPixelScaleTag are already defined"
+    // See more in https://docs.ogc.org/is/19-008r4/19-008r4.html#_geotiff_tags_for_coordinate_transformations
+    match (tie_points, transformation) {
+        // From geotiff spec: "If possible, the first tiepoint placed in this tag shall be the one establishing the location of the point (0,0) in raster space"
+        (Some(points), _) if points.len() >= 6 => Ok([points[3], points[4], points[5]]),
+
+        // coords =     matrix  * coords
+        // |- -|     |-       -|  |- -|
+        // | X |     | a b c d |  | I |
+        // | | |     |         |  |   |
+        // | Y |     | e f g h |  | J |
+        // |   |  =  |         |  |   |
+        // | Z |     | i j k l |  | K |
+        // | | |     |         |  |   |
+        // | 1 |     | m n o p |  | 1 |
+        // |- -|     |-       -|  |- -|
+
+        // The (I,J,K) of origin is (0,0,0), so:
+        //
+        //    x = I*a + J*b + K*c + 1*d => d => matrix[3]
+        //    y = I*e + J*f + k*g + 1*h => h => matrix[7]
+        //    z = I*i + J*j + K*k + 1*l => l => matrix[11]
+        (_, Some(matrix)) if matrix.len() >= 12 => Ok([matrix[3], matrix[7], matrix[11]]),
+        _ => Err(CogError::GetOriginFailed(path.to_path_buf())),
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use std::{fs::File, path::PathBuf};
@@ -521,4 +561,45 @@ mod tests {
         "###);
         assert_yaml_snapshot!(transformation, @"~");
     }
+
+    #[rstest]
+    #[case(
+        Some(vec![0.0, 0.0, 0.0, 1_620_750.250_8, 4_277_012.715_3, 0.0]),None,
+        Some([1_620_750.250_8, 4_277_012.715_3, 0.0])
+    )]
+    #[case(
+        None,Some(vec![
+            0.0, 100.0, 0.0, 400_000.0, 100.0, 0.0, 0.0, 500_000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+            0.0, 1.0,
+        ]),
+        Some([400_000.0, 500_000.0, 0.0])
+    )]
+    #[case(None, None, None)]
+    fn can_get_origin(
+        #[case] tie_point: Option<Vec<f64>>,
+        #[case] matrix: Option<Vec<f64>>,
+        #[case] expected: Option<[f64; 3]>,
+    ) {
+        use approx::assert_abs_diff_eq;
+
+        let origin = super::get_origin(
+            tie_point.as_deref(),
+            matrix.as_deref(),
+            &PathBuf::from("not_exist.tif"),
+        )
+        .ok();
+        match (origin, expected) {
+            (Some(o), Some(e)) => {
+                assert_abs_diff_eq!(o[0], e[0]);
+                assert_abs_diff_eq!(o[1], e[1]);
+                assert_abs_diff_eq!(o[2], e[2]);
+            }
+            (None, None) => {
+                // Both are None, which is expected
+            }
+            _ => {
+                panic!("Origin {:?} does not match expected {:?}", origin, expected);
+            }
+        }
+    }
 }