diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 7e281da5da..eaa5274e87 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -69,6 +69,7 @@ option(WITH_CGAL "Enable geometry interpretation using CGAL" ON) option(COLLADA_SUPPORT "Build IfcConvert with COLLADA support (requires OpenCOLLADA)." ON) option(GLTF_SUPPORT "Build IfcConvert with glTF support (requires json.hpp)." OFF) option(HDF5_SUPPORT "Enable HDF5 support (requires HDF5, zlib)" ON) +option(WITH_PROJ "Enable output of Earth-Centered Earth-Fixed glTF output using the PROJ library" OFF) option(IFCXML_SUPPORT "Build IfcParse with ifcXML support (requires libxml2)." ON) option(USD_SUPPORT "Build IfcConvert with USD support (requires pixar's USD library)." OFF) option(CITYJSON_SUPPORT "Build IfcConvert with CityJSON support (requires CityJSON library)." ON) @@ -965,6 +966,15 @@ if(BUILD_CONVERT OR BUILD_IFCPYTHON) add_library(Serializers ${SERIALIZERS_FILES}) set_target_properties(Serializers PROPERTIES COMPILE_FLAGS "-DIFC_GEOM_EXPORTS" VERSION "${PROJECT_VERSION}" SOVERSION "${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}") + if(WITH_PROJ) + target_compile_definitions(Serializers PRIVATE "WITH_PROJ") + if (PROJ_STATIC) + target_compile_definitions(Serializers PRIVATE "PROJ_DLL=") + endif() + target_include_directories(Serializers PRIVATE ${PROJ_INCLUDE_DIR} ${SQLITE_INCLUDE_DIR}) + target_link_libraries(Serializers ${PROJ_LIBRARIES}) + endif() + target_link_libraries(Serializers ${SERIALIZER_SCHEMA_LIBRARIES} ${OPENCOLLADA_LIBRARIES} ${USD_LIBRARIES}) endif(BUILD_CONVERT OR BUILD_IFCPYTHON) diff --git a/src/ifcgeom/GeometrySerializer.h b/src/ifcgeom/GeometrySerializer.h index 4c4fdded18..6a16193769 100644 --- a/src/ifcgeom/GeometrySerializer.h +++ b/src/ifcgeom/GeometrySerializer.h @@ -69,6 +69,12 @@ inline namespace settings { static constexpr bool defaultvalue = false; }; + struct WriteGltfEcef : public SettingBase { + static constexpr const char* const name = "ecef"; + static constexpr const char* const description = "Write glTF in Earth-Centered Earth-Fixed coordinates. Requires PROJ."; + static constexpr bool defaultvalue = false; + }; + struct FloatingPointDigits : public SettingBase { static constexpr const char* const name = "digits"; static constexpr const char* const description = "Sets the precision to be used to format floating-point values, 15 by default. " diff --git a/src/serializers/GltfSerializer.cpp b/src/serializers/GltfSerializer.cpp index 806f7867d6..d7bdf527c2 100644 --- a/src/serializers/GltfSerializer.cpp +++ b/src/serializers/GltfSerializer.cpp @@ -23,6 +23,10 @@ #include "../ifcparse/utils.h" +#ifdef WITH_PROJ +#include +#endif + #include static const uint32_t GLTF = 0x46546C67U; @@ -177,14 +181,26 @@ void GltfSerializer::write(const IfcGeom::TriangulationElement* o) { node_array_.push_back(json_["nodes"].size()); const auto& m = o->transformation().data()->ccomponents(); - // nb: note that this contains the Y-UP transform as well. + // @todo check - const std::array matrix_flat = { - m(0,0), m(2,0), -m(1,0), m(3,0), - m(0,1), m(2,1), -m(1,1), m(3,1), - m(0,2), m(2,2), -m(1,2), m(3,2), - m(0,3), m(2,3), -m(1,3), m(3,3) - }; + std::array matrix_flat; + if (settings_.get().get()) { + matrix_flat = { + m(0,0), m(1,0), m(2,0), m(3,0), + m(0,1), m(1,1), m(2,1), m(3,1), + m(0,2), m(1,2), m(2,2), m(3,2), + m(0,3), m(1,3), m(2,3), m(3,3) + }; + } else { + // nb: note that this contains the Y-UP transform as well. + matrix_flat = { + m(0,0), m(2,0), -m(1,0), m(3,0), + m(0,1), m(2,1), -m(1,1), m(3,1), + m(0,2), m(2,2), -m(1,2), m(3,2), + m(0,3), m(2,3), -m(1,3), m(3,3) + }; + } + static const std::array identity_matrix = {1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1}; json node; @@ -320,6 +336,22 @@ void write_block(std::ostream& fs, It begin, It end) { } void GltfSerializer::finalize() { + if (north_rotation_) { + (*north_rotation_)["children"] = json::array(); + for (int i = 0; i < json_["nodes"].size(); ++i) { + (*north_rotation_)["children"].push_back(i); + } + json_["nodes"].push_back(*north_rotation_); + } + + if (ecef_transform_) { + (*ecef_transform_)["children"] = json::array(); + for (int i = 0; i < json_["nodes"].size(); ++i) { + (*ecef_transform_)["children"].push_back(i); + } + json_["nodes"].push_back(*ecef_transform_); + } + tmp_fstream1_.close(); tmp_fstream2_.close(); @@ -338,7 +370,11 @@ void GltfSerializer::finalize() { } json scene_0; - scene_0["nodes"] = node_array_; + if (north_rotation_ || ecef_transform_) { + scene_0["nodes"] = std::array{json_["nodes"].size() - 1}; + } else { + scene_0["nodes"] = node_array_; + } json_["scenes"].push_back(scene_0); //The generated glb file will contain the indices buffer followed by the vertices buffer. @@ -378,4 +414,255 @@ void GltfSerializer::finalize() { write_padding(fstream_, binary_length); } +namespace { + void normalize(std::array& v) { + auto l = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + + v[0] /= l; + v[1] /= l; + v[2] /= l; + } + + void cross(const std::array& v1, const std::array& v2, std::array& result) { + result[0] = v1[1] * v2[2] - v1[2] * v2[1]; + result[1] = v1[2] * v2[0] - v1[0] * v2[2]; + result[2] = v1[0] * v2[1] - v1[1] * v2[0]; + } + + void proj_log(void *, int, const char* c) { + Logger::Error("PROJ: " + std::string(c)); + } +} + +void GltfSerializer::setFile(IfcParse::IfcFile* f) { + if (!settings_.get(SerializerSettings::WRITE_GLTF_ECEF)) { + return; + } + + boost::optional crs_epsg; + boost::optional> crs_x_axis; + boost::optional> eastings_northings_elevation; + + aggregate_of_instance::ptr coordops; + try { + coordops = f->instances_by_type("IfcCoordinateOperation"); + } catch (IfcParse::IfcException&) { + // Ignored. Schema likely doesn't support IfcCoordinateOperation. + } + if (coordops) { + for (auto& coordop : *coordops) { + IfcUtil::IfcBaseClass* source_crs = *coordop->as()->get("SourceCRS"); + if (source_crs->declaration().is("IfcGeometricRepresentationContext")) { + IfcUtil::IfcBaseClass* target_crs = *coordop->as()->get("TargetCRS"); + auto name_attr = target_crs->as()->get("Name"); + if (coordop->declaration().is("IfcMapConversion")) { + + if (!name_attr->isNull()) { + std::string epsg_code = *name_attr; + crs_epsg = epsg_code; + + // @todo in which unit are these? + double eastings = *coordop->as()->get("Eastings"); + double northings = *coordop->as()->get("Northings"); + double height = *coordop->as()->get("OrthogonalHeight"); + height = 0.; + + eastings_northings_elevation = { { eastings, northings, height} }; + + auto xaxis_attr = coordop->as()->get("XAxisAbscissa"); + auto yaxis_attr = coordop->as()->get("XAxisOrdinate"); + if (!xaxis_attr->isNull() && !yaxis_attr->isNull()) { + double xaxis = *xaxis_attr; + double yaxis = *yaxis_attr; + + crs_x_axis = { { xaxis, yaxis, 0. } }; + } + } + } + } + } + } + + if (!crs_epsg) { + auto sites = f->instances_by_type("IfcSite"); + + if (sites && sites->size() == 1) { + auto lat_attr = (*sites->begin())->as()->get("RefLatitude"); + auto lon_attr = (*sites->begin())->as()->get("RefLongitude"); + + if (!lat_attr->isNull() && !lon_attr->isNull()) { + std::vector lat_dms = *lat_attr; + std::vector lon_dms = *lon_attr; + + auto to_decimal = [](const std::vector& dms) { + double val = dms[0] + dms[1] / 60. + dms[2] / 3600.; + if (dms.size() == 4) { + val += dms[3] / 3600.e6; + } + return val; + }; + + auto lat = to_decimal(lat_dms); + auto lon = to_decimal(lon_dms); + double elev = 0.; + + /* + auto elev_attr = (*sites->begin())->as()->get("RefElevation"); + if (!elev_attr->isNull()) { + elev = *elev_attr; + } + */ + + crs_epsg.reset("EPSG:4326"); + eastings_northings_elevation = { { lat, lon, elev } }; + } + } + } + + auto contexts = f->instances_by_type_excl_subtypes("IfcGeometricRepresentationContext"); + + if (contexts && contexts->size() > 0) { + auto context = (*contexts->begin())->as(); + auto north_attr = context->get("TrueNorth"); + if (!north_attr->isNull()) { + IfcUtil::IfcBaseClass* north = *north_attr; + if (north->declaration().is("IfcDirection")) { + std::vector ratios = *north->as()->get("DirectionRatios"); + crs_x_axis = { { ratios[1], -ratios[0], 0. } }; + } + } + } + +#ifdef WITH_PROJ + + if (crs_epsg) { + PJ_COORD wgs84_point; + + auto C = proj_context_create(); + proj_log_func(C, nullptr, proj_log); + + // @todo a bit ugly we assume a proj.db in current working directory. + // a very simplistic but at least portable solution. + proj_context_set_database_path(C, "proj.db", nullptr, nullptr); + + if (*crs_epsg == "EPSG:4326") { + wgs84_point = proj_coord( + (*eastings_northings_elevation)[0], + (*eastings_northings_elevation)[1], + (*eastings_northings_elevation)[2], + 0); + } else { + // @todo a bit ugly we assume a proj.db in current working directory. + // a very simplistic but at least portable solution. + proj_context_set_database_path(C, "proj.db", nullptr, nullptr); + + auto P = proj_create_crs_to_crs( + C, crs_epsg->c_str(), "EPSG:4326", + NULL); + + if (!P) { + Logger::Error("Failed to create PROJ transformation object"); + return; + } + + auto a = proj_coord( + (*eastings_northings_elevation)[0], + (*eastings_northings_elevation)[1], + (*eastings_northings_elevation)[2], + 0); + + wgs84_point = proj_trans(P, PJ_FWD, a); + + Logger::Notice("Calculated latitude: " + std::to_string(wgs84_point.lp.lam) + " longitude: " + std::to_string(wgs84_point.lp.phi)); + } + + std::swap(wgs84_point.lp.phi, wgs84_point.lp.lam); + + const char *input_crs = "+proj=latlong +datum=WGS84"; + const char *output_crs = "+proj=geocent +datum=WGS84 +units=m"; + + // Create a transformation object + PJ *transform = proj_create_crs_to_crs(C, input_crs, output_crs, NULL); + + // Perform the transformation + PJ_COORD output_point = proj_trans(transform, PJ_FWD, wgs84_point); + + // Extract the ECEF coordinates + double x = output_point.xyz.x; + double y = output_point.xyz.y; + double z = output_point.xyz.z; + + const char *ellipsoid_def = "WGS84"; + + // Create a CRS object representing the ellipsoid + PJ *ellipsoid_crs = proj_create(C, ellipsoid_def); + + if (!ellipsoid_crs) { + Logger::Error("Failed to create ellipsoid CRS"); + return; + } + + auto ellipse = proj_get_ellipsoid(C, ellipsoid_crs); + + + int _; + double semi_major, semi_minor, __; + proj_ellipsoid_get_parameters(C, ellipse, &semi_major, &semi_minor, &_, &__); + + std::array dxyz = { { + x * (1. / (semi_major * semi_major)), + y * (1. / (semi_major * semi_major)), + z * (1. / (semi_minor * semi_minor)) + } }; + normalize(dxyz); + + // Oblate spheroid, so X and Y axis are equal, so rotation around Z yields east axis. + std::array east_xyz = { { + -y, + x, + 0. + } }; + normalize(east_xyz); + + std::array north; + cross(dxyz, east_xyz, north); + + std::array matrix = { + east_xyz[0], east_xyz[1], east_xyz[2], 0, + north[0], north[1], north[2], 0., + dxyz[0], dxyz[1], dxyz[2], 0, + 0,0,0,1 + }; + + ecef_transform_ = json::object({ + {"matrix", matrix } + }); + + json_["extensions"]["CESIUM_RTC"]["center"] = std::array{ {x, y, z} }; + json_["extensionsUsed"].push_back("CESIUM_RTC"); + + // Clean up + proj_destroy(ellipsoid_crs); + proj_destroy(transform); + proj_context_destroy(C); + } + + if (crs_x_axis) { + normalize(*crs_x_axis); + + auto phi = std::atan2((*crs_x_axis)[1], (*crs_x_axis)[0]); + + north_rotation_ = json::object({ + {"matrix", std::array{ + +std::cos(-phi), -std::sin(-phi), 0., 0., + +std::sin(-phi), +std::cos(-phi), 0., 0., + 0., 0., 1., 0., + 0., 0., 0., 1. + }} + }); + } +#endif +} + + #endif diff --git a/src/serializers/GltfSerializer.h b/src/serializers/GltfSerializer.h index 33b26eb696..cbdbe0dacb 100644 --- a/src/serializers/GltfSerializer.h +++ b/src/serializers/GltfSerializer.h @@ -36,9 +36,9 @@ private: std::ofstream fstream_, tmp_fstream1_, tmp_fstream2_; std::map materials_, meshes_; json json_, node_array_; + boost::optional ecef_transform_, north_rotation_; int bufferViewId; - int writeMaterial(const ifcopenshell::geometry::taxonomy::style& style); public: GltfSerializer(const std::string& filename, const ifcopenshell::geometry::Settings& geometry_settings, const ifcopenshell::geometry::SerializerSettings& settings); @@ -50,7 +50,7 @@ public: void finalize(); bool isTesselated() const { return true; } void setUnitNameAndMagnitude(const std::string& /*name*/, float /*magnitude*/) {} - void setFile(IfcParse::IfcFile*) {} + void setFile(IfcParse::IfcFile*); }; #endif diff --git a/win/build-deps.cmd b/win/build-deps.cmd index 9d5ac6f07a..a9dff10cfb 100644 --- a/win/build-deps.cmd +++ b/win/build-deps.cmd @@ -193,6 +193,41 @@ IF "%IFCOS_INSTALL_PYTHON%"=="TRUE" ( echo PYTHONHOME=%PYTHONHOME%>>"%~dp0\%BUILD_DEPS_CACHE_PATH%" ) + +:proj + +set DEPENDENCY_NAME=sqlite3 +md %INSTALL_DIR%\sqlite3\lib %INSTALL_DIR%\sqlite3\bin %INSTALL_DIR%\sqlite3\include +call :DownloadFile https://www.sqlite.org/2023/sqlite-amalgamation-3430100.zip "%DEPS_DIR%" sqlite-amalgamation-3430100.zip +IF NOT %ERRORLEVEL%==0 GOTO :Error +call :ExtractArchive sqlite-amalgamation-3430100.zip "%DEPS_DIR%" "%DEPS_DIR%\sqlite-amalgamation-3430100" +IF NOT %ERRORLEVEL%==0 GOTO :Error +pushd "%DEPS_DIR%\sqlite-amalgamation-3430100" +cl /c sqlite3.c +lib /OUT:%INSTALL_DIR%\sqlite3\lib\sqlite3.lib sqlite3.obj +cl sqlite3.c shell.c /link /out:%INSTALL_DIR%\sqlite3\bin\sqlite3.exe +copy sqlite3.h %INSTALL_DIR%\sqlite3\include +popd + +set DEPENDENCY_NAME=proj +set DEPENDENCY_DIR=%DEPS_DIR%\proj-9.2.1 +call :DownloadFile https://download.osgeo.org/proj/proj-9.2.1.zip "%DEPS_DIR%" proj-9.2.1.zip +IF NOT %ERRORLEVEL%==0 GOTO :Error +call :ExtractArchive proj-9.2.1.zip "%DEPS_DIR%" "%DEPS_DIR%\proj-9.2.1" +IF NOT %ERRORLEVEL%==0 GOTO :Error +cd "%DEPENDENCY_DIR%" +call :RunCMake -DCMAKE_INSTALL_PREFIX="%INSTALL_DIR%\proj-9.2.1" ^ + -DSQLITE3_INCLUDE_DIR=%INSTALL_DIR%\sqlite3\include -DSQLITE3_LIBRARY=%INSTALL_DIR%\sqlite3\lib\sqlite3.lib ^ + -DENABLE_TIFF=Off -DENABLE_CURL=Off -DBUILD_PROJSYNC=Off ^ + -DBUILD_SHARED_LIBS=Off +IF NOT %ERRORLEVEL%==0 GOTO :Error +call :BuildSolution "%DEPENDENCY_DIR%\%BUILD_DIR%\PROJ.sln" %BUILD_CFG% +IF NOT %ERRORLEVEL%==0 GOTO :Error +call :InstallCMakeProject "%DEPENDENCY_DIR%\%BUILD_DIR%" %BUILD_CFG% +IF NOT %ERRORLEVEL%==0 GOTO :Error + +goto :Successful + :mpir set DEPENDENCY_NAME=mpir set DEPENDENCY_DIR=%DEPS_DIR%\mpir