From dad30b0cc212e6ecd4b332cbc4f1335349f89188 Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski Date: Sat, 8 Oct 2016 08:44:44 -0400 Subject: [PATCH] alternative terrain engine - SGMesh utilizing pagedLOD --- CMakeLists.txt | 25 +- simgear/environment/CMakeLists.txt | 5 +- simgear/scene/CMakeLists.txt | 3 + simgear/scene/dem/CMakeLists.txt | 24 + simgear/scene/dem/ReaderWriterPGT.cxx | 418 +++++++++++ simgear/scene/dem/ReaderWriterPGT.hxx | 53 ++ simgear/scene/dem/SGDem.cxx | 332 +++++++++ simgear/scene/dem/SGDem.hxx | 82 +++ simgear/scene/dem/SGDemLevel.cxx | 72 ++ simgear/scene/dem/SGDemLevel.hxx | 90 +++ simgear/scene/dem/SGDemRoot.cxx | 118 +++ simgear/scene/dem/SGDemRoot.hxx | 102 +++ simgear/scene/dem/SGDemSession.cxx | 79 ++ simgear/scene/dem/SGDemSession.hxx | 78 ++ simgear/scene/dem/SGDemTile.cxx | 361 ++++++++++ simgear/scene/dem/SGDemTile.hxx | 65 ++ simgear/scene/dem/SGDemTile_gdal.cxx | 628 ++++++++++++++++ simgear/scene/dem/SGMesh.cxx | 414 +++++++++++ simgear/scene/dem/SGMesh.hxx | 718 +++++++++++++++++++ simgear/scene/dem/TriMesh.h | 257 +++++++ simgear/scene/dem/TriMesh_curvature.cc | 331 +++++++++ simgear/scene/dem/TriMesh_normals.cc | 116 +++ simgear/scene/tgdb/BucketBox.hxx | 43 +- simgear/scene/tgdb/ReaderWriterSPT.cxx | 24 +- simgear/scene/tgdb/userdata.cxx | 4 + simgear/scene/util/SGReaderWriterOptions.hxx | 18 + simgear/simgear_config_cmake.h.in | 3 +- 27 files changed, 4443 insertions(+), 20 deletions(-) create mode 100644 simgear/scene/dem/CMakeLists.txt create mode 100644 simgear/scene/dem/ReaderWriterPGT.cxx create mode 100644 simgear/scene/dem/ReaderWriterPGT.hxx create mode 100644 simgear/scene/dem/SGDem.cxx create mode 100644 simgear/scene/dem/SGDem.hxx create mode 100644 simgear/scene/dem/SGDemLevel.cxx create mode 100644 simgear/scene/dem/SGDemLevel.hxx create mode 100644 simgear/scene/dem/SGDemRoot.cxx create mode 100644 simgear/scene/dem/SGDemRoot.hxx create mode 100644 simgear/scene/dem/SGDemSession.cxx create mode 100644 simgear/scene/dem/SGDemSession.hxx create mode 100644 simgear/scene/dem/SGDemTile.cxx create mode 100644 simgear/scene/dem/SGDemTile.hxx create mode 100644 simgear/scene/dem/SGDemTile_gdal.cxx create mode 100644 simgear/scene/dem/SGMesh.cxx create mode 100644 simgear/scene/dem/SGMesh.hxx create mode 100644 simgear/scene/dem/TriMesh.h create mode 100644 simgear/scene/dem/TriMesh_curvature.cc create mode 100644 simgear/scene/dem/TriMesh_normals.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index 803d0b9e..44ff9a13 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -120,12 +120,14 @@ endif() option(SIMGEAR_HEADLESS "Set to ON to build SimGear without GUI/graphics support" OFF) option(ENABLE_RTI "Set to ON to build SimGear with RTI support" OFF) +option(ENABLE_GDAL "Set to ON to build SimGear with GDAL support" ON) option(ENABLE_TESTS "Set to OFF to disable building SimGear's test applications" ON) option(ENABLE_SOUND "Set to OFF to disable building SimGear's sound support" ON) option(USE_AEONWAVE "Set to ON to use AeonWave instead of OpenAL" OFF) option(ENABLE_PKGUTIL "Set to ON to build the sg_pkgutil application (default)" ON) option(ENABLE_DNS "Set to ON to use udns library and DNS service resolver" ON) option(ENABLE_SIMD "Enable SSE/SSE2 support for x86 compilers" ON) +option(ENABLE_OPENMP "Enable OpenMP compiler support" ON) include (DetectArch) @@ -269,6 +271,13 @@ else() message(STATUS "RTI: DISABLED") endif(ENABLE_RTI) +if(ENABLE_GDAL) + find_package(GDAL 2.0.0 REQUIRED) + if (GDAL_FOUND) + include_directories(${GDAL_INCLUDE_DIR}) + endif(GDAL_FOUND) +endif(ENABLE_GDAL) + check_function_exists(gettimeofday HAVE_GETTIMEOFDAY) check_function_exists(rint HAVE_RINT) check_function_exists(mkdtemp HAVE_MKDTEMP) @@ -381,6 +390,19 @@ if (CLANG) endif() endif() +if (ENABLE_OPENMP) + find_package(OpenMP) + if(OPENMP_FOUND) + message(STATUS "OpenMP: ENABLED") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + else() + message(STATUS "OpenMP: NOT FOUND") + endif() +else() + message(STATUS "OpenMP: DISABLED") +endif() + if(CMAKE_COMPILER_IS_GNUCXX OR CMAKE_CXX_COMPILER_ID STREQUAL "Clang") # boost goes haywire wrt static asserts check_cxx_compiler_flag(-Wno-unused-local-typedefs HAS_NOWARN_UNUSED_TYPEDEFS) @@ -464,7 +486,8 @@ set(TEST_LIBS_INTERNAL_CORE ${RT_LIBRARY} ${DL_LIBRARY} ${COCOA_LIBRARY} - ${CURL_LIBRARIES}) + ${CURL_LIBRARIES} + ${GDAL_LIBRARY}) set(TEST_LIBS SimGearCore ${TEST_LIBS_INTERNAL_CORE}) if(NOT SIMGEAR_HEADLESS) diff --git a/simgear/environment/CMakeLists.txt b/simgear/environment/CMakeLists.txt index d5255b04..4c499af6 100644 --- a/simgear/environment/CMakeLists.txt +++ b/simgear/environment/CMakeLists.txt @@ -10,13 +10,14 @@ if(ENABLE_TESTS) add_executable(test_metar test_metar.cxx) if (SIMGEAR_SHARED) - target_link_libraries(test_metar SimGearScene) + target_link_libraries(test_metar SimGearScene ${GDAL_LIBRARY}) else() target_link_libraries(test_metar SimGearScene SimGearCore ${CMAKE_THREAD_LIBS_INIT} ${ZLIB_LIBRARY} - ${RT_LIBRARY}) + ${RT_LIBRARY} + ${GDAL_LIBRARY}) endif() add_test(metar ${EXECUTABLE_OUTPUT_PATH}/test_metar) diff --git a/simgear/scene/CMakeLists.txt b/simgear/scene/CMakeLists.txt index 73429d55..ccbd9416 100644 --- a/simgear/scene/CMakeLists.txt +++ b/simgear/scene/CMakeLists.txt @@ -15,3 +15,6 @@ foreach( mylibfolder endforeach( mylibfolder ) +if(ENABLE_GDAL) + add_subdirectory(dem) +endif(ENABLE_GDAL) diff --git a/simgear/scene/dem/CMakeLists.txt b/simgear/scene/dem/CMakeLists.txt new file mode 100644 index 00000000..981a8afa --- /dev/null +++ b/simgear/scene/dem/CMakeLists.txt @@ -0,0 +1,24 @@ +include (SimGearComponent) + +set(HEADERS + ReaderWriterPGT.hxx + SGDem.hxx + SGDemLevel.hxx + SGDemRoot.hxx + SGDemSession.hxx + SGDemTile.hxx + SGMesh.hxx +) + +set(SOURCES + ReaderWriterPGT.cxx + SGDem.cxx + SGDemLevel.cxx + SGDemRoot.cxx + SGDemSession.cxx + SGDemTile.cxx + SGDemTile_gdal.cxx + SGMesh.cxx +) + +simgear_scene_component(dem scene/dem "${SOURCES}" "${HEADERS}") diff --git a/simgear/scene/dem/ReaderWriterPGT.cxx b/simgear/scene/dem/ReaderWriterPGT.cxx new file mode 100644 index 00000000..e6f94943 --- /dev/null +++ b/simgear/scene/dem/ReaderWriterPGT.cxx @@ -0,0 +1,418 @@ +// ReaderWriterPGT.cxx -- Provide a paged database for flightgear scenery. +// +// Copyright (C) 2010 - 2013 Mathias Froehlich +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include +#include "ReaderWriterPGT.hxx" + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +namespace simgear { + +// Cull away tiles that we watch from downside +struct ReaderWriterPGT::CullCallback : public osg::NodeCallback { + virtual ~CullCallback() + { } + virtual void operator()(osg::Node* node, osg::NodeVisitor* nv) + { + const osg::BoundingSphere& nodeBound = node->getBound(); + // If the bounding sphere of the node is empty, there is nothing to do + if (!nodeBound.valid()) + return; + + // Culling away tiles that we look at from the downside. + // This is done by computing the maximum distance we can + // see something from the current eyepoint. If the sphere + // that is defined by this radius does no intersects the + // nodes sphere, then this tile is culled away. + // Computing this radius happens by two rectangular triangles: + // Let r be the view point. rmin is the minimum radius we find + // a ground surface we need to look above. rmax is the + // maximum object radius we expect any object. + // + // d1 d2 + // x----x----x + // r\ rmin /rmax + // \ | / + // \ | / + // \|/ + // + // The distance from the eyepoint to the point + // where the line of sight is perpandicular to + // the radius vector with minimal height is + // d1 = sqrt(r^2 - rmin^2). + // The distance from the point where the line of sight + // is perpandicular to the radius vector with minimal height + // to the highest possible object on earth with radius rmax is + // d2 = sqrt(rmax^2 - rmin^2). + // So the maximum distance we can see something on the earth + // from a viewpoint r is + // d = d1 + d2 + + // This is the equatorial earth radius minus 450m, + // little lower than Dead Sea. + float rmin = 6378137 - 450; + float rmin2 = rmin*rmin; + // This is the equatorial earth radius plus 9000m, + // little higher than Mount Everest. + float rmax = 6378137 + 9000; + float rmax2 = rmax*rmax; + + // Check if we are looking from below any ground + osg::Vec3 viewPoint = nv->getViewPoint(); + // blow the viewpoint up to a spherical earth with equatorial radius: + osg::Vec3 sphericViewPoint = viewPoint; + sphericViewPoint[2] *= 1.0033641; + float r2 = sphericViewPoint.length2(); + if (r2 <= rmin2) + return; + + // Due to this line of sight computation, the visible tiles + // are limited to be within a sphere with radius d1 + d2. + float d1 = sqrtf(r2 - rmin2); + float d2 = sqrtf(rmax2 - rmin2); + // Note that we again base the sphere around elliptic view point, + // but use the radius from the spherical computation. + if (!nodeBound.intersects(osg::BoundingSphere(viewPoint, d1 + d2))) + return; + + traverse(node, nv); + } +}; + +struct ReaderWriterPGT::LocalOptions { + LocalOptions(const osgDB::Options* options) : + _options(options) + { + osg::ref_ptr sgOptions; + sgOptions = SGReaderWriterOptions::copyOrCreate(options); + + std::string pageLevelsString; + std::string meshResolutionString; + std::string meshTexturing; + + if (_options) { + pageLevelsString = _options->getPluginStringData("SimGear::SPT_PAGE_LEVELS"); + meshResolutionString = _options->getPluginStringData("SimGear::SPT_MESH_RESOLUTION"); + meshTexturing = _options->getPluginStringData("SimGear::SPT_LOD_TEXTURING"); + } + + // Get the default if nothing given from outside + // ignore option - pagelevels come from mesh + _pageLevels.push_back(1); + _pageLevels.push_back(2); + + // dem level 2 - 2, 4, 12 degrees + _pageLevels.push_back(3); + _pageLevels.push_back(4); + _pageLevels.push_back(5); + + // dem level 1 - 1/8, 1/4, 1/2 and 1 degree + _pageLevels.push_back(6); + _pageLevels.push_back(7); + _pageLevels.push_back(8); + _pageLevels.push_back(9); + + _dem = sgOptions->getDem(); + + // Get the default if nothing given from outside + if (meshResolutionString.empty()) { + _meshResolution = 1; + } else { + // If configured from outside + std::stringstream ss(meshResolutionString); + while (ss.good()) { + ss >> _meshResolution; + } + } + + if ( meshTexturing.empty() ) { + _textureMethod = SGMesh::TEXTURE_BLUEMARBLE; + } else if ( meshTexturing == "bluemarble" ) { + _textureMethod = SGMesh::TEXTURE_BLUEMARBLE; + } else if ( meshTexturing == "raster" ) { + _textureMethod = SGMesh::TEXTURE_RASTER; + } else if ( meshTexturing == "debug" ) { + _textureMethod = SGMesh::TEXTURE_DEBUG; + } + } + + bool isPageLevel(unsigned level) const + { + return std::find(_pageLevels.begin(), _pageLevels.end(), level) != _pageLevels.end(); + } + + std::string getLodPathForBucketBox(const BucketBox& bucketBox) const + { + std::stringstream ss; + ss << "LOD/"; + for (std::vector::const_iterator i = _pageLevels.begin(); i != _pageLevels.end(); ++i) { + if (bucketBox.getStartLevel() <= *i) + break; + ss << bucketBox.getParentBox(*i) << "/"; + } + ss << bucketBox; + return ss.str(); + } + + float getRangeMultiplier() const + { + float rangeMultiplier = 2; + if (!_options) + return rangeMultiplier; + std::stringstream ss(_options->getPluginStringData("SimGear::SPT_RANGE_MULTIPLIER")); + ss >> rangeMultiplier; + return rangeMultiplier; + } + + const osgDB::Options* _options; + std::vector _pageLevels; + + SGDemPtr _dem; + unsigned _meshResolution; + SGMesh::TextureMethod _textureMethod; +}; + +ReaderWriterPGT::ReaderWriterPGT() +{ + supportsExtension("pgt", "SimGear realtime paged terrain meta database."); +} + +ReaderWriterPGT::~ReaderWriterPGT() +{ +} + +const char* +ReaderWriterPGT::className() const +{ + return "simgear::ReaderWriterPGT"; +} + +osgDB::ReaderWriter::ReadResult +ReaderWriterPGT::readObject(const std::string& fileName, const osgDB::Options* options) const +{ + // We get called with different extensions. To make sure search continues, + // we need to return FILE_NOT_HANDLED in this case. + if (osgDB::getLowerCaseFileExtension(fileName) != "pgt") + return ReadResult(ReadResult::FILE_NOT_HANDLED); + if (fileName != "state.pgt") + return ReadResult(ReadResult::FILE_NOT_FOUND); + + osg::StateSet* stateSet = new osg::StateSet; + stateSet->setAttributeAndModes(new osg::CullFace); + + std::string imageFileName = options->getPluginStringData("SimGear::FG_WORLD_TEXTURE"); + if (imageFileName.empty()) { + imageFileName = options->getPluginStringData("SimGear::FG_ROOT"); + imageFileName = osgDB::concatPaths(imageFileName, "Textures"); + imageFileName = osgDB::concatPaths(imageFileName, "Globe"); + imageFileName = osgDB::concatPaths(imageFileName, "world.topo.bathy.200407.3x4096x2048.png"); + } + if (osg::Image* image = osgDB::readImageFile(imageFileName, options)) { + osg::Texture2D* texture = new osg::Texture2D; + texture->setImage(image); + texture->setWrap(osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT); + texture->setWrap(osg::Texture2D::WRAP_T, osg::Texture2D::CLAMP); + stateSet->setTextureAttributeAndModes(0, texture); + } + + return stateSet; +} + +osgDB::ReaderWriter::ReadResult +ReaderWriterPGT::readNode(const std::string& fileName, const osgDB::Options* options) const +{ + LocalOptions localOptions(options); + SG_LOG(SG_IO, SG_WARN, "ReaderWriterPGT::readNode - reading:" << fileName ); + + // The file name without path and without the pgt extension + std::string strippedFileName = osgDB::getStrippedName(fileName); + if (strippedFileName == "earth") + return ReadResult(createTree(BucketBox(-180, -90, 360, 180), localOptions, true)); + + std::stringstream ss(strippedFileName); + BucketBox bucketBox; + ss >> bucketBox; + if (ss.fail()) { + SG_LOG(SG_IO, SG_WARN, "error reading:" << strippedFileName ); + return ReadResult::FILE_NOT_FOUND; + } + + BucketBox bucketBoxList[2]; + unsigned bucketBoxListSize = bucketBox.periodicSplit(bucketBoxList); + if (bucketBoxListSize == 0) + return ReadResult::FILE_NOT_FOUND; + + if (bucketBoxListSize == 1) + return ReadResult(createTree(bucketBoxList[0], localOptions, true)); + + assert(bucketBoxListSize == 2); + osg::ref_ptr group = new osg::Group; + group->addChild(createTree(bucketBoxList[0], localOptions, true)); + group->addChild(createTree(bucketBoxList[1], localOptions, true)); + return ReadResult(group); +} + +osg::ref_ptr +ReaderWriterPGT::createTree(const BucketBox& bucketBox, const LocalOptions& options, bool topLevel) const +{ + if (bucketBox.getIsBucketSize()) { + std::string fileName; + fileName = bucketBox.getBucket().gen_index_str() + std::string(".stg"); + + return createTileMesh(bucketBox, options._options); + } else if (!topLevel && options.isPageLevel(bucketBox.getStartLevel())) { + return createPagedLOD(bucketBox, options); + } else { + BucketBox bucketBoxList[100]; + unsigned numTiles = bucketBox.getSubDivision(bucketBoxList, 100); + if (numTiles == 0) + return 0; + + if (numTiles == 1) + return createTree(bucketBoxList[0], options, false); + + osg::ref_ptr group = new osg::Group; + for (unsigned i = 0; i < numTiles; ++i) { + osg::ref_ptr node = createTree(bucketBoxList[i], options, false); + if (!node.valid()) + continue; + group->addChild(node.get()); + } + if (!group->getNumChildren()) + return 0; + + return group; + } +} + +osg::ref_ptr +ReaderWriterPGT::createPagedLOD(const BucketBox& bucketBox, const LocalOptions& options) const +{ + osg::PagedLOD* pagedLOD = new osg::PagedLOD; + + pagedLOD->setCenterMode(osg::PagedLOD::USER_DEFINED_CENTER); + SGSpheref sphere = bucketBox.getBoundingSphere(); + pagedLOD->setCenter(toOsg(sphere.getCenter())); + pagedLOD->setRadius(sphere.getRadius()); + + pagedLOD->setCullCallback(new CullCallback); + + osg::ref_ptr localOptions; + localOptions = static_cast(options._options->clone(osg::CopyOp())); + // FIXME: + // The particle systems have nodes with culling disabled. + // PagedLOD nodes with childnodes like this will never expire. + // So, for now switch them off. + localOptions->setPluginStringData("SimGear::PARTICLESYSTEM", "OFF"); + pagedLOD->setDatabaseOptions(localOptions.get()); + + // The break point for the low level of detail to the high level of detail + float rangeMultiplier = options.getRangeMultiplier(); + float range = rangeMultiplier*sphere.getRadius(); + + // Look for a low level of detail tile + std::string lodPath = options.getLodPathForBucketBox(bucketBox); + const char* extensions[] = { ".btg.gz", ".flt" }; + for (unsigned i = 0; i < sizeof(extensions)/sizeof(extensions[0]); ++i) { + std::string fileName = osgDB::findDataFile(lodPath + extensions[i], options._options); + if (fileName.empty()) + continue; + osg::ref_ptr node = osgDB::readRefNodeFile(fileName, options._options); + if (!node.valid()) + continue; + pagedLOD->addChild(node.get(), range, std::numeric_limits::max()); + break; + } + // Add the static sea level textured shell if there is nothing found + if (pagedLOD->getNumChildren() == 0) { + osg::ref_ptr node = createTileMesh(bucketBox, options._options); + if (node.valid()) + pagedLOD->addChild(node.get(), range, std::numeric_limits::max()); + } + + // Add the paged file name that creates the subtrees on demand + std::stringstream ss; + ss << bucketBox << ".pgt"; + pagedLOD->setFileName(pagedLOD->getNumChildren(), ss.str()); + pagedLOD->setRange(pagedLOD->getNumChildren(), 0.0, range); + + return pagedLOD; +} + +osg::ref_ptr +ReaderWriterPGT::createTileMesh(const BucketBox& bucketBox, const LocalOptions& options) const +{ + if (options._options->getPluginStringData("SimGear::FG_EARTH") != "ON") + return 0; + + SGSpheref sphere = bucketBox.getBoundingSphere(); + osg::Matrixd transform; + transform.makeTranslate(toOsg(-sphere.getCenter())); + + // TODO : return geode, not geometry - so we texture in SGMesh + // osg::Geometry* geometry = bucketBox.getTileTriangleMesh( options._dem, options._meshResolution, options._textureMethod ); + osg::Geode* geode = bucketBox.getTileTriangleMesh( options._dem, options._meshResolution, options._textureMethod, options._options ); + if ( geode ) { + transform.makeTranslate(toOsg(sphere.getCenter())); + osg::MatrixTransform* matrixTransform = new osg::MatrixTransform(transform); + matrixTransform->setDataVariance(osg::Object::STATIC); + matrixTransform->addChild(geode); + + return matrixTransform; + } else { + return 0; + } +} + +osg::ref_ptr +ReaderWriterPGT::getLowLODStateSet(const LocalOptions& options) const +{ + osg::ref_ptr localOptions; + localOptions = static_cast(options._options->clone(osg::CopyOp())); + localOptions->setObjectCacheHint(osgDB::Options::CACHE_ALL); + + osg::ref_ptr object = osgDB::readRefObjectFile("state.pgt", localOptions.get()); + if (!dynamic_cast(object.get())) + return 0; + + return static_cast(object.get()); +} + +} // namespace simgear + +// simgear::ModelRegistryCallbackProxy g_pgtCallbackProxy("pgt"); diff --git a/simgear/scene/dem/ReaderWriterPGT.hxx b/simgear/scene/dem/ReaderWriterPGT.hxx new file mode 100644 index 00000000..e7783f43 --- /dev/null +++ b/simgear/scene/dem/ReaderWriterPGT.hxx @@ -0,0 +1,53 @@ +// ReaderWriterSPT.cxx -- Provide a paged database for flightgear scenery. +// +// Copyright (C) 2010 - 2013 Mathias Froehlich +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// + +#ifndef _READERWRITERPGT_HXX +#define _READERWRITERPGT_HXX + +#include + +namespace simgear { + +class BucketBox; + +class ReaderWriterPGT : public osgDB::ReaderWriter { +public: + ReaderWriterPGT(); + virtual ~ReaderWriterPGT(); + + virtual const char* className() const; + + virtual osgDB::ReaderWriter::ReadResult readObject(const std::string& fileName, const osgDB::Options* options) const; + virtual osgDB::ReaderWriter::ReadResult readNode(const std::string& fileName, const osgDB::Options* options) const; + +protected: + struct LocalOptions; + + osg::ref_ptr createTree(const BucketBox& bucketBox, const LocalOptions& options, bool topLevel) const; + osg::ref_ptr createPagedLOD(const BucketBox& bucketBox, const LocalOptions& options) const; + osg::ref_ptr createTileMesh(const BucketBox& bucketBox, const LocalOptions& options) const; + osg::ref_ptr getLowLODStateSet(const LocalOptions& options) const; + +private: + struct CullCallback; +}; + +} // namespace simgear + +#endif diff --git a/simgear/scene/dem/SGDem.cxx b/simgear/scene/dem/SGDem.cxx new file mode 100644 index 00000000..d5c8c773 --- /dev/null +++ b/simgear/scene/dem/SGDem.cxx @@ -0,0 +1,332 @@ +// SGDem.cxx -- read, write DEM heiarchy +// +// Written by Peter Sadrozinski, started August 2016. +// +// Copyright (C) 2001 - 2003 Curtis L. Olson - http://www.flightgear.org/~curt +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ +#include +#include + +#include + +#include // for CPLMalloc() +#include "ogr_spatialref.h" + +#include +#include +#include + +#include +#include + +using namespace simgear; + +#define DEM_DEBUG (0) + +// where to add these... +int SGDem::floorWithEpsilon(double x) +{ + return static_cast(floor(x + SG_EPSILON)); +} + +unsigned SGDem::normalizeLongitude(unsigned offset) +{ + return offset - (360*8)*(offset/(360*8)); +} + +unsigned SGDem::longitudeDegToOffset(double lon) +{ + unsigned offset = (unsigned)( floorWithEpsilon( 8.0 * (lon + 180.0) ) ); + return normalizeLongitude(offset); +} + +double SGDem::offsetToLongitudeDeg(unsigned offset) +{ + return offset*0.125 - 180; +} + +unsigned SGDem::latitudeDegToOffset(double lat) +{ + if (lat < -90) + return 0; + + unsigned offset = (unsigned)( floorWithEpsilon( 8.0 * (lat + 90.0) ) ); + if (8*180 < offset) + return 8*180; + + return offset; +} + +double SGDem::offsetToLatitudeDeg(unsigned offset) +{ + return offset*0.125 - 90; +} + +int SGDem::addRoot( const SGPath& root ) +{ + GDALAllRegister(); + + SGDemRoot demRoot( root ); + + // collect subdir for each dem level - format is level_X + if ( root.isDir() ) { + Dir d(root); + if (d.exists() ) { + PathList levelPaths = d.children(Dir::TYPE_DIR); + SG_LOG( SG_TERRAIN, SG_INFO, levelPaths.size() << " Directories in " << d.path() ); + + for(const SGPath& p : levelPaths) { + std::string prefix; + int level; + + std::istringstream iss( p.file() ); + getline(iss, prefix, '_'); + if ( ((iss.rdstate() & std::ifstream::failbit ) == 0 ) && + (prefix == "level" ) ) { + iss >> level; + + // read the deminfo.txt file + SGPath infoFile = p / "deminfo.txt"; + SGPath extentsFile = p / "demextents.bin"; + + if ( infoFile.exists() ) { + std::fstream demInfo(infoFile.c_str(), std::ios_base::in); + int w, h, x, y, o; + std::string ext; + unsigned long extents[45][180]; + + demInfo >> w >> h >> x >> y >> o >> ext; + + // search level for extents + memset( (unsigned char*)extents, 0, sizeof(extents) ); + if ( extentsFile.exists() ) { + std::fstream demExtents(extentsFile.c_str(), std::ios_base::in | std::ios_base::binary ); + + for ( unsigned int i=0; i<45; i++ ) { + for ( unsigned int j=0; j<180; j++ ) { + demExtents >> extents[i][j]; + } + } + } + + SG_LOG( SG_TERRAIN, SG_INFO, " found DEM level " << level << " in directory " << p << " width : " << w << " height : " << h << " xres : " << x << " yres : " << y << " overlap : " << o << " extension: " << ext ); + demRoot.addLevel( SGDemLevel( level, p, w, h, x, y, o, &extents[0][0], ext, true ) ); + + } else { + SG_LOG( SG_TERRAIN, SG_INFO, " found DEM level " << level << " in directory " << p << " without info file " ); + } + + } else { + SG_LOG( SG_TERRAIN, SG_INFO, " invalid dem level dir " << p << " got prefix " << prefix ); + } + } + } + } + + if ( demRoot.numLevels() ) { + demRoots.push_back( demRoot ); + } + + return demRoot.numLevels(); +} + +int SGDem::createRoot( const SGPath& root ) +{ + // collect subdir for each dem level - format is level_X + // create the directory + SGPath newDir = root / "dummy"; + newDir.create_dir(); + + demRoots.push_back( SGDemRoot(root) ); + + return 0; +} + +unsigned SGDem::roundDown( unsigned offset, unsigned roundTo ) +{ + if ( roundTo == 0 ) { + return offset; + } else { + return (offset / roundTo) * roundTo; + } +} + +unsigned SGDem::roundUp( unsigned offset, unsigned roundTo ) +{ + if ( roundTo == 0 ) { + return offset; + } else { + return ((offset+roundTo-1) / roundTo) * roundTo; + } +} + +SGDemSession SGDem::openSession( unsigned wo, unsigned so, unsigned eo, unsigned no, int level, bool cache ) +{ + SGDemSession s; + + // Create the session + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::OpenSession - level " << level ); + + if ( level >= 0 ) { + SGDemRoot* demRoot = findDem( wo, so, eo, no, level ); + + if ( demRoot ) { + // traverse the demRoots, to see if any have this level + unsigned w = demRoot->getWidth( level ); + unsigned h = demRoot->getHeight( level ); + unsigned x = demRoot->getResX( level ); + unsigned y = demRoot->getResY( level ); + unsigned o = demRoot->getOverlap( level ); + + s = SGDemSession(wo, so, eo, no, level, w, h, demRoot); + + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::OpenSession - from offsets " << wo << ", " << so << " to offsets " << eo << ", " << no ); + + // calc min offsets based on tile widths and height for this level + unsigned min_lon = roundDown( wo, w ); + unsigned min_lat = roundDown( so, h ); + unsigned max_lon = roundUp( eo, w ); + unsigned max_lat = roundUp( no, h ); + + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::OpenSession - from pre level rounding offsets " << min_lon << ", " << min_lat << " to offsets " << max_lon << ", " << max_lat ); + + for (unsigned lon = min_lon; lon < max_lon; lon += w) { + for (unsigned lat = min_lat; lat < max_lat; lat += h) { + s.addTile( demRoot->getOrCreateTile( lon, lat, w, h, x, y, o, level, cache ) ); + } + } + } else { + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::OpenSession - could not find DEM for " << wo << ", " << so << " - " << eo << ", " << no << " level " << level ); + } + } else { + SG_LOG( SG_TERRAIN, SG_ALERT, "SGDem::OpenSession - invalid level " << level ); + } + + return s; +} + +SGDemSession SGDem::openSession( const SGGeod& min, const SGGeod& max, const SGPath& input ) +{ +#define FP_ROUNDOFF_OUTSIDE (0.1) + + // create a new demRoot for creation + SGDemRoot demRoot(input); + + // open all tiles between min and max + int min_lon = (int)(floor(min.getLongitudeDeg()-FP_ROUNDOFF_OUTSIDE)); + int min_lat = (int)(floor(min.getLatitudeDeg()-FP_ROUNDOFF_OUTSIDE)); + int max_lon = (int)(ceil(max.getLongitudeDeg()+FP_ROUNDOFF_OUTSIDE)); + int max_lat = (int)(ceil(max.getLatitudeDeg()+FP_ROUNDOFF_OUTSIDE)); + + // Create the session + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::OpenSession - create sesion obj - req from " << + min.getLongitudeDeg() << ", " << min.getLatitudeDeg() << " to " << + max.getLongitudeDeg() << ", " << max.getLatitudeDeg() << " - getting " << + min_lon << ", " << min_lat << " to " << max_lon << ", " << max_lat ); + + SGDemSession s(min_lon, min_lat, max_lon, max_lat, &demRoot); + + // todo - read a tile from the input dir + int w = 1; + int h = 1; + int x = 1201; + int y = 1201; + int o = 32; + + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::OpenSession - Traverse tiles"); + for (int lon = min_lon; lon < max_lon; lon += w) { + for (int lat = min_lat; lat < max_lat; lat += h) { + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::OpenSession - Create tile " << + lon << ", " << lat << " from dir " << input ); + + unsigned wo = (lon+180)*8; + unsigned so = (lat+ 90)*8; + + SGDemTile* pTile = new SGDemTile( input, wo, so, w, h, x, y, o, false ); + s.addTile( pTile ); + } + } + + return s; +} + +SGDemRoot* SGDem::findDem( unsigned wo, unsigned so, unsigned eo, unsigned no, int lvl ) +{ + SGDemRoot* dr = NULL; + + for ( unsigned int i=0; i= 0 ) { + // be careful with coordinates that lie on session boundaries - + // on min, ok. + // on max - make sure we select the tile in the session... + int lvlWidth = levels[lvlIndex].getWidth(); + int lvlHeight = levels[lvlIndex].getHeight(); + + int intLon, intLat; + + // we need to find the correct tile. + // shift by 180, 90 to use 0 based ints + intLon = (int)(round(loc.getLongitudeDeg())) + 180; + intLon = (intLon / lvlWidth) * lvlWidth; + intLon = intLon - 180; + if ( intLon == s.getMaxLon() ) { + intLon -= levels[lvlIndex].getWidth(); + } + + intLat = (int)(round(loc.getLatitudeDeg())) + 90; + intLat = (intLat / lvlHeight) * lvlHeight; + intLat = intLat - 90; + if ( intLat == s.getMaxLat() ) { + intLat -= levels[lvlIndex].getHeight(); + } + + unsigned long key = (intLon + 180) << 16 | (intLat + 90); + + SGDemCache::const_iterator it = caches[lvlIndex].find(key); + if ( it != caches[lvlIndex].end() ) { + SGDemTileRef tile = it->second; + alt = tile->getAlt( loc ); + } else { + // fprintf( stderr, " Could NOT find tile %d,%d in cache %d key is %08lx\n", intLon, intLat, lvlIndex, key ); + alt = 0; + } + } + + return alt; +} +#endif diff --git a/simgear/scene/dem/SGDem.hxx b/simgear/scene/dem/SGDem.hxx new file mode 100644 index 00000000..4fc6b2b7 --- /dev/null +++ b/simgear/scene/dem/SGDem.hxx @@ -0,0 +1,82 @@ +// SGDem.hxx -- read, write DEM heiarchy +// +// Written by Peter Sadrozinski, started August 2016. +// +// Copyright (C) 2001 - 2003 Curtis L. Olson - http://www.flightgear.org/~curt +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + + +#ifndef __SG_DEM_HXX__ +#define __SG_DEM_HXX__ + +#include + +#include +#include +#include + +class SGDem : public SGReferenced +{ +public: + SGDem() {}; + ~SGDem() {}; + + int addRoot( const SGPath& root ); + int createRoot( const SGPath& root ); + SGDemRoot* getRoot( unsigned int i ) { + if ( i < demRoots.size() ) { + return &demRoots[i]; + } else { + return NULL; + } + } + + unsigned int getNumRoots( void ) { + return demRoots.size(); + } + + // todo : move to session + // unsigned short getAlt( const SGDemSession& s, const SGGeod& loc ) const; + + // find a Dem to satisfy a session - must have the level, and extents + SGDemRoot* findDem( unsigned wo, unsigned so, unsigned eo, unsigned no, int lvl ); + + // open a session from a dem level - tiles will be read and reference counted until closed + //SGDemSession openSession( const SGGeod& min, const SGGeod& max, int level, bool cache ); + SGDemSession openSession( unsigned wo, unsigned so, unsigned eo, unsigned no, int level, bool cache ); + + // open a session from an bare directory + SGDemSession openSession( const SGGeod& min, const SGGeod& max, const SGPath& input ); + + // static helpers + static int floorWithEpsilon( double x ); + static unsigned normalizeLongitude( unsigned offset ); + static unsigned longitudeDegToOffset( double lon ); + static double offsetToLongitudeDeg( unsigned offset ); + static unsigned latitudeDegToOffset( double lat ); + static double offsetToLatitudeDeg( unsigned offset ); + static unsigned roundDown( unsigned offset, unsigned roundTo ); + static unsigned roundUp( unsigned offset, unsigned roundTo ); + +private: + std::vector demRoots; +}; + +typedef SGSharedPtr SGDemPtr; + +#endif /* __SG_DEM_HXX__ */ diff --git a/simgear/scene/dem/SGDemLevel.cxx b/simgear/scene/dem/SGDemLevel.cxx new file mode 100644 index 00000000..92abe7bf --- /dev/null +++ b/simgear/scene/dem/SGDemLevel.cxx @@ -0,0 +1,72 @@ +#include +#include + +#include +#include + +#include + +bool SGDemLevel::isValid( unsigned int wo, unsigned int so, unsigned int eo, unsigned int no ) const +{ + bool valid = true; + + return valid; +} + +void SGDemLevel::addExtent( unsigned int wo, unsigned int so, unsigned int eo, unsigned int no ) +{ + // add bitmap to level extent + unsigned int minLon = (wo+0)/8; + unsigned int maxLon = (eo+0)/8; + + unsigned int minLat = (so+0)/8; + unsigned int maxLat = (no+0)/8; + + fprintf(stderr, "addExtent %u, %u - %u, %u - minLon %u, minLat %u, maxLon %u, maxLat %u\n", wo, so, eo, no, minLon, minLat, maxLon, maxLat ); + for ( unsigned int lat=minLat; lat<=maxLat; lat++ ) { + // convert lon range to mask range + unsigned char minMask = minLon/8; + unsigned char maxMask = maxLon/8; + + fprintf(stderr, "addExtent minMask %u, maxMask %d\n", minMask, maxMask ); + for ( unsigned int msk=minMask; msk<=maxMask; msk++ ) { + if ( msk == minMask ) { + // find beginning of mask + unsigned char bitPos = minLon % 8; + unsigned char bitMask = 0; + for ( unsigned char i=0; i<=bitPos; i++ ) { + bitMask = (bitMask >> 1) | 0x80; + } + fprintf(stderr, "addExtent minMask - lat %u, lonmsk %u, bits %u\n", lat, msk, bitMask ); + extent[lat][msk] |= bitMask; + } else if ( msk == maxMask ) { + // find end of mask + unsigned char bitPos = maxLon % 8; + unsigned char bitMask = 0; + for ( unsigned int i=0; i<=bitPos; i++ ) { + bitMask = (bitMask >> 1) | 0x80; + } + fprintf(stderr, "addExtent maxMask - lat %u, lonmsk %u, bits %x\n", lat, msk, bitMask ); + extent[lat][msk] |= bitMask; + } else { + fprintf(stderr, "addExtent middlemask - lat %u, lonmsk %u, bits %x\n", lat, msk, 0xFF ); + extent[lat][msk] = 0xFF; + } + } + } +} + +void SGDemLevel::close( void ) +{ + // save extent to file + fprintf( stderr, "closing level\n" ); + SGPath extentFile = path / "demextent.txt"; + std::fstream demExtent(extentFile.c_str(), std::ios_base::out | std::ios_base::binary ); + + for ( unsigned char lat = 0; lat < 180; lat++ ) { + for ( unsigned char lon = 0; lon < 360/8; lon++ ) { + demExtent << extent[lat][lon]; + } + } + fprintf( stderr, "closing level complete\n" ); +} diff --git a/simgear/scene/dem/SGDemLevel.hxx b/simgear/scene/dem/SGDemLevel.hxx new file mode 100644 index 00000000..2cd426cd --- /dev/null +++ b/simgear/scene/dem/SGDemLevel.hxx @@ -0,0 +1,90 @@ +#ifndef __SG_DEM_LEVEL_HXX__ +#define __SG_DEM_LEVEL_HXX__ + +#include +#include + +#include + +class SGDemLevel +{ +public: + SGDemLevel( int l, const SGPath& p, int w, int h, int x, int y, int o, const unsigned long* e, std::string ext, bool r ) + { + level = l; + path = p; + width = w; + height = h; + xres = x; + yres = y; + overlap = o; + extension = ext; + ready = r; + + std::fprintf( stderr, "copying extent fle size %lu\n", sizeof(extent) ); + std::memcpy( (void *)&extent[0][0], (const void *)e, sizeof(extent) ); + }; + + SGDemLevel( int l, const SGPath& p, int w, int h, int x, int y, int o, std::string ext, bool r ) + { + level = l; + path = p; + width = w; + height = h; + xres = x; + yres = y; + overlap = o; + extension = ext; + ready = r; + + std::fprintf( stderr, "setting extent fle size %lu\n", sizeof(extent) ); + memset( (unsigned char*)extent, 0, sizeof(extent) ); + }; + + int getLevel( void ) const { + return level; + } + const SGPath& getLevelDir( void ) const { + return path; + } + bool isReady( void ) const { + return ready; + } + bool isValid( unsigned int wo, unsigned int so, unsigned int eo, unsigned int no ) const; + + void addExtent( unsigned int wo, unsigned int so, unsigned int eo, unsigned int no ); + void close( void ); + + int getWidth( void ) const { + return width; + } + int getHeight( void ) const { + return height; + } + int getResX( void ) const { + return xres; + } + int getResY( void ) const { + return yres; + } + int getOverlap( void ) const { + return overlap; + } + +private: + SGPath path; + std::string extension; + bool ready; + + int level; + int width; + int height; + int xres; + int yres; + int overlap; + + // bitmap of tiles (1x1 degree) + unsigned char extent[180][360/8]; +}; + +#endif /* __SG_DEM_LEVEL_HXX__ */ diff --git a/simgear/scene/dem/SGDemRoot.cxx b/simgear/scene/dem/SGDemRoot.cxx new file mode 100644 index 00000000..ff9b21f9 --- /dev/null +++ b/simgear/scene/dem/SGDemRoot.cxx @@ -0,0 +1,118 @@ +#include +#include +#include + +#include +#include +#include + +bool SGDemRoot::isValid( int lvl, unsigned wo, unsigned so, unsigned eo, unsigned no ) const +{ + bool valid = false; + + // check if we have requested level + if ( (unsigned int)lvl < levels.size() ) { + // check if we fit in this levels extents + valid = levels[lvl].isValid( wo, so, eo, no ); + } else { + SG_LOG( SG_TERRAIN, SG_INFO, "SGDemRoot::isValid - lvl " << lvl << " is greater than #levels in root " << levels.size() ); + } + + return valid; +} + +int SGDemRoot::createLevel( int w, int h, int x, int y, int overlap, const std::string& ext ) +{ + int lvlidx = -1; + + std::stringstream ss; + ss << "level_" << std::setw(2) << std::setfill('0') << levels.size()+1; + + // see if dir already exists + SGPath lvlPath = demRoot / ss.str(); + if ( lvlPath.isDir() ) { + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::createLevel: found existing level directory at " << lvlPath.c_str() ); + } else { + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::createLevel: creating level directory at " << lvlPath.c_str() ); + + // create the dir - needs the filename to do this !?!? + SGPath infoFile = lvlPath / "deminfo.txt"; + int res = infoFile.create_dir(); + if ( res == 0 ) { + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::createLevel: successfully created directory at " << infoFile.realpath().c_str() << " error " << res ); + + lvlidx = levels.size(); + + std::fstream demInfo(infoFile.c_str(), std::ios_base::out); + demInfo << w << " " << h << " " << x << " " << y << " " << overlap << " " << ext << std::endl; + + levels.push_back( SGDemLevel( lvlidx, lvlPath, w, h, x, y, overlap, ext, false ) ); + caches.push_back( SGDemCache() ); + } else { + SG_LOG( SG_TERRAIN, SG_INFO, "SGDem::createLevel: can't create level directory at " << lvlPath.c_str() << " error " << res ); + } + } + + return lvlidx; +} + +void SGDemRoot::closeLevel( int lvl ) +{ + if ( (unsigned int)lvl < levels.size() ) { + levels[lvl].close(); + } +} + +SGDemTileRef SGDemRoot::createTile( int lvl, int lon, int lat, int overlap, SGDemSession& s ) +{ + SGDemTileRef rTile = NULL; + bool bWritten = false; + + if ( lvl >= 0 ) { + int w = levels[lvl].getWidth(); + int h = levels[lvl].getHeight(); + int x = levels[lvl].getResX(); + int y = levels[lvl].getResY(); + + unsigned wo = (unsigned)(lon+180)*8; + unsigned so = (unsigned)(lat+ 90)*8; + + rTile = new SGDemTile( levels[lvl].getLevelDir(), + wo, so, + w, h, x, y, overlap, s, bWritten ); + + if ( bWritten ) { + fprintf( stderr, "CreateTile - add extent (lat%d+h%d+90)*8-1 is %u\n", lat, h, (lat+h+90)*8-1 ); + levels[lvl].addExtent( (lon+180)*8, (lat+90)*8, (lon+w+180)*8-1, (lat+h+90)*8-1 ); + } + } + + return rTile; +} + +void SGDemRoot::flushCaches( int lvl ) +{ + // traverse the map - delete any tiles with ref count == 1 ( us ) + //fprintf( stderr, "flush caches in lvl %d - num tiles is %lu\n", lvl, caches[lvl].size() ); + SGDemCache::iterator it = caches[lvl].begin(); + while ( it != caches[lvl].end() ) { + if ( it->second.getNumRefs() == 1 ) { + caches[lvl].erase(it++); + } else { + //fprintf( stderr, "can't flush tile - numrefs is %u\n", it->second.getNumRefs() ); + ++it; + } + } +} + +SGDemTileRef SGDemRoot::getTile( int lvlIndex, unsigned long key ) +{ + SGDemTileRef tile = NULL; + + SGDemCache::const_iterator it = caches[lvlIndex].find(key); + if ( it != caches[lvlIndex].end() ) { + tile = it->second; + } + + return tile; +} diff --git a/simgear/scene/dem/SGDemRoot.hxx b/simgear/scene/dem/SGDemRoot.hxx new file mode 100644 index 00000000..1ecc1593 --- /dev/null +++ b/simgear/scene/dem/SGDemRoot.hxx @@ -0,0 +1,102 @@ +#ifndef __SG_DEM_ROOT_HXX__ +#define __SG_DEM_ROOT_HXX__ + +#include +#include + +class SGDemRoot +{ +public: + SGDemRoot( SGPath p ) { + demRoot = p; + } + + void addLevel( const SGDemLevel& level ) { + levels.push_back( level ); + caches.push_back( SGDemCache() ); + } + + // create a new level + int createLevel( int w, int h, int x, int y, int overlap, const std::string& ext ); + SGDemTileRef createTile( int lvl, int lon, int lat, int overlap, SGDemSession& s ); + void closeLevel( int lvl ); + + unsigned int numLevels( void ) const { + return levels.size(); + } + + SGDemTileRef getTile( int lvlIndex, unsigned long key ); + + void flushCaches( int lvl ); + + bool isValid( int lvl, unsigned wo, unsigned so, unsigned eo, unsigned no ) const; + + unsigned int getWidth( unsigned int level ) { + if ( level < levels.size() ) { + return levels[level].getWidth() * 8; + } else { + return 0; + } + } + + unsigned int getHeight( unsigned int level ) { + if ( level < levels.size() ) { + return levels[level].getHeight() * 8; + } else { + return 0; + } + } + + unsigned int getResX( unsigned int level ) { + if ( level < levels.size() ) { + return levels[level].getResX(); + } else { + return 0; + } + } + + unsigned int getResY( unsigned int level ) { + if ( level < levels.size() ) { + return levels[level].getResY(); + } else { + return 0; + } + } + + unsigned int getOverlap( unsigned int level ) { + if ( level < levels.size() ) { + return levels[level].getOverlap(); + } else { + return 0; + } + } + + SGDemTile* getOrCreateTile( unsigned wo, unsigned so, + unsigned w, unsigned h, unsigned x, unsigned y, + unsigned o, int level, bool cache ) + { + unsigned long key = wo << 16 | so; + + // is this tile already loaded? + SGDemCache::const_iterator it = caches[level].find(key); + if ( it != caches[level].end() ) { + // yes - add the reference to this session + //printf( "********************** adding ref to tile at %d,%d with key %lx\n", lon, lat, key ); + return it->second; + } else { + // no load the tile, now + SGDemTile* pTile = new SGDemTile( levels[level].getLevelDir(), + wo, so, w, h, + x, y, o, cache ); + caches[level][key] = pTile; + return pTile; + } + } + +private: + SGPath demRoot; + std::vector levels; + std::vector caches; +}; + +#endif /* __SG_DEM_ROOT_HXX__ */ diff --git a/simgear/scene/dem/SGDemSession.cxx b/simgear/scene/dem/SGDemSession.cxx new file mode 100644 index 00000000..fb95caec --- /dev/null +++ b/simgear/scene/dem/SGDemSession.cxx @@ -0,0 +1,79 @@ +#include +#include + +SGDemSession::SGDemSession( int mnLon, int mnLat, int mxLon, int mxLat, int idx, int lvlW, int lvlH, SGDemRoot* root ) +{ + setOffsets( SGDem::longitudeDegToOffset((double)mnLon), + SGDem::latitudeDegToOffset((double)mnLat), + SGDem::longitudeDegToOffset((double)mxLon), + SGDem::latitudeDegToOffset((double)mxLat) ); + + pDemRoot = root; + lvlIndex = idx; + lvlWidth = lvlW; + lvlHeight = lvlH; +} + +SGDemSession::SGDemSession( int mnLon, int mnLat, int mxLon, int mxLat, SGDemRoot* root ) { + setOffsets( SGDem::longitudeDegToOffset((double)mnLon), + SGDem::latitudeDegToOffset((double)mnLat), + SGDem::longitudeDegToOffset((double)mxLon), + SGDem::latitudeDegToOffset((double)mxLat) ); + + pDemRoot = root; + lvlIndex = -1; // no level - session is raw input dir +} + +void SGDemSession::close( void ) +{ + if ( tileRefs.size() ) { + tileRefs.clear(); + if ( lvlIndex >= 0 ) { + pDemRoot->flushCaches( lvlIndex ); + } + } +} + +void SGDemSession::getGeods( unsigned wo, unsigned so, unsigned eo, unsigned no, int resx, int resy, int incx, int incy, ::std::vector& geods, bool Debug1, bool Debug2 ) +{ + // todo - store this info in deminfo + unsigned span; // smallest tile width/height in level ( in offsets ) + switch( lvlIndex ) { + case 0: span = 1; break; // 1/8 deg + case 1: span = 16; break; // 2 degrees + case 2: span = 480; break; // 60 degrees + default: + fprintf( stderr, "invalid lvlIndex %d\n", lvlIndex ); + exit(0); + } + + if ( lvlIndex >= 0 ) { + unsigned tileLon, tileLat; + unsigned meshLon, meshLat; + int subx, suby; + + meshLon = wo; + tileLon = SGDem::roundDown( meshLon, lvlWidth); + subx = (meshLon - tileLon) / span; + // fprintf(stderr, "getGeods: lon is %lf : meshLon is %u, tileLon is %u, subx is %d\n", SGDem::offsetToLongitudeDeg(wo), meshLon, tileLon, subx ); + + meshLat = so; + tileLat = SGDem::roundDown( meshLat, lvlHeight ); + suby = (meshLat - tileLat) / span; + // fprintf(stderr, "getGeods: lat is %lf : meshLat is %u, tileLat is %u, suby is %d\n", SGDem::offsetToLatitudeDeg(so), meshLat, tileLat, suby ); + + // get the tle from the tile cache + unsigned long key = tileLon << 16 | tileLat; + SGDemTileRef tile = pDemRoot->getTile( lvlIndex, key ); + if ( tile ) { + tile->getGeods(wo, so, eo, no, resx, resy, subx, suby, incx, incy, geods, Debug1, Debug2); + } else { + fprintf(stderr, " *** ERROR: tile %d,%d not in session @ (%lf,%lf) - (%lf,%lf)\n", + tileLon, tileLat, + SGDem::offsetToLongitudeDeg( west_off ), + SGDem::offsetToLatitudeDeg( south_off ), + SGDem::offsetToLongitudeDeg( east_off ), + SGDem::offsetToLatitudeDeg( north_off ) ); + } + } +} diff --git a/simgear/scene/dem/SGDemSession.hxx b/simgear/scene/dem/SGDemSession.hxx new file mode 100644 index 00000000..f0813556 --- /dev/null +++ b/simgear/scene/dem/SGDemSession.hxx @@ -0,0 +1,78 @@ +#ifndef __SG_DEM_SESSION_HXX__ +#define __SG_DEM_SESSION_HXX__ + +#include + +class SGDemSession +{ +public: + SGDemSession() { + lvlIndex = -1; + } + + SGDemSession( int mnLon, int mnLat, int mxLon, int mxLat, int idx, int lvlW, int lvlH, SGDemRoot* root ); + SGDemSession( int mnLon, int mnLat, int mxLon, int mxLat, SGDemRoot* root ); + + SGDemSession( unsigned wo, unsigned so, unsigned eo, unsigned no, int idx, unsigned lvlW, unsigned lvlH, SGDemRoot* root ) { + setOffsets( wo, so, eo, no ); + + pDemRoot = root; + lvlIndex = idx; + lvlWidth = lvlW; + lvlHeight = lvlH; + } + + SGDemSession( unsigned wo, unsigned so, unsigned eo, unsigned no, SGDemRoot* root ) { + setOffsets( wo, so, eo, no ); + + pDemRoot = root; + lvlIndex = -1; // no level - session is raw input dir + } + + ~SGDemSession() { + close(); + } + + void addTile(SGDemTileRef pTile) { + tileRefs.push_back( pTile ); + } + + const std::vector& getTiles( void ) const { + return tileRefs; + } + + unsigned int size( void ) const { + return tileRefs.size(); + } + + void getGeods( unsigned wp, unsigned so, unsigned eo, unsigned no, + int resx, int resy, int incx, int incy, + ::std::vector& geods, + bool Debug1, bool Debug2 + ); + + void close( void ); + + int getLvlIndex( void ) const { + return lvlIndex; + }; + +private: + void setOffsets( unsigned wo, unsigned so, unsigned eo, unsigned no ) { + west_off = wo; + south_off = so; + east_off = eo; + north_off = no; + } + + unsigned west_off, south_off; + unsigned east_off, north_off; + int maxLon, maxLat; + SGDemRoot* pDemRoot; + int lvlIndex; + unsigned lvlWidth, lvlHeight; + + std::vector tileRefs; +}; + +#endif /* __SG_DEM_SESSION_HXX__ */ diff --git a/simgear/scene/dem/SGDemTile.cxx b/simgear/scene/dem/SGDemTile.cxx new file mode 100644 index 00000000..ec0bd0b9 --- /dev/null +++ b/simgear/scene/dem/SGDemTile.cxx @@ -0,0 +1,361 @@ +#include +#include + +#include +#include + +#include +#include +#include + +SGDemTile::SGDemTile( const SGPath& levelDir, unsigned wo, unsigned so, int w, int h, int x, int y, int o, bool cache) +{ + // add the tile name to the level path + ref_lon = (int)wo/8 - 180; + ref_lat = (int)so/8 - 90; + path = levelDir / getTileName( ref_lon, ref_lat ); + width = w; + height = h; + resx = x; + resy = y; + overlap = o; + + pixResX = ((double)width/(double)8.0)/(double)(resx-1); + pixResY = ((double)height/(double)8.0)/(double)(resy-1); + + if ( cache ) { + raster = cacheTile( path ); + } else { + raster = NULL; + } +} + +void SGDemTile::dbgDumpDataset( GDALDataset* poDataset ) const +{ + double adfGeoTransform[6]; + + SG_LOG( SG_TERRAIN, SG_INFO, "Driver: " << poDataset->GetDriver()->GetDescription() << "/" << poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) ); + SG_LOG( SG_TERRAIN, SG_INFO, "Size is " << poDataset->GetRasterXSize() << " x " << poDataset->GetRasterYSize() << " x " << poDataset->GetRasterCount() ); + + if( poDataset->GetProjectionRef() != NULL ) { + SG_LOG( SG_TERRAIN, SG_INFO, "Projection is " << poDataset->GetProjectionRef() ); + } + + if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None ) { + SG_LOG( SG_TERRAIN, SG_INFO, "Origin = (" << adfGeoTransform[0] << ", " << adfGeoTransform[3] << ")" ); + SG_LOG( SG_TERRAIN, SG_INFO, "Pixel Size = (" << adfGeoTransform[1] << ", " << adfGeoTransform[5] << ")" ); + } +} + +void SGDemTile::dbgDumpBand( GDALRasterBand* poBand ) const +{ + int nBlockXSize, nBlockYSize; + int bGotMin, bGotMax; + double adfMinMax[2]; + + poBand->GetBlockSize( &nBlockXSize, &nBlockYSize ); + SG_LOG( SG_TERRAIN, SG_INFO, "Block=" << nBlockXSize << " x " << nBlockYSize << " Type=" << GDALGetDataTypeName(poBand->GetRasterDataType()) << "ColorInterp=" << GDALGetColorInterpretationName( poBand->GetColorInterpretation()) ); + + adfMinMax[0] = poBand->GetMinimum( &bGotMin ); + adfMinMax[1] = poBand->GetMaximum( &bGotMax ); + if( ! (bGotMin && bGotMax) ) { + GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax); + } + SG_LOG( SG_TERRAIN, SG_INFO, "Min=" << adfMinMax[0] << ", Max=" << adfMinMax[1] ); + if( poBand->GetOverviewCount() > 0 ) { + SG_LOG( SG_TERRAIN, SG_INFO, "Band has " << poBand->GetOverviewCount() << " overviews." ); + } + if( poBand->GetColorTable() != NULL ) { + SG_LOG( SG_TERRAIN, SG_INFO, "Band has a color table with " << poBand->GetColorTable()->GetColorEntryCount() << " entries." ); + } +} + +unsigned short* SGDemTile::cacheTile( const SGPath& path ) +{ + unsigned short* r = NULL; + + GDALDataset* poDataset = NULL; + GDALRasterBand* poBand = NULL; + + // check if file exists to supress GDAL errors... + if ( path.exists() ) { + poDataset = (GDALDataset *)GDALOpen( path.c_str(), GA_ReadOnly ); + } + + if ( poDataset ) { + +#if DEM_DEBUG + dbgDumpDataset( poDataset ); +#endif + + // read the bands into array + unsigned int numRasters = poDataset->GetRasterCount(); + for ( unsigned int rb=1; rb<=numRasters; rb++ ) { + poBand = poDataset->GetRasterBand(rb); + +#if DEM_DEBUG + dbgDumpBand( poBand ); +#endif + + // if this is the raster we're looking for + if ( rb == 1 ) { + int nXSize = poBand->GetXSize(); + int nYSize = poBand->GetYSize(); + // SG_LOG( SG_TERRAIN, SG_INFO, "Reading raster " << rb << " (" << nXSize << "x" << nYSize << ")" ); + + r = (unsigned short *)CPLMalloc(sizeof(unsigned short)*nXSize*nYSize); + //fprintf( stderr, "reading raster size %d x %d - buffer is %p\n", nXSize, nYSize, r ); + CPLErr err = poBand->RasterIO( GF_Read, 0, 0, nXSize, nYSize, r, nXSize, nYSize, GDT_UInt16, 0, 0 ); + if ( err ) { + fprintf(stderr, "Error reading raster\n"); + } + } + } + + GDALClose( poDataset ); + } + + return r; +} + +// Create a new DEM tile from multiple source tiles ( at an expected lower resolution ) +// code based on gdalwarp, but with most options removed. +// example dgalwarp usage this function is based on: +// gdalwarp -r cubic -ts 1201 1201 -te -85.0 32.0 -83.0 34.0 -dstnodata 0 -co "COMPRESS=DEFLATE" temp/N32W085.hgt temp/merged_cubic.tiff + +SGDemTile::SGDemTile( const SGPath& levelDir, unsigned wo, unsigned so, int w, int h, int x, int y, int overlap, const SGDemSession& s, bool& bWritten ) +{ + fprintf( stderr, "Writing tile: lon offset is %u, lat offset is %u\n", wo, so ); + + // assume failure + bWritten = false; + + // add the tile name to the level path + ref_lon = (int)wo/8 - 180; + ref_lat = (int)so/8 - 90; + + path = levelDir / getTileName( ref_lon, ref_lat ); + width = w; + height = h; + resx = x; + resy = y; + + raster = NULL; + + // use gdal warp api to generate tile from session + char** papszSrcFiles = NULL; + char** papszWarpOptions = NULL; + char** papszTO = NULL; + + double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0; // -te + double overlapw = (double)w/(double)resx * overlap; + double overlaph = (double)h/(double)resy * overlap; + int nForcePixels = resx+(2*overlap), nForceLines = resy+(2*overlap); // -ts + + // -dstnodata + papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0"); + + // target extents ( +/- 1 pixel for normals, and skirts ) + dfMinX = (double)ref_lon - overlapw; + dfMinY = (double)ref_lat - overlaph; + dfMaxX = (double)ref_lon + w + overlapw; + dfMaxY = (double)ref_lat + h + overlaph; + + SG_LOG( SG_TERRAIN, SG_INFO, "overlapw: " << overlapw << " overlapw " << overlapw << " resx " << resx << " resy " << resy << " w " << w << " h " << h ); + SG_LOG( SG_TERRAIN, SG_INFO, " dfMinX: " << dfMinX << " dfMinY: " << dfMinY << " dfMaxX: " << dfMaxX << " dfMaxY: " << dfMaxY ); + + // generate list of source files in session + const std::vector& tiles = s.getTiles(); + SG_LOG( SG_TERRAIN, SG_INFO, " create tile from " << tiles.size() << " source tiles" ); + + for ( unsigned int i=0; igetPath().exists() ) { + papszSrcFiles = CSLAddString( papszSrcFiles, tiles[i]->getPath().c_str() ); + SG_LOG( SG_TERRAIN, SG_INFO, " Adding tile " << tiles[i]->getPath() ); + } else { + // SG_LOG( SG_TERRAIN, SG_INFO, " tile " << tiles[i]->getPath() << " doesn't exist" ); + } + } + + // create output + if ( papszSrcFiles ) { + GDALDatasetH hDstDS = createTile( papszSrcFiles, path.c_str(), nForceLines, nForcePixels, + dfMinX, dfMinY, dfMaxX, dfMaxY, + papszTO ); + + if( hDstDS != NULL ) { + /* -------------------------------------------------------------------- */ + /* Loop over all source files, processing each in turn. */ + /* -------------------------------------------------------------------- */ + int iSrc; + for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ ) + { + doWarp( iSrc, papszSrcFiles[iSrc], hDstDS, papszTO, papszWarpOptions ); + } + + // manually set min/max to all tiles are consistent when viewing in qgis/grass + // we need to get the raster band + GDALDataset* poDataset = (GDALDataset *)hDstDS; + GDALRasterBand* poBand = NULL; + + // read the bands into array + unsigned int numRasters = poDataset->GetRasterCount(); + for ( unsigned int rb=1; rb<=numRasters; rb++ ) { + poBand = poDataset->GetRasterBand(rb); + + // if this is the raster we're looking for + if ( rb == 1 ) { + double pdfMin, pdfMax, pdfMean, pdfStddev; + + poBand->ComputeStatistics(false, &pdfMin, &pdfMax, &pdfMean, &pdfStddev, NULL, NULL ); + + fprintf(stderr, "Got band min as %lf, max as %lf\n", pdfMin, pdfMax ); + + // force status to sea level / round of mnt everest + pdfMin = 0; pdfMax = 9000; + poBand->SetStatistics( pdfMin, pdfMax, pdfMean, pdfStddev); + + fprintf(stderr, "Setting band min to %lf, max to %lf\n", pdfMin, pdfMax ); + } + } + + /* -------------------------------------------------------------------- */ + /* Final Cleanup. */ + /* -------------------------------------------------------------------- */ + CPLErrorReset(); + GDALFlushCache( hDstDS ); + GDALClose( hDstDS ); + + CSLDestroy( papszSrcFiles ); + CSLDestroy( papszWarpOptions ); + CSLDestroy( papszTO ); + + GDALDumpOpenDatasets( stderr ); + + bWritten = true; + } + } +} + +SGDemTile::~SGDemTile() +{ + // free the raster + if ( raster ) { + CPLFree( raster ); + } +} + +unsigned short SGDemTile::getAlt( const SGGeod& loc ) const +{ + if ( raster ) { + bool debug = false; + + // get lon and lat reletive to sw corner + double offLon = loc.getLongitudeDeg() - (double)ref_lon; + double offLat = loc.getLatitudeDeg() - (double)ref_lat; + + // raster has a 1 pixel border on all sides + // take that into account + if ( fabs( offLon ) < 0.000001 ) { + debug = true; + } + + double fractLon = offLon / (double)width; + double fractLat = offLat / (double)height; + + int l = (int)( (double)(resy-1) - ((double)(resy-1)*fractLat) + 1 ); + int p = (int)( (double)(resx-1) * fractLon + 1 ); + + if ( debug ) { + printf( "SGDemTile::getAlt at %lf,%lf: offLon is %lf, offLat is %lf, width is %d, height is %d, fractLon is %lf, fractLat is %lf, resx is %d, resy is %d, line %d, pixel %d\n", + loc.getLongitudeDeg(), loc.getLatitudeDeg(), + offLon, offLat, width, height, fractLon, fractLat, + resx, resy, l, p ); + } + + return raster[l*(resx+2)+p]; + } else { + return 0; + } +} + +void SGDemTile::getGeods( unsigned wo, unsigned so, unsigned eo, unsigned no, int grid_width, int grid_height, unsigned subx, unsigned suby, int incw, int inch, ::std::vector& geods, bool Debug1, bool Debug2 ) +{ + // grid width and height include the skirt + // sw and ne do not + // we need to find the starting and ending l and p; + int startl;//, endl; + int startp;//, endp; + double startlat, startlon; + double endlat, endlon; + + startlon = SGDem::offsetToLongitudeDeg(wo) - incw*pixResX; + startlat = SGDem::offsetToLatitudeDeg(so) - inch*pixResY; + + endlon = SGDem::offsetToLongitudeDeg(eo) + incw*pixResX; + endlat = SGDem::offsetToLatitudeDeg(no) + inch*pixResY; + + // todo : how to calculate 15- from given data + if ( Debug1 ) { + printf("resx is %d resy is %d incw is %d incy is %d\n", resx, resy, incw, inch ); + } + if ( raster ) { + int di, dj; + int p, l; + double lonPos, latPos; + double maxlon = startlon; + double maxlat = startlat; + + startl = (resy-1) + overlap - ( suby * (153-3) ) + inch; + startp = 0 + overlap + ( subx * (153-3) ) - incw; + + for ( di = 0, p = startp, lonPos = startlon; di < grid_width; di++, p+=incw, lonPos += incw*pixResX ) { + maxlon = lonPos; + + for ( dj = 0, l = startl, latPos = startlat; dj < grid_height; dj++, l-=inch, latPos += inch*pixResY ) { + maxlat = latPos; + + SGGeod pos = SGGeod::fromDeg( SGMiscd::normalizePeriodic( -180.0, 180.0, lonPos ), + SGMiscd::normalizePeriodic( -180.0, 180.0, latPos ) ); + + if ( raster ) { + pos.setElevationM( raster[l*(resx+(2*overlap))+p] ); + } + + geods[di*grid_height+dj] = pos; + } + } + + if ( fabs( endlon - maxlon ) > 0.0001 ) { + printf(" tile overlap error %lf : lon is %lf, startlon is %lf, endlon is %lf, maxlon is %lf. grid_width is %d, incw is %d, pixResX is %lf, (grid_width-1)*incw*pixResX is %lf\n", + endlon-maxlon, SGDem::offsetToLongitudeDeg(wo), startlon, endlon, maxlon, grid_width, incw, pixResX, (grid_width-1)*incw*pixResX ); + } + + if ( fabs( endlat - maxlat ) > 0.0001 ) { + printf(" tile overlap error %lf : lat is %lf, startlat is %lf, endlat is %lf, maxlat is %lf. grid_height is %d, inch is %d, pixResY is %lf, (grid_height-1)*inch*pixResY is %lf\n", + endlat-maxlat, SGDem::offsetToLatitudeDeg(so), startlat, endlat, maxlat, grid_height, inch, pixResY, (grid_height-1)*inch*pixResY ); + } + } +} + +std::string SGDemTile::getTileName( int lon, int lat ) +{ + std::stringstream ss; + + if ( lat >= 0 ) { + ss << "N" << std::setw(2) << std::setfill('0') << lat; + } else { + ss << "S" << std::setw(2) << std::setfill('0') << -lat; + } + + if ( lon >= 0 ) { + ss << "E" << std::setw(3) << std::setfill('0') << lon; + } else { + ss << "W" << std::setw(3) << std::setfill('0') << -lon; + } + ss << ".hgt"; + + printf("created tile string %s from %d,%d\n", ss.str().c_str(), lon, lat ); + return ss.str(); +} diff --git a/simgear/scene/dem/SGDemTile.hxx b/simgear/scene/dem/SGDemTile.hxx new file mode 100644 index 00000000..fb2c5d22 --- /dev/null +++ b/simgear/scene/dem/SGDemTile.hxx @@ -0,0 +1,65 @@ +#ifndef __SG_DEM_TILE_HXX__ +#define __SG_DEM_TILE_HXX__ + +#include +#include + +#include + +#include +#include + +class SGDemSession; + +class SGDemTile : public SGReferenced +{ +public: + // TODO - simple constructor, so writing a tile to disk not done in constructor. + // then - reading / writing done with tile API, not constructor. + //SGDemTile( const SGPath& path, int lon, int lat, int w, int h, int x, int y, int overlap ); + + // constructor for reading a tile + SGDemTile( const SGPath& path, unsigned wo, unsigned so, int w, int h, int x, int y, int overlap, bool cache ); + // constructor for writing a tile + SGDemTile( const SGPath& path, unsigned wo, unsigned so, int w, int h, int x, int y, int overlap, const SGDemSession& s, bool& bWritten ); + + ~SGDemTile(); + + // read / write tile from / to disk + //int read( bool cache ); + //int write( const SGDemSession& s ); + + SGPath getPath( void ) const { return path; } + unsigned short getAlt(const SGGeod& loc) const; + void getGeods(unsigned wo, unsigned so, unsigned eo, unsigned no, int grid_width, int grid_height, unsigned subx, unsigned suby, int incw, int inch, ::std::vector& geods, bool Debug1, bool Debug2); + +private: + std::string getTileName( int lon, int lat ); + void dbgDumpDataset( GDALDataset* poDataset ) const; + void dbgDumpBand( GDALRasterBand* poBand ) const; + + unsigned short* cacheTile( const SGPath& path ); + + GDALDatasetH createTile( char **papszSrcFiles, + const char *pszFilename, + int nForceLines, int nForcePixels, + double dfMinX, double dfMinY, + double dfMaxX, double dfMaxY, + char **papszTO ); + + void doWarp( int iSrc, char* pszSrcFile, GDALDatasetH hDstDS, char **papszTO, char** papszWarpOptions ); + + SGPath path; + unsigned short* raster; + + int ref_lon, ref_lat; + int width, height; + int resx, resy; + int overlap; + double pixResX, pixResY; +}; + +typedef SGSharedPtr SGDemTileRef; +typedef std::map SGDemCache; + +#endif /* #define __SG_DEM_TILE_HXX__ */ diff --git a/simgear/scene/dem/SGDemTile_gdal.cxx b/simgear/scene/dem/SGDemTile_gdal.cxx new file mode 100644 index 00000000..bcc9f08d --- /dev/null +++ b/simgear/scene/dem/SGDemTile_gdal.cxx @@ -0,0 +1,628 @@ +// dem_gdal.cxx -- perform gdal warp - based on gdalwarp.cpp + +/****************************************************************************** + * $Id: gdalwarp.cpp 29153 2015-05-04 17:51:41Z rouault $ + * + * Project: High Performance Image Reprojector + * Purpose: Test program for high performance warper API. + * Author: Frank Warmerdam + * + ****************************************************************************** + * Copyright (c) 2002, i3 - information integration and imaging + * Fort Collin, CO + * Copyright (c) 2007-2013, Even Rouault + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ +#include +#include +#include + +#include + +#include // for CPLMalloc() +#include "ogr_spatialref.h" + +#include +#include +#include + +#include +#include +#include +#include + +using namespace simgear; + +int GDALExit( int nCode ) +{ + GDALDumpOpenDatasets( stderr ); + CPLDumpSharedList( NULL ); + + GDALDestroyDriverManager(); + OGRCleanupAll(); + + exit( nCode ); +} + +GDALDatasetH SGDemTile::createTile( char **papszSrcFiles, const char *pszFilename, + int nForceLines, int nForcePixels, + double dfMinX, double dfMinY, + double dfMaxX, double dfMaxY, + char **papszTO ) +{ + GDALDriverH hDriver = NULL; + GDALDatasetH hDstDS; + GDALColorTableH hCT = NULL; + double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0; + double dfXRes=0.0, dfYRes=0.0; + int nDstBandCount = 0; + std::vector apeColorInterpretations; + + /* -------------------------------------------------------------------- */ + /* Find the output driver. */ + /* -------------------------------------------------------------------- */ + hDriver = GDALGetDriverByName( "GTiff" ); + GDALDataType eDT = GDT_Unknown; + if( hDriver == NULL || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL ) + { + int iDr; + + printf( "Output driver 'GTiff' not recognised or does not support\n" ); + printf( "direct output file creation. The following format drivers are configured\n" + "and support direct output:\n" ); + + for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ ) + { + GDALDriverH hDriver = GDALGetDriver(iDr); + + if( GDALGetMetadataItem( hDriver, GDAL_DCAP_RASTER, NULL) != NULL && + GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL ) + { + printf( " %s: %s\n", + GDALGetDriverShortName( hDriver ), + GDALGetDriverLongName( hDriver ) ); + } + } + printf( "\n" ); + } + + char *pszThisTargetSRS = (char*)CSLFetchNameValue( papszTO, "DST_SRS" ); + if( pszThisTargetSRS != NULL ) { + pszThisTargetSRS = CPLStrdup( pszThisTargetSRS ); + } else { + printf("GDALWarpCreateOutput: pszThisTargetSRS is NULL\n"); + } + + /* -------------------------------------------------------------------- */ + /* Loop over all input files to collect extents. */ + /* -------------------------------------------------------------------- */ + for( int iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ ) + { + GDALDatasetH hSrcDS; + const char *pszThisSourceSRS = CSLFetchNameValue(papszTO,"SRC_SRS"); + + hSrcDS = GDALOpenEx( papszSrcFiles[iSrc], GDAL_OF_RASTER, NULL, NULL, NULL ); + if( hSrcDS == NULL ) + GDALExit( 1 ); + + /* -------------------------------------------------------------------- */ + /* Check that there's at least one raster band */ + /* -------------------------------------------------------------------- */ + if ( GDALGetRasterCount(hSrcDS) == 0 ) + { + fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] ); + GDALExit( 1 ); + } + + if( eDT == GDT_Unknown ) { + eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1)); + } + + /* -------------------------------------------------------------------- */ + /* If we are processing the first file, and it has a color */ + /* table, then we will copy it to the destination file. */ + /* -------------------------------------------------------------------- */ + if( iSrc == 0 ) + { + nDstBandCount = GDALGetRasterCount(hSrcDS); + hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) ); + if( hCT != NULL ) + { + hCT = GDALCloneColorTable( hCT ); + printf( "Copying color table from %s to new file.\n", papszSrcFiles[iSrc] ); + } + + for(int iBand = 0; iBand < nDstBandCount; iBand++) + { + apeColorInterpretations.push_back( + GDALGetRasterColorInterpretation(GDALGetRasterBand(hSrcDS,iBand+1)) ); + } + } + + /* -------------------------------------------------------------------- */ + /* Get the sourcesrs from the dataset, if not set already. */ + /* -------------------------------------------------------------------- */ + if( pszThisSourceSRS == NULL ) + { + const char *pszMethod = CSLFetchNameValue( papszTO, "METHOD" ); + + if( GDALGetProjectionRef( hSrcDS ) != NULL + && strlen(GDALGetProjectionRef( hSrcDS )) > 0 + && (pszMethod == NULL || EQUAL(pszMethod,"GEOTRANSFORM")) ) { + + pszThisSourceSRS = GDALGetProjectionRef( hSrcDS ); + } else if( GDALGetGCPProjection( hSrcDS ) != NULL + && strlen(GDALGetGCPProjection(hSrcDS)) > 0 + && GDALGetGCPCount( hSrcDS ) > 1 + && (pszMethod == NULL || EQUALN(pszMethod,"GCP_",4)) ) { + pszThisSourceSRS = GDALGetGCPProjection( hSrcDS ); + } else if( pszMethod != NULL && EQUAL(pszMethod,"RPC") ) { + pszThisSourceSRS = SRS_WKT_WGS84; + } else { + pszThisSourceSRS = ""; + } + } + + if( pszThisTargetSRS == NULL ) { + pszThisTargetSRS = CPLStrdup( pszThisSourceSRS ); + } + + GDALClose( hSrcDS ); + } + + /* -------------------------------------------------------------------- */ + /* Did we have any usable sources? */ + /* -------------------------------------------------------------------- */ + if( nDstBandCount == 0 ) + { + CPLError( CE_Failure, CPLE_AppDefined, + "No usable source images." ); + CPLFree( pszThisTargetSRS ); + return NULL; + } + + /* -------------------------------------------------------------------- */ + /* Turn the suggested region into a geotransform and suggested */ + /* number of pixels and lines. */ + /* -------------------------------------------------------------------- */ + double adfDstGeoTransform[6] = { 0, 0, 0, 0, 0, 0 }; + int nPixels = 0, nLines = 0; + + /* -------------------------------------------------------------------- */ + /* Did the user override some parameters? */ + /* -------------------------------------------------------------------- */ + if( nForcePixels != 0 && nForceLines != 0 ) + { + if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 ) + { + dfMinX = dfWrkMinX; + dfMaxX = dfWrkMaxX; + dfMaxY = dfWrkMaxY; + dfMinY = dfWrkMinY; + } + + dfXRes = (dfMaxX - dfMinX) / nForcePixels; + dfYRes = (dfMaxY - dfMinY) / nForceLines; + + adfDstGeoTransform[0] = dfMinX; + adfDstGeoTransform[3] = dfMaxY; + adfDstGeoTransform[1] = dfXRes; + adfDstGeoTransform[5] = -dfYRes; + + nPixels = nForcePixels; + nLines = nForceLines; + } + else + { + fprintf(stderr, "UHOH - need force pixels/lines\n"); + } + + /* -------------------------------------------------------------------- */ + /* Create the output file. */ + /* -------------------------------------------------------------------- */ + char** papszCreateOptions = CSLAddString( NULL, "COMPRESS=DEFLATE" ); + hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines, + nDstBandCount, eDT, papszCreateOptions ); + + CSLDestroy( papszCreateOptions ); + papszCreateOptions = NULL; + + if( hDstDS == NULL ) + { + printf( "Error creating file - return NULL\n" ); + CPLFree( pszThisTargetSRS ); + return NULL; + } + + /* -------------------------------------------------------------------- */ + /* Write out the projection definition. */ + /* -------------------------------------------------------------------- */ + const char *pszDstMethod = CSLFetchNameValue(papszTO,"DST_METHOD"); + if( pszDstMethod == NULL || !EQUAL(pszDstMethod, "NO_GEOTRANSFORM") ) + { + if( GDALSetProjection( hDstDS, pszThisTargetSRS ) == CE_Failure || + GDALSetGeoTransform( hDstDS, adfDstGeoTransform ) == CE_Failure ) + { + printf( "Set projection of hDstDS - error\n" ); + CPLFree( pszThisTargetSRS ); + return NULL; + } + } + else + { + adfDstGeoTransform[0] = 0.0; + adfDstGeoTransform[3] = 0.0; + adfDstGeoTransform[5] = fabs(adfDstGeoTransform[5]); + } + + /* -------------------------------------------------------------------- */ + /* Copy the color table, if required. */ + /* -------------------------------------------------------------------- */ + if( hCT != NULL ) + { + printf( "copy color table\n" ); + GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT ); + GDALDestroyColorTable( hCT ); + } + + printf( "free targetSRS\n" ); + CPLFree( pszThisTargetSRS ); + + printf( "GDALWarpCreateOutput complete\n" ); + + return hDstDS; +} + +void SGDemTile::doWarp( int iSrc, char* pszSrcFile, GDALDatasetH hDstDS, char** papszTO, char** papszWarpOptions ) +{ + GDALDatasetH hSrcDS; + GDALResampleAlg eResampleAlg = GRA_NearestNeighbour; + int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE; + GDALDataType eWorkingType = GDT_Unknown; + void* hTransformArg = NULL; + GDALTransformerFunc pfnTransformer = NULL; + double dfErrorThreshold = 0.125; + + /* -------------------------------------------------------------------- */ + /* Open this file. */ + /* -------------------------------------------------------------------- */ + hSrcDS = GDALOpenEx( pszSrcFile, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, NULL, NULL, NULL ); + if( hSrcDS == NULL ) + GDALExit( 2 ); + + /* -------------------------------------------------------------------- */ + /* Check that there's at least one raster band */ + /* -------------------------------------------------------------------- */ + if ( GDALGetRasterCount(hSrcDS) == 0 ) + { + fprintf(stderr, "Input file %s has no raster bands.\n", pszSrcFile ); + GDALExit( 1 ); + } + + printf( "Processing input file %s.\n", pszSrcFile ); + +#if 0 // do we need metadata? + /* -------------------------------------------------------------------- */ + /* Get the metadata of the first source DS and copy it to the */ + /* destination DS. Copy Band-level metadata and other info, only */ + /* if source and destination band count are equal. Any values that */ + /* conflict between source datasets are set to pszMDConflictValue. */ + /* -------------------------------------------------------------------- */ + if ( true ) + { + char **papszMetadata = NULL; + const char *pszSrcInfo = NULL; + const char *pszDstInfo = NULL; + GDALRasterBandH hSrcBand = NULL; + GDALRasterBandH hDstBand = NULL; + + /* copy metadata from first dataset */ + if ( iSrc == 0 ) + { + CPLDebug("WARP", "Copying metadata from first source to destination dataset"); + /* copy dataset-level metadata */ + papszMetadata = GDALGetMetadata( hSrcDS, NULL ); + + char** papszMetadataNew = NULL; + for( int i = 0; papszMetadata != NULL && papszMetadata[i] != NULL; i++ ) + { + // Do not preserve NODATA_VALUES when the output includes an alpha band + if( bEnableDstAlpha && + EQUALN(papszMetadata[i], "NODATA_VALUES=", strlen("NODATA_VALUES=")) ) + { + continue; + } + + papszMetadataNew = CSLAddString(papszMetadataNew, papszMetadata[i]); + } + + if ( CSLCount(papszMetadataNew) > 0 ) { + if ( GDALSetMetadata( hDstDS, papszMetadataNew, NULL ) != CE_None ) + fprintf( stderr, "Warning: error copying metadata to destination dataset.\n" ); + } + + CSLDestroy(papszMetadataNew); + + /* copy band-level metadata and other info */ + if ( GDALGetRasterCount( hSrcDS ) == GDALGetRasterCount( hDstDS ) ) + { + for ( int iBand = 0; iBand < GDALGetRasterCount( hSrcDS ); iBand++ ) + { + hSrcBand = GDALGetRasterBand( hSrcDS, iBand + 1 ); + hDstBand = GDALGetRasterBand( hDstDS, iBand + 1 ); + /* copy metadata, except stats (#5319) */ + papszMetadata = GDALGetMetadata( hSrcBand, NULL); + if ( CSLCount(papszMetadata) > 0 ) + { + //GDALSetMetadata( hDstBand, papszMetadata, NULL ); + char** papszMetadataNew = NULL; + for( int i = 0; papszMetadata != NULL && papszMetadata[i] != NULL; i++ ) + { + if (strncmp(papszMetadata[i], "STATISTICS_", 11) != 0) + papszMetadataNew = CSLAddString(papszMetadataNew, papszMetadata[i]); + } + GDALSetMetadata( hDstBand, papszMetadataNew, NULL ); + CSLDestroy(papszMetadataNew); + } + /* copy other info (Description, Unit Type) - what else? */ + if ( bCopyBandInfo ) { + pszSrcInfo = GDALGetDescription( hSrcBand ); + if( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 ) + GDALSetDescription( hDstBand, pszSrcInfo ); + pszSrcInfo = GDALGetRasterUnitType( hSrcBand ); + if( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 ) + GDALSetRasterUnitType( hDstBand, pszSrcInfo ); + } + } + } + } + /* remove metadata that conflicts between datasets */ + else + { + CPLDebug("WARP", "Removing conflicting metadata from destination dataset (source #%d)", iSrc ); + /* remove conflicting dataset-level metadata */ + RemoveConflictingMetadata( hDstDS, GDALGetMetadata( hSrcDS, NULL ), pszMDConflictValue ); + + /* remove conflicting copy band-level metadata and other info */ + if ( GDALGetRasterCount( hSrcDS ) == GDALGetRasterCount( hDstDS ) ) + { + for ( int iBand = 0; iBand < GDALGetRasterCount( hSrcDS ); iBand++ ) + { + hSrcBand = GDALGetRasterBand( hSrcDS, iBand + 1 ); + hDstBand = GDALGetRasterBand( hDstDS, iBand + 1 ); + /* remove conflicting metadata */ + RemoveConflictingMetadata( hDstBand, GDALGetMetadata( hSrcBand, NULL ), pszMDConflictValue ); + /* remove conflicting info */ + if ( bCopyBandInfo ) { + pszSrcInfo = GDALGetDescription( hSrcBand ); + pszDstInfo = GDALGetDescription( hDstBand ); + if( ! ( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 && + pszDstInfo != NULL && strlen(pszDstInfo) > 0 && + EQUAL( pszSrcInfo, pszDstInfo ) ) ) + GDALSetDescription( hDstBand, "" ); + pszSrcInfo = GDALGetRasterUnitType( hSrcBand ); + pszDstInfo = GDALGetRasterUnitType( hDstBand ); + if( ! ( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 && + pszDstInfo != NULL && strlen(pszDstInfo) > 0 && + EQUAL( pszSrcInfo, pszDstInfo ) ) ) + GDALSetRasterUnitType( hDstBand, "" ); + } + } + } + } + } +#endif + + /* -------------------------------------------------------------------- */ + /* Warns if the file has a color table and something more */ + /* complicated than nearest neighbour resampling is asked */ + /* -------------------------------------------------------------------- */ + if ( eResampleAlg != GRA_NearestNeighbour && eResampleAlg != GRA_Mode && + GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL) + { + fprintf( stderr, "Warning: Input file %s has a color table, which will likely lead to " + "bad results when using a resampling method other than " + "nearest neighbour or mode. Converting the dataset prior to 24/32 bit " + "is advised.\n", pszSrcFile ); + } + + /* -------------------------------------------------------------------- */ + /* Do we have a source alpha band? */ + /* -------------------------------------------------------------------- */ + if( GDALGetRasterColorInterpretation( GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) ) == GCI_AlphaBand && !bEnableSrcAlpha ) + { + printf( "SHOULD NOT HAPPEN Using band %d of source image as alpha.\n", GDALGetRasterCount(hSrcDS) ); + bEnableSrcAlpha = TRUE; + } + + /* -------------------------------------------------------------------- */ + /* Create a transformation object from the source to */ + /* destination coordinate system. */ + /* -------------------------------------------------------------------- */ + hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, papszTO ); + if( hTransformArg == NULL ) { + printf( "SHOULD NOT HAPPEN hTransformArg is NULL\n" ); + GDALExit( 1 ); + } + + /* -------------------------------------------------------------------- */ + /* Warp the transformer with a linear approximator */ + /* -------------------------------------------------------------------- */ + hTransformArg = GDALCreateApproxTransformer( GDALGenImgProjTransform, + hTransformArg, dfErrorThreshold); + pfnTransformer = GDALApproxTransform; + GDALApproxTransformerOwnsSubtransformer(hTransformArg, TRUE); + + /* -------------------------------------------------------------------- */ + /* Clear temporary INIT_DEST settings after the first image. */ + /* -------------------------------------------------------------------- */ + if( iSrc == 1 ) + papszWarpOptions = CSLSetNameValue( papszWarpOptions, + "INIT_DEST", NULL ); + + /* -------------------------------------------------------------------- */ + /* Setup warp options. */ + /* -------------------------------------------------------------------- */ + GDALWarpOptions *psWO = GDALCreateWarpOptions(); + + psWO->papszWarpOptions = CSLDuplicate(papszWarpOptions); + psWO->eWorkingDataType = eWorkingType; + psWO->eResampleAlg = eResampleAlg; + + psWO->hSrcDS = hSrcDS; + psWO->hDstDS = hDstDS; + + psWO->pfnTransformer = pfnTransformer; + psWO->pTransformerArg = hTransformArg; + + psWO->pfnProgress = GDALTermProgress; + //psWO->pfnProgress = NULL; + + /* -------------------------------------------------------------------- */ + /* Setup band mapping. */ + /* -------------------------------------------------------------------- */ + if( bEnableSrcAlpha ) + psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1; + else + psWO->nBandCount = GDALGetRasterCount(hSrcDS); + + psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int)); + psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int)); + + for( int i = 0; i < psWO->nBandCount; i++ ) + { + psWO->panSrcBands[i] = i+1; + psWO->panDstBands[i] = i+1; + } + + /* -------------------------------------------------------------------- */ + /* Setup alpha bands used if any. */ + /* -------------------------------------------------------------------- */ + if( bEnableSrcAlpha ) + psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS); + + if( !bEnableDstAlpha + && GDALGetRasterCount(hDstDS) == psWO->nBandCount+1 + && GDALGetRasterColorInterpretation( + GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS))) + == GCI_AlphaBand ) + { + printf( "Using band %d of destination image as alpha.\n", + GDALGetRasterCount(hDstDS) ); + + bEnableDstAlpha = TRUE; + } + + if( bEnableDstAlpha ) + psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS); + + int bHaveNodata = FALSE; + double dfReal = 0.0; + + for( int i = 0; !bHaveNodata && i < psWO->nBandCount; i++ ) + { + GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 ); + dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata ); + } + + if( bHaveNodata ) + { + if (CPLIsNan(dfReal)) + printf( "Using internal nodata values (e.g. nan) for image %s.\n", + pszSrcFile ); + else + printf( "Using internal nodata values (e.g. %g) for image %s.\n", + dfReal, pszSrcFile ); + + psWO->padfSrcNoDataReal = (double *) + CPLMalloc(psWO->nBandCount*sizeof(double)); + psWO->padfSrcNoDataImag = (double *) + CPLMalloc(psWO->nBandCount*sizeof(double)); + + for( int i = 0; i < psWO->nBandCount; i++ ) + { + GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 ); + + dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata ); + + if( bHaveNodata ) + { + psWO->padfSrcNoDataReal[i] = dfReal; + psWO->padfSrcNoDataImag[i] = 0.0; + } + else + { + psWO->padfSrcNoDataReal[i] = -123456.789; + psWO->padfSrcNoDataImag[i] = 0.0; + } + } + } + + /* else try to fill dstNoData from source bands */ + if ( psWO->padfSrcNoDataReal != NULL ) + { + psWO->padfDstNoDataReal = (double *) + CPLMalloc(psWO->nBandCount*sizeof(double)); + psWO->padfDstNoDataImag = (double *) + CPLMalloc(psWO->nBandCount*sizeof(double)); + + printf( "Copying nodata values from source %s \n", pszSrcFile ); + + for( int i = 0; i < psWO->nBandCount; i++ ) + { + int bHaveNodata = FALSE; + + GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 ); + GDALGetRasterNoDataValue( hBand, &bHaveNodata ); + + CPLDebug("WARP", "band=%d bHaveNodata=%d", i, bHaveNodata); + if( bHaveNodata ) + { + psWO->padfDstNoDataReal[i] = psWO->padfSrcNoDataReal[i]; + psWO->padfDstNoDataImag[i] = psWO->padfSrcNoDataImag[i]; + CPLDebug("WARP", "srcNoData=%f dstNoData=%f", + psWO->padfSrcNoDataReal[i], psWO->padfDstNoDataReal[i] ); + } + + CPLDebug("WARP", "calling GDALSetRasterNoDataValue() for band#%d", i ); + GDALSetRasterNoDataValue( + GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ), + psWO->padfDstNoDataReal[i] ); + } + } + + /* -------------------------------------------------------------------- */ + /* Initialize and execute the warp. */ + /* -------------------------------------------------------------------- */ + GDALWarpOperation oWO; + if( oWO.Initialize( psWO ) == CE_None ) + { + oWO.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) ); + } + + /* -------------------------------------------------------------------- */ + /* Cleanup */ + /* -------------------------------------------------------------------- */ + if( hTransformArg != NULL ) + GDALDestroyTransformer( hTransformArg ); + + GDALDestroyWarpOptions( psWO ); + GDALClose( hSrcDS ); +} diff --git a/simgear/scene/dem/SGMesh.cxx b/simgear/scene/dem/SGMesh.cxx new file mode 100644 index 00000000..50e1eca8 --- /dev/null +++ b/simgear/scene/dem/SGMesh.cxx @@ -0,0 +1,414 @@ +/* +Szymon Rusinkiewicz +Princeton University + +TriMesh_normals.cc +Compute per-vertex normals for TriMeshes + +For meshes, uses average of per-face normals, weighted according to: + Max, N. + "Weights for Computing Vertex Normals from Facet Normals," + Journal of Graphics Tools, Vol. 4, No. 2, 1999. + +For raw point clouds, fits plane to k nearest neighbors. +*/ + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "SGDemSession.hxx" +#include "SGMesh.hxx" + +#define GW_LVL1 (153) +#define SKIP_MULT_LVL1 (1) + +#define GW_LVL2 (78) +#define SKIP_MULT_LVL2 (2) + +#define GW_LVL3 (53) +#define SKIP_MULT_LVL3 (3) + +#define GW_LVL4 (33) +#define SKIP_MULT_LVL4 (5) + +#define GW_LVL5 (18) +#define SKIP_MULT_LVL5 (10) + +#define GW_LVL6 (13) +#define SKIP_MULT_LVL6 (15) + +#define GW_LVL7 (9) +#define SKIP_MULT_LVL7 (25) + +#define GW_TEST GW_LVL2 +#define SKIP_TEST SKIP_MULT_LVL2 + + +using namespace std; + +SGMesh::SGMesh( const SGDemPtr dem, + unsigned wo, unsigned so, + unsigned eo, unsigned no, + unsigned heightLevel, + unsigned widthLevel, + const osg::Matrixd& transform, + TextureMethod tm, + const osgDB::Options* options ) +{ + ::std::vector skirt_geods; + ::std::vector index; + + flag_curr = 0; + + int lvl = -1; + int skipx = -1, skipy = -1; + + bool Debug1 = false; + bool Debug2 = false; + + if ( wo == 1444 && eo == 1446 && so == 1055 && no == 1056 ) { + Debug1 = true; + } else if ( wo == 1444 && eo == 1446 && so == 1056 && no == 1057 ) { + Debug2 = true; + } + + switch ( heightLevel ) { + case 9: lvl = 0; grid_height = 153; skipy = 1; + // at highest lod for height, tiles have different widths based on latitude + switch( widthLevel ) { + case 9: grid_width = 153; skipx = 1; break; // 1/8 x 1/8 + case 8: grid_width = 153; skipx = 2; break; // 1/4 x 1/8 + case 7: grid_width = 153; skipx = 4; break; // 1/2 x 1/8 + case 6: grid_width = 153; skipx = 8; break; // 1 x 1/8 + case 5: grid_width = 153; skipx = 16; break; // 2 x 1/8 + case 4: grid_width = 153; skipx = 32; break; // 4 x 1/8 + case 3: grid_width = 78; skipx = 32; break; // 12 x 1/8 + } + break; + + case 8: lvl = 0; grid_height = GW_TEST; skipy = 2*SKIP_TEST; grid_width = GW_TEST; skipx = 2*SKIP_TEST; break; // 1/4 x 1/4 + case 7: lvl = 0; grid_height = GW_TEST; skipy = 4*SKIP_TEST; grid_width = GW_TEST; skipx = 4*SKIP_TEST; break; // 1/2 x 1/2 + case 6: lvl = 0; grid_height = GW_TEST; skipy = 8*SKIP_TEST; grid_width = GW_TEST; skipx = 8*SKIP_TEST; break; // 1 x 1 + + case 5: lvl = 1; grid_height = GW_TEST; skipy = 1*SKIP_TEST; grid_width = GW_TEST; skipx = 1*SKIP_TEST; break; // 2 x 2 + case 4: lvl = 1; grid_height = GW_TEST; skipy = 2*SKIP_TEST; grid_width = GW_TEST; skipx = 2*SKIP_TEST; break; // 4 x 4 + case 3: lvl = 1; grid_height = GW_TEST; skipy = 6*SKIP_TEST; + switch ( widthLevel ) { + case 3: grid_width = GW_TEST; skipx = 6*SKIP_TEST; break; // 12 x 12 + case 2: grid_width = GW_TEST; skipx = 18*SKIP_TEST; break; // 36 x 12 + } + break; + + case 2: lvl = 2; grid_height = GW_TEST; skipy = 3*SKIP_TEST; grid_width = GW_TEST; skipx = 1*SKIP_TEST; break; // 180 x 60 + case 1: lvl = 2; grid_height = GW_TEST; skipy = 6*SKIP_TEST; grid_width = GW_TEST; skipx = 3*SKIP_TEST; break; // 360 x 180 + + case 0: printf("ERROR - no height level 0\n"); exit(0); break; + } + + double incu = ( eo*1.0/(360*8) - wo*1.0/(360*8) ) / (grid_width - 3); + double incv = ( no*1.0/(180*8) - so*1.0/(180*8) ) / (grid_height - 3); + + // try to use native mesh res + vertices = new osg::Vec3Array( grid_width*grid_height ); + normals = NULL; + texCoords = new osg::Vec2Array( grid_width*grid_height ); + + ::std::vector geodes(grid_width*grid_height); + + // session can't be paralell yet - save alts in geode array + fprintf( stderr, "SGMesh::SGMesh - create session - num dem roots is %d\n", dem->getNumRoots() ); + SGDemSession s = dem->openSession( wo, so, eo, no, lvl, true ); + s.getGeods( wo, so, eo, no, grid_width, grid_height, skipx, skipy, geodes, Debug1, Debug2 ); + s.close(); + + // save the west skirt vertices + unsigned int src_idx = 0; + unsigned int edge_idx = 0; + + // save west skirt + for ( edge_idx = 0, src_idx = 0; edge_idx < grid_height; edge_idx++, src_idx++ ) { + skirt_geods.push_back( geodes[src_idx] ); + index.push_back( src_idx ); + } + + // save the north skirt vertices + for ( edge_idx = 1, src_idx = (grid_height*2)-1; edge_idx < grid_width; edge_idx++, src_idx += grid_height ) { + skirt_geods.push_back( geodes[src_idx] ); + index.push_back( src_idx ); + } + + // save the east skirt vertices + for ( edge_idx = 0, src_idx = grid_height*(grid_width-1); edge_idx < grid_height-1; edge_idx++, src_idx++ ) { + skirt_geods.push_back( geodes[src_idx] ); + index.push_back( src_idx ); + } + + // save the south skirt vertices + for ( edge_idx = 1, src_idx = grid_height; edge_idx < grid_width-1; src_idx += grid_height, edge_idx++ ) { + skirt_geods.push_back( geodes[src_idx] ); + index.push_back( src_idx ); + } + + // we can convert to cartesian in paralell + unsigned int nv = geodes.size(); +#pragma omp parallel for + for (unsigned int i = 0; i < nv; i++) { + (*vertices)[i].set( toOsg( SGVec3f::fromGeod( geodes[i] ) ) ); + } + + need_faces(); + need_normals(); + + // lower skirt verts based on lvl and skip + unsigned int lower = 0; + switch( lvl ) { + case 0: + lower = skipx*100; + break; + + case 1: + lower = skipx*10000; + break; + + case 2: + lower = skipx*1000000; + break; + } + + for ( unsigned int i=0; isize(); v++ ) { + colors->push_back(lvlColor); + } + } else if ( tm == SGMesh::TEXTURE_BLUEMARBLE ) { + colors->push_back(osg::Vec4(1, 1, 1, 1)); + } else if ( tm == SGMesh::TEXTURE_RASTER ) { + // TODO + } + + geode = new osg::Geode; + osg::Geometry* geometry = new osg::Geometry; + + char geoName[64]; + snprintf( geoName, sizeof(geoName), "tilemesh (%u,%u)-(%u,%u):level%d,%d", + wo, so, eo, no, + widthLevel, heightLevel ); + geometry->setName(geoName); + + geometry->setDataVariance(osg::Object::STATIC); + geometry->setUseVertexBufferObjects(true); + geometry->setVertexArray( getVertices() ); + + geometry->setNormalArray( getNormals() ); + geometry->setNormalBinding(osg::Geometry::BIND_PER_VERTEX); + geometry->setColorArray(colors); + + if ( tm == SGMesh::TEXTURE_DEBUG ) { + geometry->setColorBinding(osg::Geometry::BIND_PER_VERTEX); + } else { + geometry->setColorBinding(osg::Geometry::BIND_OVERALL); + geometry->setTexCoordArray(0, getTexCoords() ); + } + + // generate triangles from the grid + osg::DrawElementsUInt* indices = getIndices(); + geometry->addPrimitiveSet(indices); + + if ( geometry ) { + geode->setDataVariance(osg::Object::STATIC); + geode->addDrawable(geometry); + + if ( tm == SGMesh::TEXTURE_DEBUG ) { + // set up the state + osg::StateSet* stateSet = geode->getOrCreateStateSet(); + stateSet->setRenderBinDetails(-10, "RenderBin"); + + osg::ShadeModel* shadeModel = new osg::ShadeModel; + shadeModel->setMode(osg::ShadeModel::SMOOTH); + stateSet->setAttributeAndModes(shadeModel); + + stateSet->setMode(GL_LIGHTING, osg::StateAttribute::ON); + stateSet->setMode(GL_FOG, osg::StateAttribute::OFF); + stateSet->setMode(GL_DEPTH_TEST, osg::StateAttribute::ON); + stateSet->setMode(GL_CULL_FACE, osg::StateAttribute::ON); + stateSet->setMode(GL_BLEND, osg::StateAttribute::OFF); + stateSet->setMode(GL_ALPHA_TEST, osg::StateAttribute::OFF); + + stateSet->setAttribute(new osg::CullFace(osg::CullFace::BACK)); + + osg::Material* material = new osg::Material; + material->setColorMode(osg::Material::DIFFUSE); + material->setDiffuse(osg::Material::FRONT_AND_BACK, + osg::Vec4(1, 1, 1, 1)); + material->setAmbient(osg::Material::FRONT_AND_BACK, + osg::Vec4(0, 0, 0, 1)); + material->setEmission(osg::Material::FRONT_AND_BACK, + osg::Vec4(0, 0, 0, 1)); + material->setSpecular(osg::Material::FRONT_AND_BACK, + osg::Vec4(0, 0, 0, 1)); + material->setShininess(osg::Material::FRONT_AND_BACK, 0); + stateSet->setAttribute(material); + + geode->setStateSet(stateSet); + } else if ( tm == SGMesh::TEXTURE_BLUEMARBLE ) { + osg::StateSet* stateSet = new osg::StateSet; + stateSet->setAttributeAndModes(new osg::CullFace); + + std::string imageFileName = options->getPluginStringData("SimGear::FG_WORLD_TEXTURE"); + if (imageFileName.empty()) { + imageFileName = options->getPluginStringData("SimGear::FG_ROOT"); + imageFileName = osgDB::concatPaths(imageFileName, "Textures"); + imageFileName = osgDB::concatPaths(imageFileName, "Globe"); + imageFileName = osgDB::concatPaths(imageFileName, "world.topo.bathy.200407.3x4096x2048.png"); + } + if (osg::Image* image = osgDB::readImageFile(imageFileName, options)) { + osg::Texture2D* texture = new osg::Texture2D; + texture->setImage(image); + texture->setWrap(osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT); + texture->setWrap(osg::Texture2D::WRAP_T, osg::Texture2D::CLAMP); + stateSet->setTextureAttributeAndModes(0, texture); + } + + geode->setStateSet(stateSet); + } else if ( tm == SGMesh::TEXTURE_RASTER ) { + // TODO + } else { + printf("texture unknown!\n"); + exit(0); + } + } +} + +void SGMesh::need_faces() +{ + for (unsigned int i = 0; i < grid_width-1; i++) { + for (unsigned int j = 0; j < grid_height-1; j++) { + // for each point - create 2 triangles with current point in the + // lower left corner + + // 0,0 - 1,0 - 0,1 + faces.push_back( Face( (i+0)*grid_height + (j+0), + (i+1)*grid_height + (j+0), + (i+0)*grid_height + (j+1) ) ); + + // 1,1 - 0,1 - 1,0 + faces.push_back( Face( (i+1)*grid_height + (j+1), + (i+0)*grid_height + (j+1), + (i+1)*grid_height + (j+0) ) ); + } + } +} + +// Compute per-vertex normals +void SGMesh::need_normals() +{ + if ( normals == NULL ) + { + normals = new osg::Vec3Array( vertices->size() ); + + // need_faces(); + if ( !faces.empty() ) { + // Compute from faces + int nf = faces.size(); +#pragma omp parallel for + for (int i = 0; i < nf; i++) { + const osg::Vec3 &p0 = (*vertices)[faces[i][0]]; + const osg::Vec3 &p1 = (*vertices)[faces[i][1]]; + const osg::Vec3 &p2 = (*vertices)[faces[i][2]]; + + osg::Vec3 a = p0-p1, b = p1-p2, c = p2-p0; + float l2a = a.length2(), l2b = b.length2(), l2c = c.length2(); + if (!l2a || !l2b || !l2c) { + continue; + } + osg::Vec3 facenormal = a ^ b; + (*normals)[faces[i][0]] += facenormal * (1.0f / (l2a * l2c)); + (*normals)[faces[i][1]] += facenormal * (1.0f / (l2b * l2a)); + (*normals)[faces[i][2]] += facenormal * (1.0f / (l2c * l2b)); + } + } + + // Make them all unit-length + unsigned int nn = normals->size(); +#pragma omp parallel for + for (unsigned int i = 0; i < nn; i++) + (*normals)[i].normalize(); + } +} + +osg::DrawElementsUInt* SGMesh::getIndices(void) +{ + osg::DrawElementsUInt* indices = new osg::DrawElementsUInt(osg::PrimitiveSet::TRIANGLES, 0); + + for (unsigned int i = 0; i < grid_width-1; i++) { + for (unsigned int j = 0; j < grid_height-1; j++) { + // for each point - create 2 triangles with current point in the + // lower left corner + indices->push_back( (i+0)*grid_height + (j+0) ); // 0,0 + indices->push_back( (i+1)*grid_height + (j+0) ); // 1,0 + indices->push_back( (i+0)*grid_height + (j+1) ); // 0,1 + + indices->push_back( (i+1)*grid_height + (j+1) ); // 1,1 + indices->push_back( (i+0)*grid_height + (j+1) ); // 0,1 + indices->push_back( (i+1)*grid_height + (j+0) ); // 1,0 + } + } + indices->setDataVariance(osg::Object::STATIC); + + return indices; +} diff --git a/simgear/scene/dem/SGMesh.hxx b/simgear/scene/dem/SGMesh.hxx new file mode 100644 index 00000000..3f9ee143 --- /dev/null +++ b/simgear/scene/dem/SGMesh.hxx @@ -0,0 +1,718 @@ +// SGMesh.hxx -- mesh normals and curvature from a regular grid +// +// Copyright (C) 2016 Peter Sadrozinski +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// + +#ifndef _SGMESH_HXX +#define _SGMESH_HXX + +/* +The Simgear Mesh object is based on the GPLv2 library trimesh by + +Szymon Rusinkiewicz +Princeton University + +TriMesh.h +Class for triangle meshes. + +from http://gfx.cs.princeton.edu/gfx/proj/trimesh2 + +Based on a paper: http://gfx.cs.princeton.edu/pubs/_2004_ECA/curvpaper.pdf + +This class has been reduced in scope, and the native types have been +converted to OpenSceneGraph for use in SimGear/FlightGear + +*/ + +#include +#include + +#include +#include +#include + +class KDtree { +private: + class Node; + Node *root; + void build(const float *ptlist, size_t n); + +public: + // Compatibility function for closest-compatible-point searches + struct CompatFunc + { + virtual bool operator () (const float *p) const = 0; + virtual ~CompatFunc() {} // To make the compiler shut up + }; + + // Constructor from an array of points + KDtree(const float *ptlist, size_t n) + { build(ptlist, n); } + + // Constructor from a vector of points + template KDtree(const ::std::vector &v) + { build((const float *) &v[0], v.size()); } + + // Destructor - recursively frees the tree + ~KDtree(); + + // The queries: returns closest point to a point or a ray, + // provided it's within sqrt(maxdist2) and is compatible + const float *closest_to_pt(const float *p, + float maxdist2 = 0.0f, + const CompatFunc *iscompat = NULL) const; + const float *closest_to_ray(const float *p, const float *dir, + float maxdist2 = 0.0f, + const CompatFunc *iscompat = NULL) const; + + // Find the k nearest neighbors + void find_k_closest_to_pt(::std::vector &knn, + int k, + const float *p, + float maxdist2 = 0.0f, + const CompatFunc *iscompat = NULL) const; +}; + +// Windows defines min and max as macros, which prevents us from +// using the type-safe versions from std:: +// Also define NOMINMAX, which prevents future bad definitions. +#ifdef min +# undef min +#endif +#ifdef max +# undef max +#endif +#ifndef NOMINMAX +# define NOMINMAX +#endif +#ifndef _USE_MATH_DEFINES + #define _USE_MATH_DEFINES +#endif + +#include +#include +#include + +// LU decomposition +template +static inline bool ludcmp(T a[N][N], int indx[N], T *d = NULL) +{ + using namespace ::std; + T vv[N]; + + if (d) + *d = T(1); + for (int i = 0; i < N; i++) { + T big = T(0); + for (int j = 0; j < N; j++) { + T tmp = fabs(a[i][j]); + if (tmp > big) + big = tmp; + } + if (big == T(0)) + return false; + vv[i] = T(1) / big; + } + for (int j = 0; j < N; j++) { + for (int i = 0; i < j; i++) { + T sum = a[i][j]; + for (int k = 0; k < i; k++) + sum -= a[i][k]*a[k][j]; + a[i][j]=sum; + } + T big = T(0); + int imax = j; + for (int i = j; i < N; i++) { + T sum = a[i][j]; + for (int k = 0; k < j; k++) + sum -= a[i][k]*a[k][j]; + a[i][j] = sum; + T tmp = vv[i] * fabs(sum); + if (tmp > big) { + big = tmp; + imax = i; + } + } + if (imax != j) { + for (int k = 0; k < N; k++) + swap(a[imax][k], a[j][k]); + if (d) + *d = -(*d); + vv[imax] = vv[j]; + } + indx[j] = imax; + if (unlikely(a[j][j] == T(0))) + return false; + if (j != N-1) { + T tmp = T(1) / a[j][j]; + for (int i = j+1; i < N; i++) + a[i][j] *= tmp; + } + } + return true; +} + + +// Backsubstitution after ludcmp +template +static inline void lubksb(T a[N][N], int indx[N], T b[N]) +{ + int ii = -1; + for (int i = 0; i < N; i++) { + int ip = indx[i]; + T sum = b[ip]; + b[ip] = b[i]; + if (ii != -1) + for (int j = ii; j < i; j++) + sum -= a[i][j] * b[j]; + else if (sum) + ii = i; + b[i] = sum; + } + for (int i = N-1; i >= 0; i--) { + T sum = b[i]; + for (int j = i+1; j < N; j++) + sum -= a[i][j] * b[j]; + b[i] = sum / a[i][i]; + } +} + + +// Perform LDL^T decomposition of a symmetric positive definite matrix. +// Like Cholesky, but no square roots. Overwrites lower triangle of matrix. +template +static inline bool ldltdc(T A[N][N], T rdiag[N]) +{ + T v[N-1]; + for (int i = 0; i < N; i++) { + for (int k = 0; k < i; k++) + v[k] = A[i][k] * rdiag[k]; + for (int j = i; j < N; j++) { + T sum = A[i][j]; + for (int k = 0; k < i; k++) + sum -= v[k] * A[j][k]; + if (i == j) { + if (unlikely(sum <= T(0))) + return false; + rdiag[i] = T(1) / sum; + } else { + A[j][i] = sum; + } + } + } + + return true; +} + + +// Solve Ax=B after ldltdc +template +static inline void ldltsl(T A[N][N], T rdiag[N], T B[N], T x[N]) +{ + for (int i = 0; i < N; i++) { + T sum = B[i]; + for (int k = 0; k < i; k++) + sum -= A[i][k] * x[k]; + x[i] = sum * rdiag[i]; + } + for (int i = N - 1; i >= 0; i--) { + T sum = 0; + for (int k = i + 1; k < N; k++) + sum += A[k][i] * x[k]; + x[i] -= sum * rdiag[i]; + } +} + + +// Eigenvector decomposition for real, symmetric matrices, +// a la Bowdler et al. / EISPACK / JAMA +// Entries of d are eigenvalues, sorted smallest to largest. +// A changed in-place to have its columns hold the corresponding eigenvectors. +// Note that A must be completely filled in on input. +template +static inline void eigdc(T A[N][N], T d[N]) +{ + using namespace ::std; + + // Householder + T e[N]; + for (int j = 0; j < N; j++) { + d[j] = A[N-1][j]; + e[j] = T(0); + } + for (int i = N-1; i > 0; i--) { + T scale = T(0); + for (int k = 0; k < i; k++) + scale += fabs(d[k]); + if (scale == T(0)) { + e[i] = d[i-1]; + for (int j = 0; j < i; j++) { + d[j] = A[i-1][j]; + A[i][j] = A[j][i] = T(0); + } + d[i] = T(0); + } else { + T h(0); + T invscale = T(1) / scale; + for (int k = 0; k < i; k++) { + d[k] *= invscale; + h += sqr(d[k]); + } + T f = d[i-1]; + T g = (f > T(0)) ? -sqrt(h) : sqrt(h); + e[i] = scale * g; + h -= f * g; + d[i-1] = f - g; + for (int j = 0; j < i; j++) + e[j] = T(0); + for (int j = 0; j < i; j++) { + f = d[j]; + A[j][i] = f; + g = e[j] + f * A[j][j]; + for (int k = j+1; k < i; k++) { + g += A[k][j] * d[k]; + e[k] += A[k][j] * f; + } + e[j] = g; + } + f = T(0); + T invh = T(1) / h; + for (int j = 0; j < i; j++) { + e[j] *= invh; + f += e[j] * d[j]; + } + T hh = f / (h + h); + for (int j = 0; j < i; j++) + e[j] -= hh * d[j]; + for (int j = 0; j < i; j++) { + f = d[j]; + g = e[j]; + for (int k = j; k < i; k++) + A[k][j] -= f * e[k] + g * d[k]; + d[j] = A[i-1][j]; + A[i][j] = T(0); + } + d[i] = h; + } + } + + for (int i = 0; i < N-1; i++) { + A[N-1][i] = A[i][i]; + A[i][i] = T(1); + T h = d[i+1]; + if (h != T(0)) { + T invh = T(1) / h; + for (int k = 0; k <= i; k++) + d[k] = A[k][i+1] * invh; + for (int j = 0; j <= i; j++) { + T g = T(0); + for (int k = 0; k <= i; k++) + g += A[k][i+1] * A[k][j]; + for (int k = 0; k <= i; k++) + A[k][j] -= g * d[k]; + } + } + for (int k = 0; k <= i; k++) + A[k][i+1] = T(0); + } + for (int j = 0; j < N; j++) { + d[j] = A[N-1][j]; + A[N-1][j] = T(0); + } + A[N-1][N-1] = T(1); + + // QL + for (int i = 1; i < N; i++) + e[i-1] = e[i]; + e[N-1] = T(0); + T f = T(0), tmp = T(0); + const T eps = ::std::numeric_limits::epsilon(); + for (int l = 0; l < N; l++) { + tmp = max(tmp, fabs(d[l]) + fabs(e[l])); + int m = l; + while (m < N) { + if (fabs(e[m]) <= eps * tmp) + break; + m++; + } + if (m > l) { + do { + T g = d[l]; + T p = (d[l+1] - g) / (e[l] + e[l]); + T r = T(hypot(p, T(1))); + if (p < T(0)) + r = -r; + d[l] = e[l] / (p + r); + d[l+1] = e[l] * (p + r); + T dl1 = d[l+1]; + T h = g - d[l]; + for (int i = l+2; i < N; i++) + d[i] -= h; + f += h; + p = d[m]; + T c = T(1), c2 = T(1), c3 = T(1); + T el1 = e[l+1], s = T(0), s2 = T(0); + for (int i = m - 1; i >= l; i--) { + c3 = c2; + c2 = c; + s2 = s; + g = c * e[i]; + h = c * p; + r = T(hypot(p, e[i])); + e[i+1] = s * r; + s = e[i] / r; + c = p / r; + p = c * d[i] - s * g; + d[i+1] = h + s * (c * g + s * d[i]); + for (int k = 0; k < N; k++) { + h = A[k][i+1]; + A[k][i+1] = s * A[k][i] + c * h; + A[k][i] = c * A[k][i] - s * h; + } + } + p = -s * s2 * c3 * el1 * e[l] / dl1; + e[l] = s * p; + d[l] = c * p; + } while (fabs(e[l]) > eps * tmp); + } + d[l] += f; + e[l] = T(0); + } + + // Sort + for (int i = 0; i < N-1; i++) { + int k = i; + T p = d[i]; + for (int j = i+1; j < N; j++) { + if (d[j] < p) { + k = j; + p = d[j]; + } + } + if (k == i) + continue; + d[k] = d[i]; + d[i] = p; + for (int j = 0; j < N; j++) + swap(A[j][i], A[j][k]); + } +} + + +// x <- A * d * A' * b +template +static inline void eigmult(T A[N][N], + T d[N], + T b[N], + T x[N]) +{ + T e[N]; + for (int i = 0; i < N; i++) { + e[i] = T(0); + for (int j = 0; j < N; j++) + e[i] += A[j][i] * b[j]; + e[i] *= d[i]; + } + for (int i = 0; i < N; i++) { + x[i] = T(0); + for (int j = 0; j < N; j++) + x[i] += A[i][j] * e[j]; + } +} + +class SGMesh { +public: + // + // Types + // + typedef enum { + TEXTURE_RASTER, + TEXTURE_BLUEMARBLE, + TEXTURE_DEBUG, + } TextureMethod; + + struct Face { + unsigned int v[3]; + + Face() {} + Face(const int &v0, const int &v1, const int &v2) + { v[0] = v0; v[1] = v1; v[2] = v2; } + Face(const int *v_) + { v[0] = v_[0]; v[1] = v_[1]; v[2] = v_[2]; } +// template explicit Face(const S &x) +// { v[0] = x[0]; v[1] = x[1]; v[2] = x[2]; } + unsigned int &operator[] (int i) { return v[i]; } + const unsigned int &operator[] (int i) const { return v[i]; } + operator const unsigned int * () const { return &(v[0]); } + operator const unsigned int * () { return &(v[0]); } + operator unsigned int * () { return &(v[0]); } + int indexof(unsigned int v_) const + { + return (v[0] == v_) ? 0 : + (v[1] == v_) ? 1 : + (v[2] == v_) ? 2 : -1; + } + }; + +#if 0 + struct BSphere { + point center; + float r; + bool valid; + BSphere() : valid(false) + {} + }; + + // + // Enums + // + enum TstripRep { TSTRIP_LENGTH, TSTRIP_TERM }; + enum { GRID_INVALID = -1 }; + enum StatOp { STAT_MIN, STAT_MAX, STAT_MEAN, STAT_MEANABS, + STAT_RMS, STAT_MEDIAN, STAT_STDEV, STAT_TOTAL }; + enum StatVal { STAT_VALENCE, STAT_FACEAREA, STAT_ANGLE, + STAT_DIHEDRAL, STAT_EDGELEN, STAT_X, STAT_Y, STAT_Z }; +#endif + + // + // Constructor + // + SGMesh() : flag_curr(0) + {} + + SGMesh( const SGDemPtr dem, + unsigned wo, unsigned so, + unsigned eo, unsigned no, + unsigned heightLevel, + unsigned widthLevel, + const osg::Matrixd& transform, + TextureMethod tm, + const osgDB::Options* options + ); + + // + // Members + // + osg::Vec3Array* getVertices(void) { + return vertices; + } + + osg::Vec3Array* getNormals(void) { + return normals; + } + + osg::Vec2Array* getTexCoords(void) { + return texCoords; + } + + osg::Geode* getGeode(void) { + return geode; + } + + osg::DrawElementsUInt* getIndices(void); + + // The basics: vertices and faces + // note: + // vertices, normals, texcoords, etc are allocated by the mesh, + // but are inserted into the scene graph. They are smart pointers + osg::Vec3Array* vertices; + osg::Vec3Array* normals; + osg::Vec2Array* texCoords; + ::std::vector faces; + +// // Triangle strips +// ::std::vector tstrips; + + // Grid, if present +// ::std::vector grid; + unsigned int grid_width, grid_height; + + // Other per-vertex properties +// ::std::vector colors; + ::std::vector confidences; + ::std::vector flags; + unsigned flag_curr; + + // Computed per-vertex properties + ::std::vector pdir1, pdir2; + ::std::vector curv1, curv2; + ::std::vector dcurv; + ::std::vector cornerareas; + ::std::vector pointareas; + + osg::Geode* geode; + + // Bounding structures + //box bbox; + //BSphere bsphere; + + // Connectivity structures: + // For each vertex, all neighboring vertices + ::std::vector< ::std::vector > neighbors; + // For each vertex, all neighboring faces + ::std::vector< ::std::vector > adjacentfaces; + // For each face, the three faces attached to its edges + // (for example, across_edge[3][2] is the number of the face + // that's touching the edge opposite vertex 2 of face 3) + ::std::vector across_edge; + + inline osg::Vec3 trinorm(const osg::Vec3 &v0, const osg::Vec3 &v1, const osg::Vec3 &v2) + { + return ( (v1 - v0) ^ (v2 - v0) ) * 0.5; + } + + static inline const float angle(const osg::Vec3 &v1, const osg::Vec3 &v2) + { + return std::atan2( (v1 ^ v2).length(), (v1 * v2) ); + } + + // + // Compute all this stuff... + // +// void need_tstrips(); +// void convert_strips(TstripRep rep); +// void unpack_tstrips(); +// void triangulate_grid(bool remove_slivers = true); + void need_faces(); +// { +// if (!faces.empty()) +// return; +// if (!tstrips.empty()) +// unpack_tstrips(); +// else if (!grid.empty()) +// triangulate_grid(); +// } + + void need_normals(); + void need_pointareas(); + void need_curvatures(); + void need_dcurv(); +// void need_bbox(); +// void need_bsphere(); + void need_neighbors(); + void need_adjacentfaces(); + void need_across_edge(); + + // + // Delete everything + // + void clear() + { + vertices = NULL; + faces.clear(); +// tstrips.clear(); +// grid.clear(); + grid_width = grid_height = 0; +// colors.clear(); + confidences.clear(); + flags.clear(); flag_curr = 0; + normals = NULL; pdir1.clear(); pdir2.clear(); + curv1.clear(); curv2.clear(); dcurv.clear(); + cornerareas.clear(); pointareas.clear(); +// bbox.valid = bsphere.valid = false; + neighbors.clear(); adjacentfaces.clear(); across_edge.clear(); + texCoords = NULL; + } + + // + // Input and output + // +protected: + // static bool read_helper(const char *filename, TriMesh *mesh); +public: + // static TriMesh *read(const char *filename); + // static TriMesh *read(const ::std::string &filename); + // bool write(const char *filename); + // bool write(const ::std::string &filename); + + + // + // Useful queries + // + + // Is vertex v on the mesh boundary? + bool is_bdy(int v) + { + if (neighbors.empty()) need_neighbors(); + if (adjacentfaces.empty()) need_adjacentfaces(); + return neighbors[v].size() != adjacentfaces[v].size(); + } + + // Centroid of face f + osg::Vec3 centroid(int f) + { + if (faces.empty()) need_faces(); + return ( (*vertices)[faces[f][0]] + (*vertices)[faces[f][1]] + (*vertices)[faces[f][2]] ) * (1.0f / 3.0f); + } + + // Normal of face f + osg::Vec3 trinorm(int f) + { + if (faces.empty()) need_faces(); + return trinorm( (*vertices)[faces[f][0]], (*vertices)[faces[f][1]], (*vertices)[faces[f][2]]); + } + + // Angle of corner j in triangle + float cornerangle(int i, int j) + { + using namespace ::std; + + if (faces.empty()) need_faces(); + const osg::Vec3 &p0 = (*vertices)[faces[i][j]]; + const osg::Vec3 &p1 = (*vertices)[faces[i][(j+1)%3]]; + const osg::Vec3 &p2 = (*vertices)[faces[i][(j+2)%3]]; + + return acos( (p1 - p0) * (p2 - p0) ); + } + + // Dihedral angle between face i and face across_edge[i][j] + float dihedral(int i, int j) + { + if (across_edge.empty()) need_across_edge(); + if (across_edge[i][j] < 0) return 0.0f; + osg::Vec3 mynorm = trinorm(i); + osg::Vec3 othernorm = trinorm(across_edge[i][j]); + float ang = angle(mynorm, othernorm); + osg::Vec3 towards = ( ((*vertices)[faces[i][(j+1)%3]] + (*vertices)[faces[i][(j+2)%3]]) - (*vertices)[faces[i][j]] ) * 0.5f; + if ( (towards * othernorm) < 0.0f) + return SG_PI + ang; + else + return SG_PI - ang; + } + + // Statistics +// float stat(StatOp op, StatVal val); +// float feature_size(); + + // + // Debugging + // + + // Debugging printout, controllable by a "verbose"ness parameter +// static int verbose; +// static void set_verbose(int); +// static void (*dprintf_hook)(const char *); +// static void set_dprintf_hook(void (*hook)(const char *)); +// static void dprintf(const char *format, ...); + + // Same as above, but fatal-error printout +// static void (*eprintf_hook)(const char *); +// static void set_eprintf_hook(void (*hook)(const char *)); +// static void eprintf(const char *format, ...); +}; + +#endif diff --git a/simgear/scene/dem/TriMesh.h b/simgear/scene/dem/TriMesh.h new file mode 100644 index 00000000..00327bd1 --- /dev/null +++ b/simgear/scene/dem/TriMesh.h @@ -0,0 +1,257 @@ +#ifndef TRIMESH_H +#define TRIMESH_H + +#include + +/* +Szymon Rusinkiewicz +Princeton University + +TriMesh.h +Class for triangle meshes. +*/ + +#include "Vec.h" +#include "Box.h" +#include "Color.h" +#include +#include +#ifndef M_PIf +# define M_PIf 3.1415927f +#endif + +namespace trimesh { + +class TriMesh { +public: + // + // Types + // + struct Face { + int v[3]; + + Face() {} + Face(const int &v0, const int &v1, const int &v2) + { v[0] = v0; v[1] = v1; v[2] = v2; } + Face(const int *v_) + { v[0] = v_[0]; v[1] = v_[1]; v[2] = v_[2]; } + template explicit Face(const S &x) + { v[0] = x[0]; v[1] = x[1]; v[2] = x[2]; } + int &operator[] (int i) { return v[i]; } + const int &operator[] (int i) const { return v[i]; } + operator const int * () const { return &(v[0]); } + operator const int * () { return &(v[0]); } + operator int * () { return &(v[0]); } + int indexof(int v_) const + { + return (v[0] == v_) ? 0 : + (v[1] == v_) ? 1 : + (v[2] == v_) ? 2 : -1; + } + }; + + struct BSphere { + point center; + float r; + bool valid; + BSphere() : valid(false) + {} + }; + + // + // Enums + // + enum TstripRep { TSTRIP_LENGTH, TSTRIP_TERM }; + enum { GRID_INVALID = -1 }; + enum StatOp { STAT_MIN, STAT_MAX, STAT_MEAN, STAT_MEANABS, + STAT_RMS, STAT_MEDIAN, STAT_STDEV, STAT_TOTAL }; + enum StatVal { STAT_VALENCE, STAT_FACEAREA, STAT_ANGLE, + STAT_DIHEDRAL, STAT_EDGELEN, STAT_X, STAT_Y, STAT_Z }; + + // + // Constructor + // + TriMesh() : grid_width(-1), grid_height(-1), flag_curr(0) + {} + + // + // Members + // + + // The basics: vertices and faces + ::std::vector vertices; + ::std::vector faces; + + // Triangle strips + ::std::vector tstrips; + + // Grid, if present + ::std::vector grid; + int grid_width, grid_height; + + // Other per-vertex properties + ::std::vector colors; + ::std::vector confidences; + ::std::vector flags; + unsigned flag_curr; + + // Computed per-vertex properties + ::std::vector normals; + ::std::vector pdir1, pdir2; + ::std::vector curv1, curv2; + ::std::vector< Vec<4,float> > dcurv; + ::std::vector cornerareas; + ::std::vector pointareas; + + // Bounding structures + box bbox; + BSphere bsphere; + + // Connectivity structures: + // For each vertex, all neighboring vertices + ::std::vector< ::std::vector > neighbors; + // For each vertex, all neighboring faces + ::std::vector< ::std::vector > adjacentfaces; + // For each face, the three faces attached to its edges + // (for example, across_edge[3][2] is the number of the face + // that's touching the edge opposite vertex 2 of face 3) + ::std::vector across_edge; + + // + // Compute all this stuff... + // + void need_tstrips(); + void convert_strips(TstripRep rep); + void unpack_tstrips(); + void triangulate_grid(bool remove_slivers = true); + void need_faces() + { + if (!faces.empty()) + return; + if (!tstrips.empty()) + unpack_tstrips(); + else if (!grid.empty()) + triangulate_grid(); + } + void need_normals(); + void need_pointareas(); + void need_curvatures(); + void need_dcurv(); + void need_bbox(); + void need_bsphere(); + void need_neighbors(); + void need_adjacentfaces(); + void need_across_edge(); + + // + // Delete everything + // + void clear() + { + vertices.clear(); faces.clear(); tstrips.clear(); + grid.clear(); grid_width = grid_height = -1; + colors.clear(); confidences.clear(); + flags.clear(); flag_curr = 0; + normals.clear(); pdir1.clear(); pdir2.clear(); + curv1.clear(); curv2.clear(); dcurv.clear(); + cornerareas.clear(); pointareas.clear(); + bbox.valid = bsphere.valid = false; + neighbors.clear(); adjacentfaces.clear(); across_edge.clear(); + } + + // + // Input and output + // +protected: + static bool read_helper(const char *filename, TriMesh *mesh); +public: + static TriMesh *read(const char *filename); + static TriMesh *read(const ::std::string &filename); + bool write(const char *filename); + bool write(const ::std::string &filename); + + + // + // Useful queries + // + + // Is vertex v on the mesh boundary? + bool is_bdy(int v) + { + if (neighbors.empty()) need_neighbors(); + if (adjacentfaces.empty()) need_adjacentfaces(); + return neighbors[v].size() != adjacentfaces[v].size(); + } + + // Centroid of face f + vec centroid(int f) + { + if (faces.empty()) need_faces(); + return (1.0f / 3.0f) * + (vertices[faces[f][0]] + + vertices[faces[f][1]] + + vertices[faces[f][2]]); + } + + // Normal of face f + vec trinorm(int f) + { + if (faces.empty()) need_faces(); + return trimesh::trinorm(vertices[faces[f][0]], vertices[faces[f][1]], + vertices[faces[f][2]]); + } + + // Angle of corner j in triangle i + float cornerangle(int i, int j) + { + using namespace ::std; + + if (faces.empty()) need_faces(); + const point &p0 = vertices[faces[i][j]]; + const point &p1 = vertices[faces[i][(j+1)%3]]; + const point &p2 = vertices[faces[i][(j+2)%3]]; + return acos((p1 - p0) DOT (p2 - p0)); + } + + // Dihedral angle between face i and face across_edge[i][j] + float dihedral(int i, int j) + { + if (across_edge.empty()) need_across_edge(); + if (across_edge[i][j] < 0) return 0.0f; + vec mynorm = trinorm(i); + vec othernorm = trinorm(across_edge[i][j]); + float ang = angle(mynorm, othernorm); + vec towards = 0.5f * (vertices[faces[i][(j+1)%3]] + + vertices[faces[i][(j+2)%3]]) - + vertices[faces[i][j]]; + if ((towards DOT othernorm) < 0.0f) + return M_PIf + ang; + else + return M_PIf - ang; + } + + // Statistics + float stat(StatOp op, StatVal val); + float feature_size(); + + // + // Debugging + // + + // Debugging printout, controllable by a "verbose"ness parameter + static int verbose; + static void set_verbose(int); + static void (*dprintf_hook)(const char *); + static void set_dprintf_hook(void (*hook)(const char *)); + static void dprintf(const char *format, ...); + + // Same as above, but fatal-error printout + static void (*eprintf_hook)(const char *); + static void set_eprintf_hook(void (*hook)(const char *)); + static void eprintf(const char *format, ...); + +}; + +}; // namespace trimesh + +#endif diff --git a/simgear/scene/dem/TriMesh_curvature.cc b/simgear/scene/dem/TriMesh_curvature.cc new file mode 100644 index 00000000..cbf514bf --- /dev/null +++ b/simgear/scene/dem/TriMesh_curvature.cc @@ -0,0 +1,331 @@ +/* +Szymon Rusinkiewicz +Princeton University + +TriMesh_curvature.cc +Computation of per-vertex principal curvatures and directions. + +Uses algorithm from + Rusinkiewicz, Szymon. + "Estimating Curvatures and Their Derivatives on Triangle Meshes," + Proc. 3DPVT, 2004. +*/ + +#include "TriMesh.h" +#include "TriMesh_algo.h" +#include "lineqn.h" +using namespace std; + + +// i+1 and i-1 modulo 3 +// This way of computing it tends to be faster than using % +#define NEXT(i) ((i)<2 ? (i)+1 : (i)-2) +#define PREV(i) ((i)>0 ? (i)-1 : (i)+2) + + +namespace trimesh { + +// Rotate a coordinate system to be perpendicular to the given normal +static void rot_coord_sys(const vec &old_u, const vec &old_v, + const vec &new_norm, + vec &new_u, vec &new_v) +{ + new_u = old_u; + new_v = old_v; + vec old_norm = old_u CROSS old_v; + float ndot = old_norm DOT new_norm; + if (unlikely(ndot <= -1.0f)) { + new_u = -new_u; + new_v = -new_v; + return; + } + vec perp_old = new_norm - ndot * old_norm; + vec dperp = 1.0f / (1 + ndot) * (old_norm + new_norm); + new_u -= dperp * (new_u DOT perp_old); + new_v -= dperp * (new_v DOT perp_old); +} + + +// Reproject a curvature tensor from the basis spanned by old_u and old_v +// (which are assumed to be unit-length and perpendicular) to the +// new_u, new_v basis. +void proj_curv(const vec &old_u, const vec &old_v, + float old_ku, float old_kuv, float old_kv, + const vec &new_u, const vec &new_v, + float &new_ku, float &new_kuv, float &new_kv) +{ + vec r_new_u, r_new_v; + rot_coord_sys(new_u, new_v, old_u CROSS old_v, r_new_u, r_new_v); + + float u1 = r_new_u DOT old_u; + float v1 = r_new_u DOT old_v; + float u2 = r_new_v DOT old_u; + float v2 = r_new_v DOT old_v; + new_ku = old_ku * u1*u1 + old_kuv * (2.0f * u1*v1) + old_kv * v1*v1; + new_kuv = old_ku * u1*u2 + old_kuv * (u1*v2 + u2*v1) + old_kv * v1*v2; + new_kv = old_ku * u2*u2 + old_kuv * (2.0f * u2*v2) + old_kv * v2*v2; +} + + +// Like the above, but for dcurv +void proj_dcurv(const vec &old_u, const vec &old_v, + const Vec<4> old_dcurv, + const vec &new_u, const vec &new_v, + Vec<4> &new_dcurv) +{ + vec r_new_u, r_new_v; + rot_coord_sys(new_u, new_v, old_u CROSS old_v, r_new_u, r_new_v); + + float u1 = r_new_u DOT old_u; + float v1 = r_new_u DOT old_v; + float u2 = r_new_v DOT old_u; + float v2 = r_new_v DOT old_v; + + new_dcurv[0] = old_dcurv[0]*u1*u1*u1 + + old_dcurv[1]*3.0f*u1*u1*v1 + + old_dcurv[2]*3.0f*u1*v1*v1 + + old_dcurv[3]*v1*v1*v1; + new_dcurv[1] = old_dcurv[0]*u1*u1*u2 + + old_dcurv[1]*(u1*u1*v2 + 2.0f*u2*u1*v1) + + old_dcurv[2]*(u2*v1*v1 + 2.0f*u1*v1*v2) + + old_dcurv[3]*v1*v1*v2; + new_dcurv[2] = old_dcurv[0]*u1*u2*u2 + + old_dcurv[1]*(u2*u2*v1 + 2.0f*u1*u2*v2) + + old_dcurv[2]*(u1*v2*v2 + 2.0f*u2*v2*v1) + + old_dcurv[3]*v1*v2*v2; + new_dcurv[3] = old_dcurv[0]*u2*u2*u2 + + old_dcurv[1]*3.0f*u2*u2*v2 + + old_dcurv[2]*3.0f*u2*v2*v2 + + old_dcurv[3]*v2*v2*v2; +} + + +// Given a curvature tensor, find principal directions and curvatures +// Makes sure that pdir1 and pdir2 are perpendicular to normal +void diagonalize_curv(const vec &old_u, const vec &old_v, + float ku, float kuv, float kv, + const vec &new_norm, + vec &pdir1, vec &pdir2, float &k1, float &k2) +{ + vec r_old_u, r_old_v; + rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v); + + float c = 1, s = 0, tt = 0; + if (likely(kuv != 0.0f)) { + // Jacobi rotation to diagonalize + float h = 0.5f * (kv - ku) / kuv; + tt = (h < 0.0f) ? + 1.0f / (h - sqrt(1.0f + h*h)) : + 1.0f / (h + sqrt(1.0f + h*h)); + c = 1.0f / sqrt(1.0f + tt*tt); + s = tt * c; + } + + k1 = ku - tt * kuv; + k2 = kv + tt * kuv; + + if (fabs(k1) >= fabs(k2)) { + pdir1 = c*r_old_u - s*r_old_v; + } else { + swap(k1, k2); + pdir1 = s*r_old_u + c*r_old_v; + } + pdir2 = new_norm CROSS pdir1; +} + + +// Compute principal curvatures and directions. +void TriMesh::need_curvatures() +{ + if (curv1.size() == vertices.size()) + return; + need_faces(); + need_normals(); + need_pointareas(); + + dprintf("Computing curvatures... "); + + // Resize the arrays we'll be using + int nv = vertices.size(), nf = faces.size(); + curv1.clear(); curv1.resize(nv); curv2.clear(); curv2.resize(nv); + pdir1.clear(); pdir1.resize(nv); pdir2.clear(); pdir2.resize(nv); + vector curv12(nv); + + // Set up an initial coordinate system per vertex + for (int i = 0; i < nf; i++) { + pdir1[faces[i][0]] = vertices[faces[i][1]] - + vertices[faces[i][0]]; + pdir1[faces[i][1]] = vertices[faces[i][2]] - + vertices[faces[i][1]]; + pdir1[faces[i][2]] = vertices[faces[i][0]] - + vertices[faces[i][2]]; + } +#pragma omp parallel for + for (int i = 0; i < nv; i++) { + pdir1[i] = pdir1[i] CROSS normals[i]; + normalize(pdir1[i]); + pdir2[i] = normals[i] CROSS pdir1[i]; + } + + // Compute curvature per-face +#pragma omp parallel for + for (int i = 0; i < nf; i++) { + // Edges + vec e[3] = { vertices[faces[i][2]] - vertices[faces[i][1]], + vertices[faces[i][0]] - vertices[faces[i][2]], + vertices[faces[i][1]] - vertices[faces[i][0]] }; + + // N-T-B coordinate system per face + vec t = e[0]; + normalize(t); + vec n = e[0] CROSS e[1]; + vec b = n CROSS t; + normalize(b); + + // Estimate curvature based on variation of normals + // along edges + float m[3] = { 0, 0, 0 }; + float w[3][3] = { {0,0,0}, {0,0,0}, {0,0,0} }; + for (int j = 0; j < 3; j++) { + float u = e[j] DOT t; + float v = e[j] DOT b; + w[0][0] += u*u; + w[0][1] += u*v; + //w[1][1] += v*v + u*u; + //w[1][2] += u*v; + w[2][2] += v*v; + vec dn = normals[faces[i][PREV(j)]] - + normals[faces[i][NEXT(j)]]; + float dnu = dn DOT t; + float dnv = dn DOT b; + m[0] += dnu*u; + m[1] += dnu*v + dnv*u; + m[2] += dnv*v; + } + w[1][1] = w[0][0] + w[2][2]; + w[1][2] = w[0][1]; + + // Least squares solution + float diag[3]; + if (!ldltdc(w, diag)) { + //dprintf("ldltdc failed!\n"); + continue; + } + ldltsl(w, diag, m, m); + + // Push it back out to the vertices + for (int j = 0; j < 3; j++) { + int vj = faces[i][j]; + float c1, c12, c2; + proj_curv(t, b, m[0], m[1], m[2], + pdir1[vj], pdir2[vj], c1, c12, c2); + float wt = cornerareas[i][j] / pointareas[vj]; +#pragma omp atomic + curv1[vj] += wt * c1; +#pragma omp atomic + curv12[vj] += wt * c12; +#pragma omp atomic + curv2[vj] += wt * c2; + } + } + + // Compute principal directions and curvatures at each vertex +#pragma omp parallel for + for (int i = 0; i < nv; i++) { + diagonalize_curv(pdir1[i], pdir2[i], + curv1[i], curv12[i], curv2[i], + normals[i], pdir1[i], pdir2[i], + curv1[i], curv2[i]); + } + dprintf("Done.\n"); +} + + +// Compute derivatives of curvature. +void TriMesh::need_dcurv() +{ + if (dcurv.size() == vertices.size()) + return; + need_curvatures(); + + dprintf("Computing dcurv... "); + + // Resize the arrays we'll be using + int nv = vertices.size(), nf = faces.size(); + dcurv.clear(); dcurv.resize(nv); + + // Compute dcurv per-face +#pragma omp parallel for + for (int i = 0; i < nf; i++) { + // Edges + vec e[3] = { vertices[faces[i][2]] - vertices[faces[i][1]], + vertices[faces[i][0]] - vertices[faces[i][2]], + vertices[faces[i][1]] - vertices[faces[i][0]] }; + + // N-T-B coordinate system per face + vec t = e[0]; + normalize(t); + vec n = e[0] CROSS e[1]; + vec b = n CROSS t; + normalize(b); + + // Project curvature tensor from each vertex into this + // face's coordinate system + vec fcurv[3]; + for (int j = 0; j < 3; j++) { + int vj = faces[i][j]; + proj_curv(pdir1[vj], pdir2[vj], curv1[vj], 0, curv2[vj], + t, b, fcurv[j][0], fcurv[j][1], fcurv[j][2]); + + } + + // Estimate dcurv based on variation of curvature along edges + float m[4] = { 0, 0, 0, 0 }; + float w[4][4] = { {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} }; + for (int j = 0; j < 3; j++) { + // Variation of curvature along each edge + vec dfcurv = fcurv[PREV(j)] - fcurv[NEXT(j)]; + float u = e[j] DOT t; + float v = e[j] DOT b; + float u2 = u*u, v2 = v*v, uv = u*v; + w[0][0] += u2; + w[0][1] += uv; + //w[1][1] += 2.0f*u2 + v2; + //w[1][2] += 2.0f*uv; + //w[2][2] += u2 + 2.0f*v2; + //w[2][3] += uv; + w[3][3] += v2; + m[0] += u*dfcurv[0]; + m[1] += v*dfcurv[0] + 2.0f*u*dfcurv[1]; + m[2] += 2.0f*v*dfcurv[1] + u*dfcurv[2]; + m[3] += v*dfcurv[2]; + } + w[1][1] = 2.0f * w[0][0] + w[3][3]; + w[1][2] = 2.0f * w[0][1]; + w[2][2] = w[0][0] + 2.0f * w[3][3]; + w[2][3] = w[0][1]; + + // Least squares solution + float d[4]; + if (!ldltdc(w, d)) { + //dprintf("ldltdc failed!\n"); + continue; + } + ldltsl(w, d, m, m); + Vec<4> face_dcurv(m); + + // Push it back out to each vertex + for (int j = 0; j < 3; j++) { + int vj = faces[i][j]; + Vec<4> this_vert_dcurv; + proj_dcurv(t, b, face_dcurv, + pdir1[vj], pdir2[vj], this_vert_dcurv); + float wt = cornerareas[i][j] / pointareas[vj]; + dcurv[vj] += wt * this_vert_dcurv; + } + } + + dprintf("Done.\n"); +} + +}; // namespace trimesh diff --git a/simgear/scene/dem/TriMesh_normals.cc b/simgear/scene/dem/TriMesh_normals.cc new file mode 100644 index 00000000..0458ebd8 --- /dev/null +++ b/simgear/scene/dem/TriMesh_normals.cc @@ -0,0 +1,116 @@ +/* +Szymon Rusinkiewicz +Princeton University + +TriMesh_normals.cc +Compute per-vertex normals for TriMeshes + +For meshes, uses average of per-face normals, weighted according to: + Max, N. + "Weights for Computing Vertex Normals from Facet Normals," + Journal of Graphics Tools, Vol. 4, No. 2, 1999. + +For raw point clouds, fits plane to k nearest neighbors. +*/ + +#include "TriMesh.h" +#include "KDtree.h" +#include "lineqn.h" +using namespace std; + + +namespace trimesh { + +// Compute per-vertex normals +void TriMesh::need_normals() +{ + // Nothing to do if we already have normals + int nv = vertices.size(); + if (int(normals.size()) == nv) + return; + + dprintf("Computing normals... "); + normals.clear(); + normals.resize(nv); + + // TODO: direct handling of grids + if (!tstrips.empty()) { + // Compute from tstrips + const int *t = &tstrips[0], *end = t + tstrips.size(); + while (likely(t < end)) { + int striplen = *t - 2; + t += 3; + bool flip = false; + for (int i = 0; i < striplen; i++, t++, flip = !flip) { + const point &p0 = vertices[*(t-2)]; + const point &p1 = vertices[*(t-1)]; + const point &p2 = vertices[* t ]; + vec a = p0-p1, b = p1-p2, c = p2-p0; + float l2a = len2(a), l2b = len2(b), l2c = len2(c); + if (!l2a || !l2b || !l2c) + continue; + vec facenormal = flip ? (b CROSS a) : (a CROSS b); + normals[*(t-2)] += facenormal * (1.0f / (l2a * l2c)); + normals[*(t-1)] += facenormal * (1.0f / (l2b * l2a)); + normals[* t ] += facenormal * (1.0f / (l2c * l2b)); + } + } + } else if (need_faces(), !faces.empty()) { + // Compute from faces + int nf = faces.size(); +#pragma omp parallel for + for (int i = 0; i < nf; i++) { + const point &p0 = vertices[faces[i][0]]; + const point &p1 = vertices[faces[i][1]]; + const point &p2 = vertices[faces[i][2]]; + vec a = p0-p1, b = p1-p2, c = p2-p0; + float l2a = len2(a), l2b = len2(b), l2c = len2(c); + if (!l2a || !l2b || !l2c) + continue; + vec facenormal = a CROSS b; + normals[faces[i][0]] += facenormal * (1.0f / (l2a * l2c)); + normals[faces[i][1]] += facenormal * (1.0f / (l2b * l2a)); + normals[faces[i][2]] += facenormal * (1.0f / (l2c * l2b)); + } + } else { + // Find normals of a point cloud + const int k = 6; + const vec ref(0, 0, 1); + KDtree kd(vertices); +#pragma omp parallel for + for (int i = 0; i < nv; i++) { + vector knn; + kd.find_k_closest_to_pt(knn, k, vertices[i]); + int actual_k = knn.size(); + if (actual_k < 3) { + dprintf("Warning: not enough points for vertex %d\n", i); + normals[i] = ref; + continue; + } + // Compute covariance + float C[3][3] = { {0,0,0}, {0,0,0}, {0,0,0} }; + // The below loop starts at 1, since element 0 + // is just vertices[i] itself + for (int j = 1; j < actual_k; j++) { + vec d = point(knn[j]) - vertices[i]; + for (int l = 0; l < 3; l++) + for (int m = 0; m < 3; m++) + C[l][m] += d[l] * d[m]; + } + float e[3]; + eigdc(C, e); + normals[i] = vec(C[0][0], C[1][0], C[2][0]); + if ((normals[i] DOT ref) < 0.0f) + normals[i] = -normals[i]; + } + } + + // Make them all unit-length +#pragma omp parallel for + for (int i = 0; i < nv; i++) + normalize(normals[i]); + + dprintf("Done.\n"); +} + +}; // namespace trimesh diff --git a/simgear/scene/tgdb/BucketBox.hxx b/simgear/scene/tgdb/BucketBox.hxx index 2e259f89..8cc616ac 100644 --- a/simgear/scene/tgdb/BucketBox.hxx +++ b/simgear/scene/tgdb/BucketBox.hxx @@ -26,8 +26,16 @@ #include #include +#include +#include + #include #include +#include + +#ifdef ENABLE_GDAL +#include +#endif namespace simgear { @@ -293,7 +301,7 @@ public: } return numTiles; } - + unsigned getTileTriangles(unsigned i, unsigned j, unsigned width, unsigned height, SGVec3f points[6], SGVec3f normals[6], SGVec2f texCoords[6]) const { @@ -304,7 +312,7 @@ public: unsigned y0 = _offset[1] + j; unsigned y1 = y0 + height; - + SGGeod p00 = _offsetToGeod(x0, y0, 0); SGVec3f v00 = SGVec3f::fromGeod(p00); SGVec3f n00 = SGQuatf::fromLonLat(p00).backTransform(SGVec3f(0, 0, -1)); @@ -324,7 +332,7 @@ public: SGVec3f v01 = SGVec3f::fromGeod(p01); SGVec3f n01 = SGQuatf::fromLonLat(p01).backTransform(SGVec3f(0, 0, -1)); SGVec2f t01(x0*1.0/(360*8), y1*1.0/(180*8)); - + if (y0 != 0) { points[numPoints] = v00; normals[numPoints] = n00; @@ -360,6 +368,33 @@ public: return numPoints; } +#ifdef ENABLE_GDAL + osg::Geode* getTileTriangleMesh( const SGDemPtr dem, unsigned res, SGMesh::TextureMethod tm, const osgDB::Options* options ) const + { + unsigned widthLevel = getWidthLevel(); + unsigned heightLevel = getHeightLevel(); + + unsigned x0 = _offset[0]; + unsigned x1 = x0 + _size[0]; + + unsigned y0 = _offset[1]; + unsigned y1 = y0 + _size[1]; + + SGSpheref sphere = getBoundingSphere(); + osg::Matrixd transform; + transform.makeTranslate(toOsg(-sphere.getCenter())); + + // create a mesh of this dimension + if ( (y0 != 0) && (y1 != 180*8) ) { + SGMesh mesh( dem, x0, y0, x1, y1, heightLevel, widthLevel, transform, tm, options ); + return mesh.getGeode(); + } else { + // todo - handle poles + return 0; + } + } +#endif + private: static unsigned _normalizeLongitude(unsigned offset) { return offset - (360*8)*(offset/(360*8)); } @@ -398,7 +433,7 @@ private: static SGGeod _offsetToGeod(unsigned offset0, unsigned offset1, double elev) { return SGGeod::fromDegM(_offsetToLongitudeDeg(offset0), _offsetToLatitudeDeg(offset1), elev); } - + static unsigned _getLevel(const unsigned factors[], unsigned nFactors, unsigned begin, unsigned end) { unsigned rbegin = end - 1; diff --git a/simgear/scene/tgdb/ReaderWriterSPT.cxx b/simgear/scene/tgdb/ReaderWriterSPT.cxx index 3116b2c4..cb6f4df1 100644 --- a/simgear/scene/tgdb/ReaderWriterSPT.cxx +++ b/simgear/scene/tgdb/ReaderWriterSPT.cxx @@ -6,7 +6,7 @@ // modify it under the terms of the GNU General Public License as // published by the Free Software Foundation; either version 2 of the // License, or (at your option) any later version. -// +// // This program is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU @@ -112,7 +112,7 @@ struct ReaderWriterSPT::CullCallback : public osg::NodeCallback { }; struct ReaderWriterSPT::LocalOptions { - LocalOptions(const osgDB::Options* options) : + LocalOptions(const osgDB::Options* options) : _options(options) { std::string pageLevelsString; @@ -213,7 +213,7 @@ ReaderWriterSPT::readObject(const std::string& fileName, const osgDB::Options* o texture->setImage(image); texture->setWrap(osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT); texture->setWrap(osg::Texture2D::WRAP_T, osg::Texture2D::CLAMP); - stateSet->setTextureAttributeAndModes(0, texture); + stateSet->setTextureAttributeAndModes(0, texture); } return stateSet; @@ -234,7 +234,7 @@ ReaderWriterSPT::readNode(const std::string& fileName, const osgDB::Options* opt ss >> bucketBox; if (ss.fail()) return ReadResult::FILE_NOT_FOUND; - + BucketBox bucketBoxList[2]; unsigned bucketBoxListSize = bucketBox.periodicSplit(bucketBoxList); if (bucketBoxListSize == 0) @@ -265,7 +265,7 @@ ReaderWriterSPT::createTree(const BucketBox& bucketBox, const LocalOptions& opti if (numTiles == 0) return 0; - if (numTiles == 1) + if (numTiles == 1) return createTree(bucketBoxList[0], options, false); osg::ref_ptr group = new osg::Group; @@ -277,7 +277,7 @@ ReaderWriterSPT::createTree(const BucketBox& bucketBox, const LocalOptions& opti } if (!group->getNumChildren()) return 0; - + return group; } } @@ -332,7 +332,7 @@ ReaderWriterSPT::createPagedLOD(const BucketBox& bucketBox, const LocalOptions& ss << bucketBox << ".spt"; pagedLOD->setFileName(pagedLOD->getNumChildren(), ss.str()); pagedLOD->setRange(pagedLOD->getNumChildren(), 0.0, range); - + return pagedLOD; } @@ -349,7 +349,7 @@ ReaderWriterSPT::createSeaLevelTile(const BucketBox& bucketBox, const LocalOptio osg::Vec3Array* vertices = new osg::Vec3Array; osg::Vec3Array* normals = new osg::Vec3Array; osg::Vec2Array* texCoords = new osg::Vec2Array; - + unsigned widthLevel = bucketBox.getWidthLevel(); unsigned heightLevel = bucketBox.getHeightLevel(); @@ -373,10 +373,10 @@ ReaderWriterSPT::createSeaLevelTile(const BucketBox& bucketBox, const LocalOptio i += incx; incx = std::min(incx, bucketBox.getSize(0) - i); } - + osg::Vec4Array* colors = new osg::Vec4Array; colors->push_back(osg::Vec4(1, 1, 1, 1)); - + osg::Geometry* geometry = new osg::Geometry; geometry->setDataVariance(osg::Object::STATIC); geometry->setUseVertexBufferObjects(true); @@ -386,11 +386,11 @@ ReaderWriterSPT::createSeaLevelTile(const BucketBox& bucketBox, const LocalOptio geometry->setColorArray(colors); geometry->setColorBinding(osg::Geometry::BIND_OVERALL); geometry->setTexCoordArray(0, texCoords); - + osg::DrawArrays* drawArrays = new osg::DrawArrays(osg::DrawArrays::TRIANGLES, 0, vertices->size()); drawArrays->setDataVariance(osg::Object::STATIC); geometry->addPrimitiveSet(drawArrays); - + osg::Geode* geode = new osg::Geode; geode->setDataVariance(osg::Object::STATIC); geode->addDrawable(geometry); diff --git a/simgear/scene/tgdb/userdata.cxx b/simgear/scene/tgdb/userdata.cxx index 2482f033..3f81bd81 100644 --- a/simgear/scene/tgdb/userdata.cxx +++ b/simgear/scene/tgdb/userdata.cxx @@ -38,6 +38,7 @@ #include "SGReaderWriterBTG.hxx" #include "ReaderWriterSPT.hxx" #include "ReaderWriterSTG.hxx" +#include // the following are static values needed by the runtime object // loader. However, the loading is done via a call back so these @@ -62,6 +63,9 @@ simgear::ModelRegistryCallbackProxy g_stgCallbackProx osgDB::RegisterReaderWriterProxy g_readerWriterSPTProxy; simgear::ModelRegistryCallbackProxy g_sptCallbackProxy("spt"); + +osgDB::RegisterReaderWriterProxy g_readerWriterPGTProxy; +simgear::ModelRegistryCallbackProxy g_pgtCallbackProxy("pgt"); } void sgUserDataInit( SGPropertyNode *p ) { diff --git a/simgear/scene/util/SGReaderWriterOptions.hxx b/simgear/scene/util/SGReaderWriterOptions.hxx index a67a6cb4..ce48dce9 100644 --- a/simgear/scene/util/SGReaderWriterOptions.hxx +++ b/simgear/scene/util/SGReaderWriterOptions.hxx @@ -26,6 +26,10 @@ #include +#ifdef ENABLE_GDAL +#include +#endif + class SGPropertyNode; typedef std::vector < std::string > string_list; @@ -60,6 +64,9 @@ public: osgDB::Options(options, copyop), _propertyNode(options._propertyNode), _materialLib(options._materialLib), +#ifdef ENABLE_GDAL + _dem(options._dem), +#endif _load_panel(options._load_panel), _model_data(options._model_data), _instantiateEffects(options._instantiateEffects), @@ -78,6 +85,13 @@ public: void setMaterialLib(SGMaterialLib* materialLib) { _materialLib = materialLib; } +#ifdef ENABLE_GDAL + SGDemPtr getDem() const + { return _dem; } + void setDem(SGDem* dem) + { _dem = dem; } +#endif + typedef osg::Node *(*panel_func)(SGPropertyNode *); panel_func getLoadPanel() const @@ -110,6 +124,10 @@ protected: private: SGSharedPtr _propertyNode; SGSharedPtr _materialLib; +#ifdef ENABLE_GDAL + SGSharedPtr _dem; +#endif + osg::Node *(*_load_panel)(SGPropertyNode *); osg::ref_ptr _model_data; bool _instantiateEffects; diff --git a/simgear/simgear_config_cmake.h.in b/simgear/simgear_config_cmake.h.in index 461fa243..ba28e761 100644 --- a/simgear/simgear_config_cmake.h.in +++ b/simgear/simgear_config_cmake.h.in @@ -1,7 +1,7 @@ #cmakedefine HAVE_SYS_TIME_H #cmakedefine HAVE_SYS_TIMEB_H -#cmakedefine HAVE_UNISTD_H +#cmakedefine HAVE_UNISTD_H 1 #cmakedefine HAVE_GETTIMEOFDAY @@ -21,3 +21,4 @@ #cmakedefine SYSTEM_EXPAT #cmakedefine ENABLE_SOUND #cmakedefine ENABLE_SIMD +#cmakedefine ENABLE_GDAL