From Jason Baverage, A

dded GeospatialExtents bounding box class which used doubles
in place of the original usage of osg::BoundingBox.

Added path for computing interpolation elevation data being read from GDAL.
This commit is contained in:
Robert Osfield 2006-05-15 13:12:55 +00:00
parent 4431b381d3
commit 2e0865d0c9
2 changed files with 184 additions and 32 deletions

View File

@ -25,9 +25,9 @@
#ifndef DATASET_H
#define DATASET_H 1
#include <osg/BoundingSphere>
#include <osg/Node>
#include <osg/Matrixd>
#include <osg/BoundingBox>
#include <osg/Image>
#include <osg/Shape>
#include <osg/CoordinateSystemNode>
@ -47,6 +47,101 @@ namespace osgTerrain
#define MAXIMUM_NUMBER_OF_LEVELS 30
class GeospatialExtents
{
public:
osg::Vec2d _min;
osg::Vec2d _max;
inline GeospatialExtents() :
_min(DBL_MAX,DBL_MAX),
_max(DBL_MIN,DBL_MIN) {}
inline GeospatialExtents(double xmin, double ymin, double xmax, double ymax) :
_min(xmin, ymin),
_max(xmax, ymax){}
inline double& xMin() { return _min.x(); }
inline double xMin() const { return _min.x(); }
inline double& yMin() { return _min.y(); }
inline double yMin() const { return _min.y(); }
inline double& xMax() { return _max.x(); }
inline double xMax() const { return _max.x(); }
inline double& yMax() { return _max.y(); }
inline double yMax() const { return _max.y(); }
inline void init()
{
_min.set(DBL_MAX,DBL_MAX);
_max.set(-DBL_MAX,-DBL_MAX);
}
inline bool valid() const
{
return _max.x()>=_min.x() && _max.y()>=_min.y();
}
inline double radius() const
{
return sqrt((radius2()));
}
inline double radius2() const
{
return 0.25f*((_max-_min).length2());
}
GeospatialExtents intersect(const GeospatialExtents& e) const
{
return GeospatialExtents(osg::maximum(xMin(),e.xMin()),osg::maximum(yMin(),e.yMin()),
osg::minimum(xMax(),e.xMax()),osg::minimum(yMax(),e.yMax()));
}
/** Return true if this bounding box intersects the specified bounding box. */
bool intersects(const GeospatialExtents& bb) const
{
return osg::maximum(xMin(),bb.xMin()) <= osg::minimum(xMax(),bb.xMax()) &&
osg::maximum(yMin(),bb.yMin()) <= osg::minimum(yMax(),bb.yMax());
}
void expandBy(const osg::BoundingSphere& sh)
{
if (!sh.valid()) return;
if(sh._center.x()-sh._radius<_min.x()) _min.x() = sh._center.x()-sh._radius;
if(sh._center.x()+sh._radius>_max.x()) _max.x() = sh._center.x()+sh._radius;
if(sh._center.y()-sh._radius<_min.y()) _min.y() = sh._center.y()-sh._radius;
if(sh._center.y()+sh._radius>_max.y()) _max.y() = sh._center.y()+sh._radius;
}
inline void expandBy(const osg::Vec3& v)
{
if(v.x()<_min.x()) _min.x() = v.x();
if(v.x()>_max.x()) _max.x() = v.x();
if(v.y()<_min.y()) _min.y() = v.y();
if(v.y()>_max.y()) _max.y() = v.y();
}
void expandBy(const GeospatialExtents& e)
{
if (!e.valid()) return;
if(e._min.x()<_min.x()) _min.x() = e._min.x();
if(e._max.x()>_max.x()) _max.x() = e._max.x();
if(e._min.y()<_min.y()) _min.y() = e._min.y();
if(e._max.y()>_max.y()) _max.y() = e._max.y();
}
};
class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
{
public:
@ -70,7 +165,7 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
_numValuesY(sp._numValuesY),
_numValuesZ(sp._numValuesZ) {}
SpatialProperties(osg::CoordinateSystemNode* cs, const osg::BoundingBox& extents):
SpatialProperties(osg::CoordinateSystemNode* cs, const GeospatialExtents& extents):
_cs(cs),
_extents(extents),
_numValuesX(0),
@ -100,7 +195,7 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
osg::ref_ptr<osg::CoordinateSystemNode> _cs;
osg::Matrixd _geoTransform;
osg::BoundingBox _extents;
GeospatialExtents _extents;
unsigned int _numValuesX;
unsigned int _numValuesY;
unsigned int _numValuesZ;
@ -141,7 +236,7 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
static SourceData* readData(Source* source);
osg::BoundingBox getExtents(const osg::CoordinateSystemNode* cs) const;
GeospatialExtents getExtents(const osg::CoordinateSystemNode* cs) const;
const SpatialProperties& computeSpatialProperties(const osg::CoordinateSystemNode* cs) const;
@ -149,9 +244,9 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
void read(DestinationData& destination);
void readImage(DestinationData& destination);
void readHeightField(DestinationData& destination);
void readModels(DestinationData& destination);
virtual void readImage(DestinationData& destination);
virtual void readHeightField(DestinationData& destination);
virtual void readModels(DestinationData& destination);
Source* _source;
@ -831,7 +926,7 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
_maxVisibleDistance(FLT_MAX),
_subTileGenerated(false) {}
CompositeDestination(osg::CoordinateSystemNode* cs, const osg::BoundingBox& extents):
CompositeDestination(osg::CoordinateSystemNode* cs, const GeospatialExtents& extents):
SpatialProperties(cs,extents),
_dataSet(0),
_parent(0),
@ -966,7 +1061,7 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
bool mapLatLongsToXYZ() const { return getConvertFromGeographicToGeocentric() && getEllipsoidModel(); }
void setDestinationExtents(const osg::BoundingBox& extents) { _extents = extents; }
void setDestinationExtents(const GeospatialExtents& extents) { _extents = extents; }
void setDestinationGeoTransform(const osg::Matrixd& geoTransform) { _geoTransform = geoTransform; }
@ -1071,7 +1166,7 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
CompositeDestination* createDestinationGraph(CompositeDestination* parent,
osg::CoordinateSystemNode* cs,
const osg::BoundingBox& extents,
const GeospatialExtents& extents,
unsigned int maxImageSize,
unsigned int maxTerrainSize,
unsigned int currentLevel,
@ -1138,7 +1233,7 @@ class OSGTERRAIN_EXPORT DataSet : public osg::Referenced
osg::ref_ptr<osg::EllipsoidModel> _ellipsoidModel;
osg::Matrixd _geoTransform;
osg::BoundingBox _extents;
GeospatialExtents _extents;
std::string _archiveName;
osg::ref_ptr<osgDB::Archive> _archive;
std::string _directory;

View File

@ -333,7 +333,7 @@ DataSet::SourceData* DataSet::SourceData::readData(Source* source)
return 0;
}
osg::BoundingBox DataSet::SourceData::getExtents(const osg::CoordinateSystemNode* cs) const
GeospatialExtents DataSet::SourceData::getExtents(const osg::CoordinateSystemNode* cs) const
{
return computeSpatialProperties(cs)._extents;
}
@ -446,11 +446,11 @@ void DataSet::SourceData::readImage(DestinationData& destination)
{
if (destination._image.valid())
{
osg::BoundingBox s_bb = getExtents(destination._cs.get());
GeospatialExtents s_bb = getExtents(destination._cs.get());
osg::BoundingBox d_bb = destination._extents;
GeospatialExtents d_bb = destination._extents;
osg::BoundingBox intersect_bb(s_bb.intersect(d_bb));
GeospatialExtents intersect_bb(s_bb.intersect(d_bb));
if (!intersect_bb.valid())
{
my_notify(osg::INFO)<<"Reading image but it does not intesection destination - ignoring"<<std::endl;
@ -789,11 +789,11 @@ void DataSet::SourceData::readHeightField(DestinationData& destination)
{
my_notify(osg::INFO)<<"Reading height field"<<std::endl;
osg::BoundingBox s_bb = getExtents(destination._cs.get());
GeospatialExtents s_bb = getExtents(destination._cs.get());
osg::BoundingBox d_bb = destination._extents;
GeospatialExtents d_bb = destination._extents;
osg::BoundingBox intersect_bb(s_bb.intersect(d_bb));
GeospatialExtents intersect_bb(s_bb.intersect(d_bb));
if (!intersect_bb.valid())
{
@ -887,15 +887,71 @@ void DataSet::SourceData::readHeightField(DestinationData& destination)
my_notify(osg::INFO)<<"********* getLinearUnits = "<<getLinearUnits(_cs.get())<<std::endl;
int readWidth = destWidth;
int readHeight = destHeight;
bool doResample = false;
bool interpolateTerrain = true;
if (interpolateTerrain &&
((windowWidth != readWidth) || (windowHeight != readHeight)))
{
readWidth = windowWidth;
readHeight = windowHeight;
doResample = true;
}
// read data into temporary array
float* heightData = new float [ destWidth*destHeight ];
float* heightData = new float [ readWidth*readHeight ];
// raad the data.
osg::HeightField* hf = destination._heightField.get();
my_notify(osg::INFO)<<"reading height field"<<std::endl;
//bandSelected->RasterIO(GF_Read,windowX,_numValuesY-(windowY+windowHeight),windowWidth,windowHeight,floatdata,destWidth,destHeight,GDT_Float32,numBytesPerZvalue,lineSpace);
bandSelected->RasterIO(GF_Read,windowX,_numValuesY-(windowY+windowHeight),windowWidth,windowHeight,heightData,destWidth,destHeight,GDT_Float32,0,0);
//bandSelected->RasterIO(GF_Read,windowX,_numValuesY-(windowY+windowHeight),windowWidth,windowHeight,heightData,destWidth,destHeight,GDT_Float32,0,0);
bandSelected->RasterIO(GF_Read,windowX,_numValuesY-(windowY+windowHeight),windowWidth,windowHeight,heightData,readWidth,readHeight,GDT_Float32,0,0);
if (doResample)
{
my_notify(osg::INFO)<<"Resampling height field."<<std::endl;
float* resampledHeightData = new float[destWidth * destHeight];
for (int j=0; j<destHeight;++j)
{
double t_d = (double)j/((double)destHeight-1);
for (int i=0;i<destWidth;++i)
{
double s_d = (double)i/((double)destWidth-1);
double read_i = s_d * ((double)readWidth-1);
double read_j = t_d * ((double)readHeight-1);
int right = osg::minimum((int)ceil(read_i), readWidth-1);
int left = osg::maximum((int)floor(read_i), 0);
int top = osg::minimum((int)ceil(read_j), readHeight-1);
int bottom = osg::maximum((int)floor(read_j),0);
double v00 = (double)heightData[left+bottom*readWidth];
double v01 = (double)heightData[left+top*readWidth];
double v10 = (double)heightData[right+bottom*readWidth];
double v11 = (double)heightData[right+top*readWidth];
double x_rem = read_i - (int)read_i;
double y_rem = read_j - (int)read_j;
double w00 = (1.0 - y_rem) * (1.0 - x_rem) * v00;
double w01 = (1.0 - y_rem) * x_rem * v01;
double w10 = y_rem * (1.0 - x_rem) * v10;
double w11 = y_rem * x_rem * v11;
float result = (float)(w00 + w01 + w10 + w11);
resampledHeightData[i+j*destWidth] = result;
}
}
delete [] heightData;
heightData = resampledHeightData;
}
my_notify(osg::INFO)<<" scaling height field"<<std::endl;
@ -916,6 +972,7 @@ void DataSet::SourceData::readHeightField(DestinationData& destination)
}
}
delete [] heightData;
}
@ -3828,7 +3885,7 @@ void DataSet::loadSources()
DataSet::CompositeDestination* DataSet::createDestinationGraph(CompositeDestination* parent,
osg::CoordinateSystemNode* cs,
const osg::BoundingBox& extents,
const GeospatialExtents& extents,
unsigned int maxImageSize,
unsigned int maxTerrainSize,
unsigned int currentLevel,
@ -3942,10 +3999,10 @@ DataSet::CompositeDestination* DataSet::createDestinationGraph(CompositeDestinat
{
my_notify(osg::INFO)<<"Need to Divide X + Y for level "<<currentLevel<<std::endl;
// create four tiles.
osg::BoundingBox bottom_left(extents.xMin(),extents.yMin(),extents.zMin(),xCenter,yCenter,extents.zMax());
osg::BoundingBox bottom_right(xCenter,extents.yMin(),extents.zMin(),extents.xMax(),yCenter,extents.zMax());
osg::BoundingBox top_left(extents.xMin(),yCenter,extents.zMin(),xCenter,extents.yMax(),extents.zMax());
osg::BoundingBox top_right(xCenter,yCenter,extents.zMin(),extents.xMax(),extents.yMax(),extents.zMax());
GeospatialExtents bottom_left(extents.xMin(),extents.yMin(),xCenter,yCenter);
GeospatialExtents bottom_right(xCenter,extents.yMin(),extents.xMax(),yCenter);
GeospatialExtents top_left(extents.xMin(),yCenter,xCenter,extents.yMax());
GeospatialExtents top_right(xCenter,yCenter,extents.xMax(),extents.yMax());
destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
cs,
@ -4003,8 +4060,8 @@ DataSet::CompositeDestination* DataSet::createDestinationGraph(CompositeDestinat
my_notify(osg::INFO)<<"Need to Divide X only"<<std::endl;
// create two tiles.
osg::BoundingBox left(extents.xMin(),extents.yMin(),extents.zMin(),xCenter,extents.yMax(),extents.zMax());
osg::BoundingBox right(xCenter,extents.yMin(),extents.zMin(),extents.xMax(),extents.yMax(),extents.zMax());
GeospatialExtents left(extents.xMin(),extents.yMin(),xCenter,extents.yMax());
GeospatialExtents right(xCenter,extents.yMin(),extents.xMax(),extents.yMax());
destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
cs,
@ -4043,8 +4100,8 @@ DataSet::CompositeDestination* DataSet::createDestinationGraph(CompositeDestinat
my_notify(osg::INFO)<<"Need to Divide Y only"<<std::endl;
// create two tiles.
osg::BoundingBox top(extents.xMin(),yCenter,extents.zMin(),extents.xMax(),extents.yMax(),extents.zMax());
osg::BoundingBox bottom(extents.xMin(),extents.yMin(),extents.zMin(),extents.xMax(),yCenter,extents.zMax());
GeospatialExtents top(extents.xMin(),yCenter,extents.xMax(),extents.yMax());
GeospatialExtents bottom(extents.xMin(),extents.yMin(),extents.xMax(),yCenter);
destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
cs,
@ -4140,7 +4197,7 @@ void DataSet::computeDestinationGraphFromSources(unsigned int numLevels)
}
// get the extents of the sources and
osg::BoundingBox extents(_extents);
GeospatialExtents extents(_extents);
if (!extents.valid())
{
for(CompositeSource::source_iterator itr(_sourceGraph.get());itr.valid();++itr)
@ -4148,7 +4205,7 @@ void DataSet::computeDestinationGraphFromSources(unsigned int numLevels)
SourceData* sd = (*itr)->getSourceData();
if (sd)
{
osg::BoundingBox local_extents(sd->getExtents(_intermediateCoordinateSystem.get()));
GeospatialExtents local_extents(sd->getExtents(_intermediateCoordinateSystem.get()));
my_notify(osg::INFO)<<"local_extents = xMin()"<<local_extents.xMin()<<" "<<local_extents.xMax()<<std::endl;
my_notify(osg::INFO)<<" yMin()"<<local_extents.yMin()<<" "<<local_extents.yMax()<<std::endl;
extents.expandBy(local_extents);