diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 05b9c9cfb5..32f3debde2 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -58,6 +58,7 @@ option(BUILD_PACKAGE "" OFF) 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) @@ -934,6 +935,17 @@ 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}") + message(STATUS "WITH_PROJ ${WITH_PROJ}") + + 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/ifcconvert/IfcConvert.cpp b/src/ifcconvert/IfcConvert.cpp index 059aebd45c..f73b137921 100644 --- a/src/ifcconvert/IfcConvert.cpp +++ b/src/ifcconvert/IfcConvert.cpp @@ -461,6 +461,7 @@ int main(int argc, char** argv) { ("space-name-transform", po::value(), "Additional transform to the space labels in SVG") ("edge-arrows", "Adds arrow heads to edge segments to signify edge direction") + ("ecef", "Write glTF in Earth-Centered Earth-Fixed coordinates. Requires PROJ") ; po::options_description cmdline_options; @@ -527,6 +528,7 @@ int main(int argc, char** argv) { const bool generate_uvs = vmap.count("generate-uvs") != 0; const bool validate = vmap.count("validate") != 0; const bool edge_arrows = vmap.count("edge-arrows") != 0; + const bool write_gltf_ecef = vmap.count("ecef") != 0; const bool no_wire_intersection_check = vmap.count("no-wire-intersection-check") != 0; const bool no_wire_intersection_tolerance = vmap.count("no-wire-intersection-tolerance") != 0; const bool strict_tolerance = vmap.count("strict-tolerance") != 0; @@ -840,7 +842,8 @@ int main(int argc, char** argv) { settings.set(SerializerSettings::USE_MATERIAL_NAMES, use_material_names); settings.set(SerializerSettings::USE_ELEMENT_TYPES, use_element_types); settings.set(SerializerSettings::USE_ELEMENT_HIERARCHY, use_element_hierarchy); - settings.set_deflection_tolerance(deflection_tolerance); + settings.set(SerializerSettings::WRITE_GLTF_ECEF, write_gltf_ecef); + settings.set_deflection_tolerance(deflection_tolerance); settings.set_angular_tolerance(angular_tolerance); settings.precision = precision; diff --git a/src/ifcgeom_schema_agnostic/GeometrySerializer.h b/src/ifcgeom_schema_agnostic/GeometrySerializer.h index 58fb7c4d64..63b1507d97 100644 --- a/src/ifcgeom_schema_agnostic/GeometrySerializer.h +++ b/src/ifcgeom_schema_agnostic/GeometrySerializer.h @@ -49,8 +49,11 @@ public: /// Use Y UP . /// Applicable for OBJ output. USE_Y_UP = 1ULL << (IfcGeom::IteratorSettings::NUM_SETTINGS + 7ULL), + /// Write in ECEF coordinates. + /// Applicable to glTF output + WRITE_GLTF_ECEF = 1ULL << (IfcGeom::IteratorSettings::NUM_SETTINGS + 8ULL), /// Number of different setting flags. - NUM_SETTINGS = 7 + NUM_SETTINGS = 8 }; SerializerSettings() diff --git a/src/serializers/GltfSerializer.cpp b/src/serializers/GltfSerializer.cpp index 3c285b7c03..f0a4f62737 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,13 +181,23 @@ void GltfSerializer::write(const IfcGeom::TriangulationElement* o) { node_array_.push_back(json_["nodes"].size()); const std::vector& m = o->transformation().matrix().data(); - // nb: note that this contains the Y-UP transform as well. - const std::array matrix_flat = { - m[0], m[ 2], -m[ 1], 0, - m[3], m[ 5], -m[ 4], 0, - m[6], m[ 8], -m[ 7], 0, - m[9], m[11], -m[10], 1 - }; + std::array matrix_flat; + if (settings_.get(SerializerSettings::WRITE_GLTF_ECEF)) { + matrix_flat = { + m[0], m[1], m[2], 0, + m[3], m[4], m[5], 0, + m[6], m[7], m[8], 0, + m[9], m[10], m[11], 1 + }; + } else { + // nb: note that this contains the Y-UP transform as well. + matrix_flat = { + m[0], m[2], -m[1], 0, + m[3], m[5], -m[4], 0, + m[6], m[8], -m[7], 0, + m[9], m[11], -m[10], 1 + }; + } static const std::array identity_matrix = {1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1}; json node; @@ -317,6 +331,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(); @@ -335,7 +365,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. @@ -375,4 +409,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 7450fa2afb..5f7265730f 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 IfcGeom::Material& style); public: GltfSerializer(const std::string& filename, const 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 4c7d1aa591..238d3df0d1 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