OBJ sinkの作成 ジオメトリ出力のみ

app/src-tauri/src/main.rs

@@ -15,8 +15,9 @@
     sink::{
         cesiumtiles::CesiumTilesSinkProvider, czml::CzmlSinkProvider, geojson::GeoJsonSinkProvider,
         gltf::GltfSinkProvider, gpkg::GpkgSinkProvider, kml::KmlSinkProvider,
-        minecraft::MinecraftSinkProvider, mvt::MvtSinkProvider, ply::StanfordPlySinkProvider,
-        serde::SerdeSinkProvider, shapefile::ShapefileSinkProvider, DataSinkProvider,
+        minecraft::MinecraftSinkProvider, mvt::MvtSinkProvider, obj::ObjSinkProvider,
+        ply::StanfordPlySinkProvider, serde::SerdeSinkProvider, shapefile::ShapefileSinkProvider,
+        DataSinkProvider,
     },
     source::{citygml::CityGmlSourceProvider, DataSourceProvider},
     transformer::{
@@ -123,6 +124,7 @@
         "ply" => Some(Box::new(StanfordPlySinkProvider {})),
         "cesiumtiles" => Some(Box::new(CesiumTilesSinkProvider {})),
         "minecraft" => Some(Box::new(MinecraftSinkProvider {})),
+        "obj" => Some(Box::new(ObjSinkProvider {})),
         _ => None,
     }
 }

nusamai/src/lib.rs

@@ -17,4 +17,5 @@ pub static BUILTIN_SINKS: &[&dyn sink::DataSinkProvider] = &[
     &sink::shapefile::ShapefileSinkProvider {},
     &sink::noop::NoopSinkProvider {},
     &sink::minecraft::MinecraftSinkProvider {},
+    &sink::obj::ObjSinkProvider {},
 ];

nusamai/src/sink/mod.rs

@@ -7,6 +7,7 @@ pub mod gltf;
 pub mod gpkg;
 pub mod kml;
 pub mod minecraft;
+pub mod obj;
 pub mod mvt;
 pub mod noop;
 pub mod ply;

nusamai/src/sink/obj/mod.rs

@@ -0,0 +1,355 @@
+//! obj sink
+
+use std::{f64::consts::FRAC_PI_2, io::Write, path::PathBuf, sync::Mutex};
+
+use ahash::{HashMap, RandomState};
+use earcut::{utils3d::project3d_to_2d, Earcut};
+use flatgeom::MultiPolygon;
+use indexmap::IndexSet;
+use nusamai_citygml::{
+    object::{ObjectStereotype, Value},
+    schema::Schema,
+    GeometryType,
+};
+
+use serde::{Deserialize, Serialize};
+
+use glam::{DMat4, DVec3, DVec4};
+
+use crate::{
+    get_parameter_value,
+    parameters::*,
+    pipeline::{Feedback, PipelineError, Receiver, Result},
+    sink::{DataRequirements, DataSink, DataSinkProvider, SinkInfo},
+    transformer::{TransformerOption, TransformerRegistry},
+};
+
+use itertools::Itertools;
+use nusamai_projection::cartesian::geodetic_to_geocentric;
+use rayon::iter::{IntoParallelIterator, ParallelBridge, ParallelIterator};
+
+pub struct ObjSinkProvider {}
+
+impl DataSinkProvider for ObjSinkProvider {
+    fn info(&self) -> SinkInfo {
+        SinkInfo {
+            id_name: "obj".to_string(),
+            name: "OBJ".to_string(),
+        }
+    }
+
+    fn parameters(&self) -> Parameters {
+        let mut params = Parameters::new();
+        params.define(
+            "@output".into(),
+            ParameterEntry {
+                description: "Output file path".into(),
+                required: true,
+                parameter: ParameterType::FileSystemPath(FileSystemPathParameter {
+                    value: None,
+                    must_exist: false,
+                }),
+                label: None,
+            },
+        );
+        params
+    }
+
+    fn available_transformer(&self) -> TransformerRegistry {
+        let settings: TransformerRegistry = TransformerRegistry::new();
+
+        settings
+    }
+    fn create(&self, params: &Parameters) -> Box<dyn DataSink> {
+        let output_path = get_parameter_value!(params, "@output", FileSystemPath);
+        let transform_options = self.available_transformer();
+
+        Box::<ObjSink>::new(ObjSink {
+            output_path: output_path.as_ref().unwrap().into(),
+            transform_settings: transform_options,
+        })
+    }
+}
+
+pub struct ObjSink {
+    output_path: PathBuf,
+    transform_settings: TransformerRegistry,
+}
+
+#[derive(Debug)]
+pub struct BoundingVolume {
+    pub min_lng: f64,
+    pub max_lng: f64,
+    pub min_lat: f64,
+    pub max_lat: f64,
+    pub min_height: f64,
+    pub max_height: f64,
+}
+
+impl BoundingVolume {
+    fn update(&mut self, other: &Self) {
+        self.min_lng = self.min_lng.min(other.min_lng);
+        self.max_lng = self.max_lng.max(other.max_lng);
+        self.min_lat = self.min_lat.min(other.min_lat);
+        self.max_lat = self.max_lat.max(other.max_lat);
+        self.min_height = self.min_height.min(other.min_height);
+        self.max_height = self.max_height.max(other.max_height);
+    }
+}
+
+impl Default for BoundingVolume {
+    fn default() -> Self {
+        Self {
+            min_lng: f64::MAX,
+            max_lng: f64::MIN,
+            min_lat: f64::MAX,
+            max_lat: f64::MIN,
+            min_height: f64::MAX,
+            max_height: f64::MIN,
+        }
+    }
+}
+
+#[derive(Serialize, Deserialize, Debug, Clone)]
+pub struct Feature {
+    // polygons [x, y, z, u, v]
+    pub polygons: MultiPolygon<'static, [f64; 5]>,
+    // feature_id
+    pub feature_id: Option<u32>,
+}
+
+type ClassifiedFeatures = HashMap<String, ClassFeatures>;
+
+#[derive(Default)]
+struct ClassFeatures {
+    features: Vec<Feature>,
+    bounding_volume: BoundingVolume,
+}
+
+impl DataSink for ObjSink {
+    fn make_requirements(&mut self, properties: Vec<TransformerOption>) -> DataRequirements {
+        let default_requirements: DataRequirements = DataRequirements {
+            resolve_appearance: true,
+            key_value: crate::transformer::KeyValueSpec::JsonifyObjectsAndArrays,
+            ..Default::default()
+        };
+
+        for prop in properties {
+            let _ = &self
+                .transform_settings
+                .update_transformer(&prop.key, prop.is_enabled);
+        }
+
+        self.transform_settings.build(default_requirements)
+    }
+
+    fn run(&mut self, upstream: Receiver, feedback: &Feedback, schema: &Schema) -> Result<()> {
+        let ellipsoid = nusamai_projection::ellipsoid::wgs84();
+
+        let classified_features: Mutex<ClassifiedFeatures> = Default::default();
+
+        // Construct a Feature classified by typename from Entity
+        // The coordinates of polygon store the actual coordinate values (WGS84) and UV coordinates, not the index.
+        let _ = upstream.into_iter().par_bridge().try_for_each(|parcel| {
+            feedback.ensure_not_canceled()?;
+
+            let entity = parcel.entity;
+
+            // entity must be a Feature
+            let Value::Object(obj) = &entity.root else {
+                return Ok(());
+            };
+            let ObjectStereotype::Feature { geometries, .. } = &obj.stereotype else {
+                return Ok(());
+            };
+
+            let geom_store = entity.geometry_store.read().unwrap();
+            if geom_store.multipolygon.is_empty() {
+                return Ok(());
+            }
+
+            let mut feature = Feature {
+                polygons: MultiPolygon::new(),
+                feature_id: None,
+            };
+
+            let mut local_bvol = BoundingVolume::default();
+
+            geometries.iter().for_each(|entry| {
+                match entry.ty {
+                    GeometryType::Solid | GeometryType::Surface | GeometryType::Triangle => {
+                        // extract the polygon and UV
+                        for (idx_poly, poly_uv) in
+                            geom_store
+                                .multipolygon
+                                .iter_range(entry.pos as usize..(entry.pos + entry.len) as usize)
+                                .zip_eq(geom_store.polygon_uvs.iter_range(
+                                    entry.pos as usize..(entry.pos + entry.len) as usize,
+                                ))
+                        {
+                            // convert to idx_poly to polygon
+                            let poly = idx_poly.transform(|c| geom_store.vertices[*c as usize]);
+
+                            let mut ring_buffer: Vec<[f64; 5]> = Vec::new();
+
+                            poly.rings()
+                                .zip_eq(poly_uv.rings())
+                                .for_each(|(ring, uv_ring)| {
+                                    ring.iter_closed().zip_eq(uv_ring.iter_closed()).for_each(
+                                        |(c, uv)| {
+                                            let [lng, lat, height] = c;
+                                            ring_buffer.push([lng, lat, height, uv[0], uv[1]]);
+
+                                            local_bvol.min_lng = local_bvol.min_lng.min(lng);
+                                            local_bvol.max_lng = local_bvol.max_lng.max(lng);
+                                            local_bvol.min_lat = local_bvol.min_lat.min(lat);
+                                            local_bvol.max_lat = local_bvol.max_lat.max(lat);
+                                            local_bvol.min_height =
+                                                local_bvol.min_height.min(height);
+                                            local_bvol.max_height =
+                                                local_bvol.max_height.max(height);
+                                        },
+                                    );
+
+                                    feature.polygons.add_exterior(ring_buffer.drain(..));
+                                });
+                        }
+                    }
+                    GeometryType::Curve => {
+                        // TODO: implement
+                    }
+                    GeometryType::Point => {
+                        // TODO: implement
+                    }
+                }
+            });
+
+            {
+                let mut locked_features = classified_features.lock().unwrap();
+                let feats = locked_features.entry(obj.typename.to_string()).or_default();
+                feats.features.push(feature);
+                feats.bounding_volume.update(&local_bvol);
+            }
+
+            Ok::<(), PipelineError>(())
+        });
+
+        let classified_features = classified_features.into_inner().unwrap();
+
+        // Bounding volume for the entire dataset
+        let global_bvol = {
+            let mut global_bvol = BoundingVolume::default();
+            for features in classified_features.values() {
+                global_bvol.update(&features.bounding_volume);
+            }
+            global_bvol
+        };
+
+        let transform_matrix = {
+            let bounds = &global_bvol;
+            let center_lng = (bounds.min_lng + bounds.max_lng) / 2.0;
+            let center_lat = (bounds.min_lat + bounds.max_lat) / 2.0;
+
+            let psi = ((1. - ellipsoid.e_sq()) * center_lat.to_radians().tan()).atan();
+
+            let (tx, ty, tz) = geodetic_to_geocentric(&ellipsoid, center_lng, center_lat, 0.);
+            let h = (tx * tx + ty * ty + tz * tz).sqrt();
+
+            DMat4::from_translation(DVec3::new(0., -h, 0.))
+                * DMat4::from_rotation_x(-(FRAC_PI_2 - psi))
+                * DMat4::from_rotation_y((-center_lng - 90.).to_radians())
+        };
+        let _ = transform_matrix.inverse();
+
+        classified_features
+            .into_par_iter()
+            .try_for_each(|(typename, mut features)| {
+                feedback.ensure_not_canceled()?;
+
+                // Triangulation
+                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();
+                let mut triangles = Vec::new();
+
+                for feature in features.features.iter_mut() {
+                    feedback.ensure_not_canceled()?;
+
+                    feature
+                        .polygons
+                        .transform_inplace(|&[lng, lat, height, u, v]| {
+                            let (x, y, z) = geodetic_to_geocentric(&ellipsoid, lng, lat, height);
+                            let v_xyz = DVec4::new(x, z, -y, 1.0);
+                            let v_enu = transform_matrix * v_xyz;
+                            [v_enu[0], v_enu[1], v_enu[2], u, 1.0 - v]
+                        });
+
+                    for poly in feature.polygons.iter() {
+                        let num_outer = match poly.hole_indices().first() {
+                            Some(&v) => v as usize,
+                            None => poly.raw_coords().len(),
+                        };
+
+                        buf3d.clear();
+                        buf3d.extend(poly.raw_coords().iter().map(|c| [c[0], c[1], c[2]]));
+
+                        if project3d_to_2d(&buf3d, num_outer, &mut buf2d) {
+                            // earcut
+                            earcutter.earcut(
+                                buf2d.iter().cloned(),
+                                poly.hole_indices(),
+                                &mut index_buf,
+                            );
+                            triangles.extend(index_buf.iter().map(|&idx| buf3d[idx as usize]));
+                        }
+                    }
+                }
+
+                // make vertices and indices
+                let mut vertices: IndexSet<[u64; 3], RandomState> = IndexSet::default();
+                let indices: Vec<_> = triangles
+                    .iter()
+                    .map(|[x, y, z]| {
+                        let vbits = [(x).to_bits(), (y).to_bits(), (z).to_bits()];
+                        let (index, _) = vertices.insert_full(vbits);
+                        index as u32
+                    })
+                    .collect();
+
+                println!("{:?} {:?}", vertices.len(), indices.len());
+
+                feedback.ensure_not_canceled()?;
+
+                // Write to file
+                let file_path = {
+                    let filename = format!("{}.obj", typename.replace(':', "_"));
+                    self.output_path.join(filename)
+                };
+                let file = std::fs::File::create(file_path)?;
+                let mut writer = std::io::BufWriter::new(file);
+
+                // Writing vertex data
+                for vertex in &vertices {
+                    let [vx, vy, vz] = vertex;
+                    let vx = f64::from_bits(*vx);
+                    let vy = f64::from_bits(*vy);
+                    let vz = f64::from_bits(*vz);
+                    writeln!(writer, "v {} {} {}", vx, vy, vz)?;
+                }
+
+                // Writing of surface data (index starts at 1, so +1 is used)
+                for face in indices.chunks(3) {
+                    if let [i1, i2, i3] = face {
+                        writeln!(writer, "f {} {} {}", i1 + 1, i2 + 1, i3 + 1)?;
+                    }
+                }
+
+                writer.flush()?;
+
+                Ok::<(), PipelineError>(())
+            })?;
+
+        Ok(())
+    }
+}