3D Tilesのアトラス化 (#620)

<!-- Close or Related Issues -->
Close #303
Close #550
Close #568
Close #569

### What I did（変更内容）
<!-- Please describe the motivation behind this PR and the changes it
introduces. -->
<!-- どのような変更をしますか？ 目的は？ -->

- 3D Tiles Sinkの拡張
  - アトラス化
    - 利用されるテクスチャのサイズに応じて作成するアトラスの大きさを変更
    - ズームレベルに応じてテクスチャの解像度を変更
- 処理速度は既存のSinkより10倍程度になっていますが、Cesiumでの描画パフォーマンスは大きく向上しています

### Notes（連絡事項）
<!-- If manual testing is required, please describe the procedure. -->
<!-- 手動での動作確認が必要なら手順を簡単に伝えてください。そのほか連絡事項など。 -->

- 以下のように変換します
```bash
cargo run -- <input_path> --sink 3dtiles -o transform=use_texture --output <output_path> --release
```
<img width="1718" alt="image"
src="https://github.com/user-attachments/assets/0d31300f-59e6-460e-97ee-1ced0bcbd261">

nusamai/src/sink/cesiumtiles/mod.rs

@@ -16,6 +16,12 @@ use std::{
 };
 
 use ahash::RandomState;
+use atlas_packer::{
+    export::{AtlasExporter as _, JpegAtlasExporter},
+    pack::TexturePacker,
+    place::{GuillotineTexturePlacer, PlacedTextureInfo, TexturePlacerConfig},
+    texture::{CroppedTexture, DownsampleFactor, TextureCache, TextureSizeCache},
+};
 use bytemuck::Zeroable;
 use earcut::{utils3d::project3d_to_2d, Earcut};
 use gltf::write_gltf_glb;
@@ -26,7 +32,9 @@ use nusamai_mvt::tileid::TileIdMethod;
 use nusamai_projection::cartesian::geodetic_to_geocentric;
 use rayon::prelude::*;
 use slice::{slice_to_tiles, SlicedFeature};
+use tempfile::tempdir;
 use tiling::{TileContent, TileTree};
+use url::Url;
 
 use crate::{
     get_parameter_value,
@@ -306,17 +314,28 @@ fn tile_writing_stage(
     let contents: Arc<Mutex<Vec<TileContent>>> = Default::default();
     let bincode_config = bincode::config::standard();
 
+    // Texture cache
+    // use default cache size
+    let texture_cache = TextureCache::new(100_000_000);
+    let texture_size_cache = TextureSizeCache::new();
+
+    // Use a temporary directory for embedding in glb.
+    let binding = tempdir().unwrap();
+    let folder_path = binding.path();
+    let texture_folder_name = "textures";
+    let atlas_dir = folder_path.join(texture_folder_name);
+    std::fs::create_dir_all(&atlas_dir)?;
+
     // Make a glTF (.glb) file for each tile
     receiver_sorted
         .into_iter()
         .par_bridge()
         .try_for_each(|(tile_id, typename, feats)| {
             feedback.ensure_not_canceled()?;
+            let (tile_zoom, tile_x, tile_y) = tile_id_conv.id_to_zxy(tile_id);
 
             // Tile information
             let (mut content, translation) = {
-                let zxy = tile_id_conv.id_to_zxy(tile_id);
-                let (tile_zoom, tile_x, tile_y) = zxy;
                 let (min_lat, max_lat) = tiling::y_slice_range(tile_zoom, tile_y);
                 let (min_lng, max_lng) = tiling::x_slice_range(
                     tile_zoom,
@@ -348,7 +367,7 @@ fn tile_writing_stage(
                     format!("{tile_zoom}/{tile_x}/{tile_y}_{normalized_typename}.glb")
                 };
                 let content = TileContent {
-                    zxy,
+                    zxy: (tile_zoom, tile_x, tile_y),
                     content_path,
                     min_lng: f64::MAX,
                     max_lng: f64::MIN,
@@ -361,16 +380,60 @@ fn tile_writing_stage(
                 (content, translation)
             };
 
-            let mut earcutter = Earcut::new();
-            let mut buf3d: Vec<[f64; 3]> = Vec::new();
-            let mut buf2d: Vec<[f64; 2]> = Vec::new(); // 2d-projected [x, y]
-            let mut index_buf: Vec<u32> = Vec::new();
-
             let mut vertices: IndexSet<[u32; 9], RandomState> = IndexSet::default(); // [x, y, z, u, v, feature_id]
             let mut primitives: gltf::Primitives = Default::default();
 
             let mut metadata_encoder = metadata::MetadataEncoder::new(schema);
 
+            // Check the size of all the textures and calculate the power of 2 of the largest size
+            let mut max_width = 0;
+            let mut max_height = 0;
+            for serialized_feat in feats.iter() {
+                feedback.ensure_not_canceled()?;
+
+                let feature = {
+                    let (feature, _): (SlicedFeature, _) =
+                        bincode::serde::decode_from_slice(serialized_feat, bincode_config)
+                            .map_err(|err| {
+                                PipelineError::Other(format!(
+                                    "Failed to deserialize a sliced feature: {:?}",
+                                    err
+                                ))
+                            })?;
+
+                    feature
+                };
+
+                for (_, orig_mat_id) in feature
+                    .polygons
+                    .iter()
+                    .zip_eq(feature.polygon_material_ids.iter())
+                {
+                    let mat = feature.materials[*orig_mat_id as usize].clone();
+                    let t = mat.base_texture.clone();
+                    if let Some(base_texture) = t {
+                        let texture_uri = base_texture.uri.to_file_path().unwrap();
+                        let texture_size = texture_size_cache.get_or_insert(&texture_uri);
+                        max_width = max_width.max(texture_size.0);
+                        max_height = max_height.max(texture_size.1);
+                    }
+                }
+            }
+            let max_width = max_width.next_power_of_two();
+            let max_height = max_height.next_power_of_two();
+
+            // initialize texture packer
+            // In rare cases, large textures also exist, so the maximum texture size is set to 8096x8096.
+            let config = TexturePlacerConfig {
+                width: max_width,
+                height: max_height,
+                padding: 0,
+            };
+            let placer = GuillotineTexturePlacer::new(config.clone());
+            let exporter = JpegAtlasExporter::default();
+            let ext = exporter.clone().get_extension().to_string();
+            let packer = Mutex::new(TexturePacker::new(placer, exporter));
+
             // For each feature
             let mut feature_id = 0;
             for serialized_feat in feats.into_iter() {
@@ -401,14 +464,14 @@ fn tile_writing_stage(
                             // - geographic to geocentric
                             // - z-up to y-up
                             // - subtract the translation
-                            // - flip the texture v-coordinate
+                            // - The origin of atlas-packer is in the lower right.
                             let (x, y, z) = geodetic_to_geocentric(&ellipsoid, lng, lat, height);
                             [
                                 x - translation[0],
                                 z - translation[1],
                                 -y - translation[2],
                                 u,
-                                1.0 - v,
+                                v,
                             ]
                         });
 
@@ -425,23 +488,100 @@ fn tile_writing_stage(
                 }
 
                 // Triangulation, etc.
-                for (poly, orig_mat_id) in feature
+                for (poly_count, (mut poly, orig_mat_id)) in feature
                     .polygons
                     .iter()
                     .zip_eq(feature.polygon_material_ids.iter())
+                    .enumerate()
                 {
-                    let num_outer_points = match poly.hole_indices().first() {
-                        Some(&v) => v as usize,
-                        None => poly.raw_coords().len(),
-                    };
+                    let mut mat = feature.materials[*orig_mat_id as usize].clone();
+                    let t = mat.base_texture.clone();
+                    if let Some(base_texture) = t {
+                        // texture packing
+                        let original_vertices = poly
+                            .raw_coords()
+                            .iter()
+                            .map(|[x, y, z, u, v]| (*x, *y, *z, *u, *v))
+                            .collect::<Vec<(f64, f64, f64, f64, f64)>>();
+                        let uv_coords = original_vertices
+                            .iter()
+                            .map(|(_, _, _, u, v)| (*u, *v))
+                            .collect::<Vec<(f64, f64)>>();
+
+                        let texture_uri = base_texture.uri.to_file_path().unwrap();
+                        let texture_size = texture_size_cache.get_or_insert(&texture_uri);
+                        let factor = apply_downsample_factor(tile_zoom);
+                        let downsample_factor = DownsampleFactor::new(&factor);
+                        let cropped_texture = CroppedTexture::new(
+                            &texture_uri,
+                            texture_size,
+                            &uv_coords,
+                            downsample_factor,
+                        );
+
+                        // Unique id required for placement in atlas
+                        let (z, x, y) = tile_id_conv.id_to_zxy(tile_id);
+                        let texture_id = format!("{}_{}_{}_{}_{}", z, x, y, feature_id, poly_count);
+                        let info: PlacedTextureInfo = packer
+                            .lock()
+                            .unwrap()
+                            .add_texture(texture_id, cropped_texture);
+
+                        let atlas_placed_uv_coords = info
+                            .placed_uv_coords
+                            .iter()
+                            .map(|(u, v)| ({ *u }, { *v }))
+                            .collect::<Vec<(f64, f64)>>();
+                        let updated_vertices = original_vertices
+                            .iter()
+                            .zip(atlas_placed_uv_coords.iter())
+                            .map(|((x, y, z, _, _), (u, v))| (*x, *y, *z, *u, *v))
+                            .collect::<Vec<(f64, f64, f64, f64, f64)>>();
+
+                        // Apply the UV coordinates placed in the atlas to the original polygon
+                        poly.transform_inplace(|&[x, y, z, _, _]| {
+                            let (u, v) = updated_vertices
+                                .iter()
+                                .find(|(x_, y_, z_, _, _)| {
+                                    (*x_ - x).abs() < 1e-6
+                                        && (*y_ - y).abs() < 1e-6
+                                        && (*z_ - z).abs() < 1e-6
+                                })
+                                .map(|(_, _, _, u, v)| (*u, *v))
+                                .unwrap();
+                            [x, y, z, u, v]
+                        });
+
+                        let atlas_file_name = info.atlas_id.to_string();
+
+                        let atlas_uri = atlas_dir
+                            .join(format!("{}/{}/{}/{}", z, x, y, atlas_file_name))
+                            .with_extension(ext.clone());
+
+                        // update material
+                        mat = material::Material {
+                            base_color: mat.base_color,
+                            base_texture: Some(material::Texture {
+                                uri: Url::from_file_path(atlas_uri).unwrap(),
+                            }),
+                        };
+                    }
 
-                    let mat = feature.materials[*orig_mat_id as usize].clone();
                     let primitive = primitives.entry(mat).or_default();
                     primitive.feature_ids.insert(feature_id as u32);
 
                     if let Some((nx, ny, nz)) =
                         calculate_normal(poly.exterior().iter().map(|v| [v[0], v[1], v[2]]))
                     {
+                        let num_outer_points = match poly.hole_indices().first() {
+                            Some(&v) => v as usize,
+                            None => poly.raw_coords().len(),
+                        };
+                        let mut earcutter = Earcut::new();
+                        let mut buf3d: Vec<[f64; 3]> = Vec::new();
+                        let mut buf2d: Vec<[f64; 2]> = Vec::new();
+                        let mut index_buf: Vec<u32> = Vec::new();
+
                         buf3d.clear();
                         buf3d.extend(poly.raw_coords().iter().map(|c| [c[0], c[1], c[2]]));
 
@@ -464,7 +604,8 @@ fn tile_writing_stage(
                                     (ny as f32).to_bits(),
                                     (nz as f32).to_bits(),
                                     (u as f32).to_bits(),
-                                    (v as f32).to_bits(),
+                                    // flip the texture v-coordinate
+                                    ((1.0 - v) as f32).to_bits(),
                                     (feature_id as f32).to_bits(), // UNSIGNED_INT can't be used for vertex attribute
                                 ];
                                 let (index, _) = vertices.insert_full(vbits);
@@ -473,10 +614,18 @@ fn tile_writing_stage(
                         }
                     }
                 }
-
                 feature_id += 1;
             }
 
+            let mut packer = packer.into_inner().unwrap();
+            packer.finalize();
+
+            // Write to atlas
+            let (z, x, y) = tile_id_conv.id_to_zxy(tile_id);
+            let atlas_path = atlas_dir.join(format!("{}/{}/{}", z, x, y));
+            fs::create_dir_all(&atlas_path)?;
+            packer.export(&atlas_path, &texture_cache, config.width, config.height);
+
             // Write to file
             let path_glb = output_path.join(Path::new(&content.content_path));
             if let Some(dir) = path_glb.parent() {
@@ -526,3 +675,12 @@ fn tile_writing_stage(
 
     Ok(())
 }
+
+fn apply_downsample_factor(z: u8) -> f32 {
+    match z {
+        0..=14 => 0.0,
+        15..=16 => 0.25,
+        17 => 0.5,
+        _ => 1.0,
+    }
+}