alternative terrain engine - SGMesh utilizing pagedLOD
This commit is contained in:
Peter Sadrozinski
@@ -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)
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -15,3 +15,6 @@ foreach( mylibfolder
|
||||
|
||||
endforeach( mylibfolder )
|
||||
|
||||
if(ENABLE_GDAL)
|
||||
add_subdirectory(dem)
|
||||
endif(ENABLE_GDAL)
|
||||
|
||||
24
simgear/scene/dem/CMakeLists.txt
Normal file
24
simgear/scene/dem/CMakeLists.txt
Normal file
@@ -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}")
|
||||
418
simgear/scene/dem/ReaderWriterPGT.cxx
Normal file
418
simgear/scene/dem/ReaderWriterPGT.cxx
Normal file
@@ -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 <simgear_config.h>
|
||||
#endif
|
||||
|
||||
#include <simgear/scene/util/SGReaderWriterOptions.hxx>
|
||||
#include "ReaderWriterPGT.hxx"
|
||||
|
||||
#include <cassert>
|
||||
|
||||
#include <osg/PagedLOD>
|
||||
#include <osg/MatrixTransform>
|
||||
#include <osg/Texture2D>
|
||||
|
||||
#include <osgDB/FileNameUtils>
|
||||
#include <osgDB/FileUtils>
|
||||
#include <osgDB/ReadFile>
|
||||
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#include <simgear/scene/dem/SGMesh.hxx>
|
||||
#include <simgear/scene/util/OsgMath.hxx>
|
||||
|
||||
#include <simgear/scene/tgdb/BucketBox.hxx>
|
||||
#include <simgear/scene/model/ModelRegistry.hxx>
|
||||
|
||||
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<SGReaderWriterOptions> 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<unsigned>::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<unsigned> _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<osg::Group> group = new osg::Group;
|
||||
group->addChild(createTree(bucketBoxList[0], localOptions, true));
|
||||
group->addChild(createTree(bucketBoxList[1], localOptions, true));
|
||||
return ReadResult(group);
|
||||
}
|
||||
|
||||
osg::ref_ptr<osg::Node>
|
||||
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<osg::Group> group = new osg::Group;
|
||||
for (unsigned i = 0; i < numTiles; ++i) {
|
||||
osg::ref_ptr<osg::Node> node = createTree(bucketBoxList[i], options, false);
|
||||
if (!node.valid())
|
||||
continue;
|
||||
group->addChild(node.get());
|
||||
}
|
||||
if (!group->getNumChildren())
|
||||
return 0;
|
||||
|
||||
return group;
|
||||
}
|
||||
}
|
||||
|
||||
osg::ref_ptr<osg::Node>
|
||||
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<osgDB::Options> localOptions;
|
||||
localOptions = static_cast<osgDB::Options*>(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<osg::Node> node = osgDB::readRefNodeFile(fileName, options._options);
|
||||
if (!node.valid())
|
||||
continue;
|
||||
pagedLOD->addChild(node.get(), range, std::numeric_limits<float>::max());
|
||||
break;
|
||||
}
|
||||
// Add the static sea level textured shell if there is nothing found
|
||||
if (pagedLOD->getNumChildren() == 0) {
|
||||
osg::ref_ptr<osg::Node> node = createTileMesh(bucketBox, options._options);
|
||||
if (node.valid())
|
||||
pagedLOD->addChild(node.get(), range, std::numeric_limits<float>::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<osg::Node>
|
||||
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<osg::StateSet>
|
||||
ReaderWriterPGT::getLowLODStateSet(const LocalOptions& options) const
|
||||
{
|
||||
osg::ref_ptr<osgDB::Options> localOptions;
|
||||
localOptions = static_cast<osgDB::Options*>(options._options->clone(osg::CopyOp()));
|
||||
localOptions->setObjectCacheHint(osgDB::Options::CACHE_ALL);
|
||||
|
||||
osg::ref_ptr<osg::Object> object = osgDB::readRefObjectFile("state.pgt", localOptions.get());
|
||||
if (!dynamic_cast<osg::StateSet*>(object.get()))
|
||||
return 0;
|
||||
|
||||
return static_cast<osg::StateSet*>(object.get());
|
||||
}
|
||||
|
||||
} // namespace simgear
|
||||
|
||||
// simgear::ModelRegistryCallbackProxy<simgear::LoadOnlyCallback> g_pgtCallbackProxy("pgt");
|
||||
53
simgear/scene/dem/ReaderWriterPGT.hxx
Normal file
53
simgear/scene/dem/ReaderWriterPGT.hxx
Normal file
@@ -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 <osgDB/ReaderWriter>
|
||||
|
||||
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<osg::Node> createTree(const BucketBox& bucketBox, const LocalOptions& options, bool topLevel) const;
|
||||
osg::ref_ptr<osg::Node> createPagedLOD(const BucketBox& bucketBox, const LocalOptions& options) const;
|
||||
osg::ref_ptr<osg::Node> createTileMesh(const BucketBox& bucketBox, const LocalOptions& options) const;
|
||||
osg::ref_ptr<osg::StateSet> getLowLODStateSet(const LocalOptions& options) const;
|
||||
|
||||
private:
|
||||
struct CullCallback;
|
||||
};
|
||||
|
||||
} // namespace simgear
|
||||
|
||||
#endif
|
||||
332
simgear/scene/dem/SGDem.cxx
Normal file
332
simgear/scene/dem/SGDem.cxx
Normal file
@@ -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 <fstream>
|
||||
#include <limits>
|
||||
|
||||
#include <boost/foreach.hpp>
|
||||
|
||||
#include <cpl_conv.h> // for CPLMalloc()
|
||||
#include "ogr_spatialref.h"
|
||||
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <simgear/misc/stdint.hxx>
|
||||
#include <simgear/misc/sg_dir.hxx>
|
||||
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#include <simgear/scene/dem/SGDemSession.hxx>
|
||||
|
||||
using namespace simgear;
|
||||
|
||||
#define DEM_DEBUG (0)
|
||||
|
||||
// where to add these...
|
||||
int SGDem::floorWithEpsilon(double x)
|
||||
{
|
||||
return static_cast<int>(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<demRoots.size(); i++ )
|
||||
{
|
||||
if ( demRoots[i].isValid( lvl, wo, so, eo, no ) ) {
|
||||
dr = &demRoots[i];
|
||||
break;
|
||||
} else {
|
||||
SG_LOG( SG_TERRAIN, SG_ALERT, "SGDem::findDem - dem " << i << " is not vald." );
|
||||
}
|
||||
}
|
||||
|
||||
return dr;
|
||||
}
|
||||
|
||||
#if 0 // todo - move to session
|
||||
unsigned short SGDem::getAlt( const SGDemSession& s, const SGGeod& loc ) const
|
||||
{
|
||||
int alt = 0;
|
||||
|
||||
// which level index are we querying?
|
||||
int lvlIndex = s.getLvlIndex();
|
||||
|
||||
if ( lvlIndex >= 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
|
||||
82
simgear/scene/dem/SGDem.hxx
Normal file
82
simgear/scene/dem/SGDem.hxx
Normal file
@@ -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 <map>
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/misc/sg_path.hxx>
|
||||
#include <simgear/scene/dem/SGDemRoot.hxx>
|
||||
|
||||
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<SGDemRoot> demRoots;
|
||||
};
|
||||
|
||||
typedef SGSharedPtr<SGDem> SGDemPtr;
|
||||
|
||||
#endif /* __SG_DEM_HXX__ */
|
||||
72
simgear/scene/dem/SGDemLevel.cxx
Normal file
72
simgear/scene/dem/SGDemLevel.cxx
Normal file
@@ -0,0 +1,72 @@
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
|
||||
#include <simgear/misc/sg_dir.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <simgear/scene/dem/SGDemLevel.hxx>
|
||||
|
||||
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" );
|
||||
}
|
||||
90
simgear/scene/dem/SGDemLevel.hxx
Normal file
90
simgear/scene/dem/SGDemLevel.hxx
Normal file
@@ -0,0 +1,90 @@
|
||||
#ifndef __SG_DEM_LEVEL_HXX__
|
||||
#define __SG_DEM_LEVEL_HXX__
|
||||
|
||||
#include <cstring>
|
||||
#include <cstdio>
|
||||
|
||||
#include <simgear/misc/sg_dir.hxx>
|
||||
|
||||
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__ */
|
||||
118
simgear/scene/dem/SGDemRoot.cxx
Normal file
118
simgear/scene/dem/SGDemRoot.cxx
Normal file
@@ -0,0 +1,118 @@
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
#include <sstream>
|
||||
|
||||
#include <simgear/misc/sg_dir.hxx>
|
||||
#include <simgear/scene/dem/SGDemRoot.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
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;
|
||||
}
|
||||
102
simgear/scene/dem/SGDemRoot.hxx
Normal file
102
simgear/scene/dem/SGDemRoot.hxx
Normal file
@@ -0,0 +1,102 @@
|
||||
#ifndef __SG_DEM_ROOT_HXX__
|
||||
#define __SG_DEM_ROOT_HXX__
|
||||
|
||||
#include <simgear/scene/dem/SGDemLevel.hxx>
|
||||
#include <simgear/scene/dem/SGDemTile.hxx>
|
||||
|
||||
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<SGDemLevel> levels;
|
||||
std::vector<SGDemCache> caches;
|
||||
};
|
||||
|
||||
#endif /* __SG_DEM_ROOT_HXX__ */
|
||||
79
simgear/scene/dem/SGDemSession.cxx
Normal file
79
simgear/scene/dem/SGDemSession.cxx
Normal file
@@ -0,0 +1,79 @@
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#include <simgear/scene/dem/SGDemSession.hxx>
|
||||
|
||||
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<SGGeod>& 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 ) );
|
||||
}
|
||||
}
|
||||
}
|
||||
78
simgear/scene/dem/SGDemSession.hxx
Normal file
78
simgear/scene/dem/SGDemSession.hxx
Normal file
@@ -0,0 +1,78 @@
|
||||
#ifndef __SG_DEM_SESSION_HXX__
|
||||
#define __SG_DEM_SESSION_HXX__
|
||||
|
||||
#include <simgear/scene/dem/SGDemRoot.hxx>
|
||||
|
||||
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<SGDemTileRef>& 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<SGGeod>& 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<SGDemTileRef> tileRefs;
|
||||
};
|
||||
|
||||
#endif /* __SG_DEM_SESSION_HXX__ */
|
||||
361
simgear/scene/dem/SGDemTile.cxx
Normal file
361
simgear/scene/dem/SGDemTile.cxx
Normal file
@@ -0,0 +1,361 @@
|
||||
#include <iomanip>
|
||||
#include <sstream>
|
||||
|
||||
#include <simgear/misc/sg_dir.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#include <simgear/scene/dem/SGDemTile.hxx>
|
||||
#include <simgear/scene/dem/SGDemSession.hxx>
|
||||
|
||||
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<SGDemTileRef>& tiles = s.getTiles();
|
||||
SG_LOG( SG_TERRAIN, SG_INFO, " create tile from " << tiles.size() << " source tiles" );
|
||||
|
||||
for ( unsigned int i=0; i<tiles.size(); i++ ) {
|
||||
// check if the file exists
|
||||
if ( tiles[i]->getPath().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<SGGeod>& 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();
|
||||
}
|
||||
65
simgear/scene/dem/SGDemTile.hxx
Normal file
65
simgear/scene/dem/SGDemTile.hxx
Normal file
@@ -0,0 +1,65 @@
|
||||
#ifndef __SG_DEM_TILE_HXX__
|
||||
#define __SG_DEM_TILE_HXX__
|
||||
|
||||
#include <gdal.h>
|
||||
#include <gdal_priv.h>
|
||||
|
||||
#include <map>
|
||||
|
||||
#include <simgear/structure/SGSharedPtr.hxx>
|
||||
#include <simgear/math/SGGeod.hxx>
|
||||
|
||||
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<SGGeod>& 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<SGDemTile> SGDemTileRef;
|
||||
typedef std::map<unsigned long, SGDemTileRef> SGDemCache;
|
||||
|
||||
#endif /* #define __SG_DEM_TILE_HXX__ */
|
||||
628
simgear/scene/dem/SGDemTile_gdal.cxx
Normal file
628
simgear/scene/dem/SGDemTile_gdal.cxx
Normal file
@@ -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 <warmerdam@pobox.com>
|
||||
*
|
||||
******************************************************************************
|
||||
* Copyright (c) 2002, i3 - information integration and imaging
|
||||
* Fort Collin, CO
|
||||
* Copyright (c) 2007-2013, Even Rouault <even dot rouault at mines-paris dot org>
|
||||
*
|
||||
* 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 <iostream>
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
|
||||
#include <boost/foreach.hpp>
|
||||
|
||||
#include <cpl_conv.h> // for CPLMalloc()
|
||||
#include "ogr_spatialref.h"
|
||||
|
||||
#include <gdal.h>
|
||||
#include <gdalwarper.h>
|
||||
#include <gdal_priv.h>
|
||||
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <simgear/misc/stdint.hxx>
|
||||
#include <simgear/misc/sg_dir.hxx>
|
||||
|
||||
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<GDALColorInterp> 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 );
|
||||
}
|
||||
414
simgear/scene/dem/SGMesh.cxx
Normal file
414
simgear/scene/dem/SGMesh.cxx
Normal file
@@ -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 <osg/Material>
|
||||
#include <osg/ShadeModel>
|
||||
#include <osg/CullFace>
|
||||
#include <osg/Texture2D>
|
||||
|
||||
#include <osgDB/FileNameUtils>
|
||||
#include <osgDB/FileUtils>
|
||||
#include <osgDB/ReadFile>
|
||||
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#include <simgear/scene/util/OsgMath.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#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<SGGeod> skirt_geods;
|
||||
::std::vector<unsigned int> 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<SGGeod> 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; i<skirt_geods.size(); i++ ) {
|
||||
skirt_geods[i].setElevationM( skirt_geods[i].getElevationM() - lower );
|
||||
}
|
||||
|
||||
for ( unsigned int i=0; i<index.size(); i++ ) {
|
||||
(*vertices)[index[i]].set( toOsg( SGVec3f::fromGeod( skirt_geods[i] ) ) );
|
||||
}
|
||||
|
||||
// translate pos after normals computed
|
||||
#pragma omp parallel for
|
||||
for ( unsigned int i=0; i < nv; i++ ) {
|
||||
(*vertices)[i].set( transform.preMult( (*vertices)[i]) );
|
||||
}
|
||||
|
||||
float startu = SGMiscd::normalizePeriodic( 0.0, 1.0, wo*1.0/(360*8) - incu);
|
||||
float startv = SGMiscd::normalizePeriodic( 0.0, 1.0, so*1.0/(180*8) - incv);
|
||||
|
||||
for ( unsigned int di = 0; di < grid_width; di++ ) {
|
||||
for ( unsigned int dj = 0; dj < grid_height; dj++ ) {
|
||||
(*texCoords)[di*grid_height+dj].set( toOsg(SGVec2f(startu + di*incu, startv + dj*incv)) );
|
||||
}
|
||||
}
|
||||
|
||||
osg::Vec4Array* colors = new osg::Vec4Array;
|
||||
if ( tm == SGMesh::TEXTURE_DEBUG ) {
|
||||
osg::Vec4 lvlColor;
|
||||
switch( heightLevel ) {
|
||||
case 9: lvlColor = osg::Vec4(0.5, 0.7, 0.9, 1); break;
|
||||
case 8: lvlColor = osg::Vec4(1, 1, 0, 1); break;
|
||||
case 7: lvlColor = osg::Vec4(0, 1, 0, 1); break;
|
||||
case 6: lvlColor = osg::Vec4(0, 1, 1, 1); break;
|
||||
|
||||
case 5: lvlColor = osg::Vec4(0, 0, 1, 1); break;
|
||||
case 4: lvlColor = osg::Vec4(1, 1, 1, 1); break;
|
||||
case 3: lvlColor = osg::Vec4(1, 0, 1, 1); break;
|
||||
|
||||
case 2: lvlColor = osg::Vec4(0.5, 0.5, 0.5, 1); break;
|
||||
case 1: lvlColor = osg::Vec4(0.5, 0, 0.5, 1); break;
|
||||
}
|
||||
|
||||
// special colors - green, and ref
|
||||
if ( Debug1 ) {
|
||||
lvlColor = osg::Vec4(0.0, 0.7, 0.0, 1);
|
||||
} else if ( Debug2 ) {
|
||||
lvlColor = osg::Vec4(0.7, 0.0, 0.0, 1);
|
||||
}
|
||||
|
||||
for ( unsigned int v=0; v<getVertices()->size(); 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;
|
||||
}
|
||||
718
simgear/scene/dem/SGMesh.hxx
Normal file
718
simgear/scene/dem/SGMesh.hxx
Normal file
@@ -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 <osg/Geometry>
|
||||
#include <osg/Geode>
|
||||
|
||||
#include <simgear/scene/util/SGReaderWriterOptions.hxx>
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#include <simgear/math/SGMath.hxx>
|
||||
|
||||
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 <class T> KDtree(const ::std::vector<T> &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<const float *> &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 <cmath>
|
||||
#include <algorithm>
|
||||
#include <limits>
|
||||
|
||||
// LU decomposition
|
||||
template <class T, int N>
|
||||
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 <class T, int N>
|
||||
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 <class T, int N>
|
||||
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 <class T, int N>
|
||||
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 <class T, int N>
|
||||
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<T>::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 <class T, int N>
|
||||
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 <class S> 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<Face> faces;
|
||||
|
||||
// // Triangle strips
|
||||
// ::std::vector<int> tstrips;
|
||||
|
||||
// Grid, if present
|
||||
// ::std::vector<int> grid;
|
||||
unsigned int grid_width, grid_height;
|
||||
|
||||
// Other per-vertex properties
|
||||
// ::std::vector<Color> colors;
|
||||
::std::vector<float> confidences;
|
||||
::std::vector<unsigned> flags;
|
||||
unsigned flag_curr;
|
||||
|
||||
// Computed per-vertex properties
|
||||
::std::vector<osg::Vec3> pdir1, pdir2;
|
||||
::std::vector<float> curv1, curv2;
|
||||
::std::vector<osg::Vec4> dcurv;
|
||||
::std::vector<osg::Vec3> cornerareas;
|
||||
::std::vector<float> pointareas;
|
||||
|
||||
osg::Geode* geode;
|
||||
|
||||
// Bounding structures
|
||||
//box bbox;
|
||||
//BSphere bsphere;
|
||||
|
||||
// Connectivity structures:
|
||||
// For each vertex, all neighboring vertices
|
||||
::std::vector< ::std::vector<unsigned int> > neighbors;
|
||||
// For each vertex, all neighboring faces
|
||||
::std::vector< ::std::vector<unsigned int> > 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<Face> 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
|
||||
257
simgear/scene/dem/TriMesh.h
Normal file
257
simgear/scene/dem/TriMesh.h
Normal file
@@ -0,0 +1,257 @@
|
||||
#ifndef TRIMESH_H
|
||||
#define TRIMESH_H
|
||||
|
||||
#include <simgear_config.h>
|
||||
|
||||
/*
|
||||
Szymon Rusinkiewicz
|
||||
Princeton University
|
||||
|
||||
TriMesh.h
|
||||
Class for triangle meshes.
|
||||
*/
|
||||
|
||||
#include "Vec.h"
|
||||
#include "Box.h"
|
||||
#include "Color.h"
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#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 <class S> 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<point> vertices;
|
||||
::std::vector<Face> faces;
|
||||
|
||||
// Triangle strips
|
||||
::std::vector<int> tstrips;
|
||||
|
||||
// Grid, if present
|
||||
::std::vector<int> grid;
|
||||
int grid_width, grid_height;
|
||||
|
||||
// Other per-vertex properties
|
||||
::std::vector<Color> colors;
|
||||
::std::vector<float> confidences;
|
||||
::std::vector<unsigned> flags;
|
||||
unsigned flag_curr;
|
||||
|
||||
// Computed per-vertex properties
|
||||
::std::vector<vec> normals;
|
||||
::std::vector<vec> pdir1, pdir2;
|
||||
::std::vector<float> curv1, curv2;
|
||||
::std::vector< Vec<4,float> > dcurv;
|
||||
::std::vector<vec> cornerareas;
|
||||
::std::vector<float> pointareas;
|
||||
|
||||
// Bounding structures
|
||||
box bbox;
|
||||
BSphere bsphere;
|
||||
|
||||
// Connectivity structures:
|
||||
// For each vertex, all neighboring vertices
|
||||
::std::vector< ::std::vector<int> > neighbors;
|
||||
// For each vertex, all neighboring faces
|
||||
::std::vector< ::std::vector<int> > 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<Face> 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
|
||||
331
simgear/scene/dem/TriMesh_curvature.cc
Normal file
331
simgear/scene/dem/TriMesh_curvature.cc
Normal file
@@ -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<float> 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<float,3>(w, diag)) {
|
||||
//dprintf("ldltdc failed!\n");
|
||||
continue;
|
||||
}
|
||||
ldltsl<float,3>(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<float,4>(w, d)) {
|
||||
//dprintf("ldltdc failed!\n");
|
||||
continue;
|
||||
}
|
||||
ldltsl<float,4>(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
|
||||
116
simgear/scene/dem/TriMesh_normals.cc
Normal file
116
simgear/scene/dem/TriMesh_normals.cc
Normal file
@@ -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<const float *> 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<float,3>(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
|
||||
@@ -26,8 +26,16 @@
|
||||
#include <ostream>
|
||||
#include <sstream>
|
||||
|
||||
#include <osg/Geode>
|
||||
#include <osg/CullFace>
|
||||
|
||||
#include <simgear/bucket/newbucket.hxx>
|
||||
#include <simgear/math/SGGeometry.hxx>
|
||||
#include <simgear/scene/util/OsgMath.hxx>
|
||||
|
||||
#ifdef ENABLE_GDAL
|
||||
#include <simgear/scene/dem/SGMesh.hxx>
|
||||
#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;
|
||||
|
||||
@@ -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<osg::Group> 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);
|
||||
|
||||
@@ -38,6 +38,7 @@
|
||||
#include "SGReaderWriterBTG.hxx"
|
||||
#include "ReaderWriterSPT.hxx"
|
||||
#include "ReaderWriterSTG.hxx"
|
||||
#include <simgear/scene/dem/ReaderWriterPGT.hxx>
|
||||
|
||||
// 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<simgear::LoadOnlyCallback> g_stgCallbackProx
|
||||
|
||||
osgDB::RegisterReaderWriterProxy<simgear::ReaderWriterSPT> g_readerWriterSPTProxy;
|
||||
simgear::ModelRegistryCallbackProxy<simgear::LoadOnlyCallback> g_sptCallbackProxy("spt");
|
||||
|
||||
osgDB::RegisterReaderWriterProxy<simgear::ReaderWriterPGT> g_readerWriterPGTProxy;
|
||||
simgear::ModelRegistryCallbackProxy<simgear::LoadOnlyCallback> g_pgtCallbackProxy("pgt");
|
||||
}
|
||||
|
||||
void sgUserDataInit( SGPropertyNode *p ) {
|
||||
|
||||
@@ -26,6 +26,10 @@
|
||||
|
||||
#include <simgear/props/props.hxx>
|
||||
|
||||
#ifdef ENABLE_GDAL
|
||||
#include <simgear/scene/dem/SGDem.hxx>
|
||||
#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<SGPropertyNode> _propertyNode;
|
||||
SGSharedPtr<SGMaterialLib> _materialLib;
|
||||
#ifdef ENABLE_GDAL
|
||||
SGSharedPtr<SGDem> _dem;
|
||||
#endif
|
||||
|
||||
osg::Node *(*_load_panel)(SGPropertyNode *);
|
||||
osg::ref_ptr<SGModelData> _model_data;
|
||||
bool _instantiateEffects;
|
||||
|
||||
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user