unproject test用

sevendev · sevendev · commit 973ba4a8437c · 2025-04-15T17:08:15.000+09:00
diff --git a/include/plateau/geometry/geo_coordinate.h b/include/plateau/geometry/geo_coordinate.h
@@ -37,7 +37,6 @@ namespace plateau::geometry {
         GeoCoordinate operator/(GeoCoordinate op) const;
     };
 
-
     /**
      * @enum CoordinateSystem
      *
@@ -99,4 +98,48 @@ namespace plateau::geometry {
             };
         }
     };
+
+    struct ReferencePointFactory {
+        // EPSGごとの基準点取得
+        static void GetReferencePoint(double epsg, GeoCoordinate& point) {
+            // EPSG:10169は、日本測地系2011（JGD2011）に基づく平面直角座標系第VIII系を指します。この座標系の基準点は、緯度36度、経度138.5度に設定されています。
+            if (epsg == 10169) {
+                point = GeoCoordinate(36, 138.5, 0);
+                return;
+            }
+            else if (epsg == 10162) {
+                point = GeoCoordinate(33, 129.5, 0);
+                return;
+            }
+            else if (epsg == 10163) {
+                point = GeoCoordinate(33, 131, 0);
+                return;
+            }
+            point = GeoCoordinate();
+        }
+
+        static GeoCoordinate GetReferencePoint(double epsg) {
+            // EPSG:10169は、日本測地系2011（JGD2011）に基づく平面直角座標系第VIII系を指します。この座標系の基準点は、緯度36度、経度138.5度に設定されています。
+            if (epsg == 10169) {
+                return GeoCoordinate(36, 138.5, 0);
+            }
+            else if (epsg == 10162) {
+                return GeoCoordinate(33, 129.5, 0);
+            }
+            else if (epsg == 10163) {
+                return GeoCoordinate(33, 131, 0);
+            }
+            return GeoCoordinate();
+        }
+
+        // 極座標系・平面直角座標系判定
+        static bool IsPolarCoordinateSystem(double epsg) {
+            // 平面直角座標系の区分についてはこちらを参照してください :
+            // https://www.mlit.go.jp/plateaudocument/toc9/toc9_08/toc9_08_04/
+            if (epsg >= 10162 && epsg <= 10174) {
+                return false;
+            }
+            return true;
+        }
+    };
 }
diff --git a/include/plateau/geometry/geo_reference.h b/include/plateau/geometry/geo_reference.h
@@ -27,6 +27,7 @@ namespace plateau::geometry {
         TVec3d projectWithoutAxisConvert(const TVec3d& lat_lon) const;
         TVec3d convertAxisToENU(const TVec3d& vertex) const;
         TVec3d convert(const TVec3d& lat_lon, const bool convert_axis = true, const bool project = true) const;
+        TVec3d convert(const TVec3d& lat_lon, const bool convert_axis = true, const double epsg = 6697) const;
         static TVec3d convertAxisFromENUTo(CoordinateSystem axis, const TVec3d& vertex);
         static TVec3d convertAxisToENU(CoordinateSystem axis, const TVec3d& vertex); 
 
@@ -38,6 +39,7 @@ namespace plateau::geometry {
         void setZoneID(int value);
         float getUnitScale() const;
         CoordinateSystem getCoordinateSystem() const;
+        TVec3d getOffset(const double epsg) const;
 
     private:
         TVec3d reference_point_;
diff --git a/include/plateau/polygon_mesh/mesh_extract_options.h b/include/plateau/polygon_mesh/mesh_extract_options.h
@@ -38,7 +38,8 @@ namespace plateau::polygonMesh {
                 attach_map_tile(true),
                 map_tile_zoom_level(15),
                 map_tile_url("https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg"),
-                is_polar_coordinate_system(true)
+                is_polar_coordinate_system(true),
+                epsg_code(6697)
                 {}
 
     public:
@@ -117,5 +118,10 @@ namespace plateau::polygonMesh {
 		 * 極座標系か、平面直角座標系かを指定します。
 		 */
         bool is_polar_coordinate_system;
+
+        /**
+         * 平面直角座標系は、EPSGコードに応じて基準点を取得します。
+         */
+        double epsg_code;
     };
 }
diff --git a/include/plateau/polygon_mesh/polygon_mesh_utils.h b/include/plateau/polygon_mesh/polygon_mesh_utils.h
@@ -37,6 +37,8 @@ namespace plateau::polygonMesh {
          */
         static TVec3d cityObjPos(const citygml::CityObject& city_obj);
 
+        //static TVec3d cityObjPos(const citygml::CityObject& city_obj, const double epsg, int coordinate_zone_id);
+
         /**
          * cityObjのポリゴンであり、頂点数が1以上であるものを検索します。
          * 最初に見つかったポリゴンを返します。なければ nullptr を返します。
diff --git a/src/geometry/geo_reference.cpp b/src/geometry/geo_reference.cpp
@@ -41,6 +41,41 @@ namespace plateau::geometry {
         return converted_point;
     }
 
+    TVec3d GeoReference::convert(const TVec3d& lat_lon, const bool convert_axis, const double epsg) const {
+        //平面直角座標変換、座標軸変換をフラグに応じてスキップします。
+        TVec3d point = lat_lon;
+
+        // 極座標系・平面直角座標系判定
+        const bool project = ReferencePointFactory::IsPolarCoordinateSystem(epsg);
+        // 平面直角座標系に変換
+        if (project) {
+            PolarToPlaneCartesian().project(point, zone_id_);
+            if (!convert_axis) {
+                // 座標軸変換をしない場合
+                TVec3 converted_point = point / unit_scale_ - convertAxisToENU(coordinate_system_, reference_point_);
+                return converted_point;
+            }
+            // 座標軸変換をする場合
+            TVec3 converted_point = convertAxisFromENUTo(coordinate_system_, point);
+            converted_point = converted_point / unit_scale_ - reference_point_;
+            return converted_point;
+        }
+
+        // 平面直角座標系の場合基準座標値をオフセット
+        //GeoCoordinate ref_point;
+        //ReferencePointFactory::GetReferencePoint(epsg, ref_point);
+        //plateau::geometry::GeoReference geo_ref(zone_id_);
+        //const auto& offset = geo_ref.project(ref_point);
+        const auto& offset = getOffset(epsg);
+        return point + offset;
+    }
+
+    TVec3d GeoReference::getOffset(const double epsg) const {        
+        GeoCoordinate ref_point = ReferencePointFactory::GetReferencePoint(epsg);
+        plateau::geometry::GeoReference geo_ref(zone_id_);
+        return geo_ref.project(ref_point);
+    }
+
     TVec3d GeoReference::convertAxisToENU(const TVec3d& vertex) const {
         return convertAxisToENU(getCoordinateSystem(), vertex);
     }
diff --git a/src/polygon_mesh/area_mesh_factory.cpp b/src/polygon_mesh/area_mesh_factory.cpp
@@ -24,16 +24,25 @@ namespace {
             if (!options.is_polar_coordinate_system) {
                 // 平面直角座標系の判定
                 plateau::geometry::GeoReference geo_ref(options.coordinate_zone_id, options.reference_point, options.unit_scale, options.mesh_axes);
+                //plateau::geometry::GeoReference geo_ref(options.coordinate_zone_id);
                 try {
                     auto pos = PolygonMeshUtils::cityObjPos(city_obj);
-                    if (extent.contains(geo_ref.unproject(pos)))
+                    const auto unprojected = geo_ref.unproject(pos);
+                    if (extent.contains(unprojected))
                         return false;
                 }
                 catch (std::invalid_argument& e) {}
-            }
 
-            if (extent.contains(city_obj))
-                return false;
+                //Temp Debug =======================
+
+                //return false;
+
+                //Temp Debug =======================
+            }
+            else {
+                if (extent.contains(city_obj))
+                    return false;
+            }
         }
 
         return true;
diff --git a/src/polygon_mesh/mesh_extractor.cpp b/src/polygon_mesh/mesh_extractor.cpp
@@ -28,15 +28,70 @@ namespace {
 
         if (!options.is_polar_coordinate_system) {
             // 平面直角座標系の判定
-            plateau::geometry::GeoReference geo_ref(options.coordinate_zone_id, options.reference_point, options.unit_scale, options.mesh_axes);  
+
+            //plateau::geometry::GeoCoordinate refPoint = plateau::geometry::ReferencePointFactory::GetReferencePoint(options.epsg_code);
+            //plateau::geometry::GeoReference geo_ref1(options.coordinate_zone_id);
+            //const auto epsg_offset = geo_ref1.getOffset(options.epsg_code);
+            //geo_ref1.setReferencePoint(epsg_offset);
+
+            plateau::geometry::GeoReference geo_ref2(options.coordinate_zone_id, options.reference_point, options.unit_scale, options.mesh_axes); 
+            const auto epsg_offset = geo_ref2.getOffset(options.epsg_code);
+
+            //const auto diff = options.reference_point - epsg_offset;
+            //plateau::geometry::GeoReference geo_ref3(options.coordinate_zone_id, options.reference_point - diff , options.unit_scale, options.mesh_axes);
+
+            const auto meshcode_offset = GeoReference::convertAxisToENU(options.mesh_axes, options.reference_point);
+            //const auto original_ref = geo_ref2.unproject(options.reference_point);
+
+            plateau::geometry::GeoReference geo_ref3(options.coordinate_zone_id, TVec3d(), options.unit_scale, options.mesh_axes);
+            const auto original_ref = geo_ref3.unproject(options.reference_point);
+
+            plateau::geometry::GeoReference geo(options.coordinate_zone_id);
+            GeoCoordinate ref_point = ReferencePointFactory::GetReferencePoint(options.epsg_code); //EPSGの基準点
+            const auto prj = geo.project(ref_point);
+            plateau::geometry::GeoReference geo_ref4(options.coordinate_zone_id, prj, options.unit_scale); 
+
+            plateau::geometry::GeoReference geo2(options.coordinate_zone_id);
+            GeoCoordinate ref_point2(37.4258, 138.7378, 0); //08EE751の中心
+            const auto prj2 = geo2.project(ref_point2);
+            plateau::geometry::GeoReference geo_ref5(options.coordinate_zone_id, prj2, options.unit_scale);
+
+            plateau::geometry::GeoReference geo3(options.coordinate_zone_id);
+            GeoCoordinate ref_point3(36, 138.5, 0); //EPSG:10169 の基準点
+            const auto prj3 = geo3.project(ref_point3);
+            plateau::geometry::GeoReference geo_ref6(options.coordinate_zone_id, prj3, options.unit_scale);
+
             try {
-                auto pos = geo_ref.unproject(PolygonMeshUtils::cityObjPos(city_obj));
+                const auto pos = PolygonMeshUtils::cityObjPos(city_obj);   
+
+                const auto unprojected = geo_ref4.unproject(pos);
+
+                const auto unprojected1 = geo_ref2.unproject(pos + epsg_offset);
+                //const auto unprojected2 = geo_ref2.unproject(epsg_offset - pos);
+                //const auto unprojected2 = geo_ref2.unproject(pos);
+                //const auto unprojected3 = geo_ref4.unproject(pos + epsg_offset);
+
+                const auto unprojected4 = geo_ref5.unproject(pos + epsg_offset);
+                //const auto unprojected5 = geo_ref5.unproject(epsg_offset- pos);
+
+                const auto unprojected5 = geo_ref6.unproject(pos + epsg_offset);
+
                 for (const auto& extent : extents) {
-                    if (extent.contains(pos))
+                    if (extent.contains(unprojected1))
                         return false;
                 }
+
+                //auto pos = PolygonMeshUtils::cityObjPos(city_obj);
+                //for (const auto& extent : extents) {
+                //    if (extent.contains(geo_ref.unproject(pos)))
+                //        return false;
+                //}
             }
             catch (std::invalid_argument& e) {}
+
+            //Temp Debug =======================
+            //return false;
+            //Temp Debug =======================
         }
         else {
             for (const auto& extent : extents) {
diff --git a/src/polygon_mesh/mesh_factory.cpp b/src/polygon_mesh/mesh_factory.cpp
@@ -32,7 +32,7 @@ namespace plateau::polygonMesh {
          */
         void cityGmlPolygonToMesh(
             const Polygon& polygon, const std::string& gml_path,
-            const GeoReference& geo_reference, Mesh& out_mesh, bool is_polar_coordinate_system) {
+            const GeoReference& geo_reference, Mesh& out_mesh, double epsg) {
 
             // マージ対象の情報を取得します。ここでの頂点は極座標です。
             const auto& vertices_lat_lon = polygon.getVertices();
@@ -54,7 +54,19 @@ namespace plateau::polygonMesh {
             auto& out_vertices = out_mesh.getVertices();
             out_vertices.reserve(vertices_lat_lon.size());
             for (const auto& lat_lon : vertices_lat_lon) {
-                auto xyz = geo_reference.convert(lat_lon, false, is_polar_coordinate_system);
+
+                //auto xyz = geo_reference.convert(lat_lon, false, is_polar_coordinate_system);
+                //if (!is_polar_coordinate_system) {
+                //    plateau::geometry::GeoCoordinate refPoint;
+                //    plateau::geometry::ReferencePointFactory::GetReferencePoint(epsg, refPoint);
+                //    plateau::geometry::GeoReference geo_ref(geo_reference.getZoneID());
+                //    const auto& offset = geo_ref.project(refPoint);
+                //    //xyz = xyz + GeoReference::convertAxisToENU(geo_reference.getCoordinateSystem(), geo_reference.getReferencePoint());
+                //    xyz = lat_lon + offset;
+                //}
+
+                auto xyz = geo_reference.convert(lat_lon, false, epsg);
+
                 out_vertices.push_back(xyz);
             }
             assert(out_vertices.size() == vertices_lat_lon.size());
@@ -150,10 +162,17 @@ namespace plateau::polygonMesh {
                         if (!options.is_polar_coordinate_system) {
                             // 平面直角座標系の判定
                             plateau::geometry::GeoReference geo_ref(options.coordinate_zone_id, options.reference_point, options.unit_scale, options.mesh_axes);
-                            if (extent.contains(geo_ref.unproject(vertex))) {
+                            //plateau::geometry::GeoReference geo_ref(options.coordinate_zone_id);
+                            const auto unprojected = geo_ref.unproject(vertex);
+                            if (extent.contains(unprojected)) {
                                 is_in_extent = true;
                                 break;
                             }
+
+                            //Temp Debug =======================
+                            //is_in_extent = true;
+                            //break;
+                            //Temp Debug =======================
                         }
                         else {
                             if (extent.contains(vertex)) {
@@ -215,7 +234,7 @@ namespace plateau::polygonMesh {
             return;
 
         Mesh mesh;
-        cityGmlPolygonToMesh(polygon, gml_path, geo_reference_, mesh, options_.is_polar_coordinate_system);
+        cityGmlPolygonToMesh(polygon, gml_path, geo_reference_, mesh, options_.epsg_code);
 
         const auto from_axis = geometry::CoordinateSystem::ENU;
         const auto to_axis = options_.mesh_axes;
diff --git a/src/polygon_mesh/polygon_mesh_utils.cpp b/src/polygon_mesh/polygon_mesh_utils.cpp
@@ -1,6 +1,7 @@
 #include <plateau/polygon_mesh/polygon_mesh_utils.h>
 #include "plateau/polygon_mesh/mesh.h"
 #include "plateau/geometry/geo_reference.h"
+#include "plateau/geometry/geo_coordinate.h"
 #include "citygml/citymodel.h"
 #include <plateau/dataset/gml_file.h>
 
@@ -55,6 +56,16 @@ namespace plateau::polygonMesh {
         throw std::invalid_argument("Could not find position of CityObject.");
     }
 
+    //TVec3d PolygonMeshUtils::cityObjPos(const citygml::CityObject& city_obj, const double epsg, int coordinate_zone_id) {
+    //    const auto is_polar = plateau::geometry::ReferencePointFactory::IsPolarCoordinateSystem(epsg);
+    //    const auto pos = cityObjPos(city_obj);
+    //    if (is_polar) {
+    //        return pos;
+    //    }
+    //    plateau::geometry::GeoReference geo_ref(coordinate_zone_id);
+    //    return geo_ref.unproject(pos);
+    //}
+
     TVec3d PolygonMeshUtils::getCenterPoint(const CityModel& city_model, int coordinate_zone_id) {
         auto& envelope = city_model.getEnvelope();
         if (!envelope.validBounds()) {
diff --git a/wrappers/csharp/LibPLATEAU.NET/CSharpPLATEAU/PolygonMesh/MeshExtractOptions.cs b/wrappers/csharp/LibPLATEAU.NET/CSharpPLATEAU/PolygonMesh/MeshExtractOptions.cs
@@ -72,7 +72,7 @@ public static ConvertGranularity ToConvertGranularity(this MeshGranularity g)
     [StructLayout(LayoutKind.Sequential, CharSet = CharSet.Ansi)]
     public struct MeshExtractOptions
     {
-        public MeshExtractOptions(PlateauVector3d referencePoint, CoordinateSystem meshAxes, MeshGranularity meshGranularity, uint minLOD, uint maxLOD, bool exportAppearance, int gridCountOfSide, float unitScale, int coordinateZoneID, bool excludeCityObjectOutsideExtent, bool excludePolygonsOutsideExtent, bool enableTexturePacking, uint texturePackingResolution, bool attachMapTile, int mapTileZoomLevel, string mapTileURL, bool isPolarCoordinateSystem)
+        public MeshExtractOptions(PlateauVector3d referencePoint, CoordinateSystem meshAxes, MeshGranularity meshGranularity, uint minLOD, uint maxLOD, bool exportAppearance, int gridCountOfSide, float unitScale, int coordinateZoneID, bool excludeCityObjectOutsideExtent, bool excludePolygonsOutsideExtent, bool enableTexturePacking, uint texturePackingResolution, bool attachMapTile, int mapTileZoomLevel, string mapTileURL, bool isPolarCoordinateSystem, double epsgCode)
         {
             this.ReferencePoint = referencePoint;
             this.MeshAxes = meshAxes;
@@ -91,6 +91,7 @@ public MeshExtractOptions(PlateauVector3d referencePoint, CoordinateSystem meshA
             this.MapTileZoomLevel = mapTileZoomLevel;
             this.mapTileURL = mapTileURL;
             this.IsPolarCoordinateSystem = isPolarCoordinateSystem;
+            this.epsgCode = epsgCode;
 
             // 上で全てのメンバー変数を設定できてますが、バリデーションをするため念のためメソッドやプロパティも呼びます。
             SetLODRange(minLOD, maxLOD);
@@ -231,6 +232,8 @@ public string MapTileURL
         /// </summary>
         [MarshalAs(UnmanagedType.U1)] public bool IsPolarCoordinateSystem;
 
+        public double epsgCode;
+
         /// <summary> デフォルト値の設定を返します。 </summary>
         internal static MeshExtractOptions DefaultValue()
         {
