976 lines
31 KiB
C++
Executable File
976 lines
31 KiB
C++
Executable File
/**
|
||
* @RDE_HEADER@
|
||
* adapting GeISiR source to SituationMap
|
||
* $Id$
|
||
* $HeadURL$
|
||
* $Author$
|
||
*/
|
||
|
||
|
||
#include "elevation_map.h"
|
||
|
||
#include <QtConcurrent>
|
||
#include <QPolygonF>
|
||
#include <MathUtils/math_const.h>
|
||
#include <QElapsedTimer>
|
||
|
||
#include <gdal.h>
|
||
|
||
#include <iostream>
|
||
#include <limits>
|
||
#include <memory>
|
||
#include <math.h>
|
||
|
||
namespace SituationMap
|
||
{
|
||
|
||
class BoundingRect
|
||
{
|
||
double m_top = 0.0;
|
||
double m_bottom = 0.0;
|
||
double m_left = 0.0;
|
||
double m_right = 0.0;
|
||
|
||
void swapBounds()
|
||
{
|
||
if (m_top < m_bottom)
|
||
std::swap(m_bottom, m_top);
|
||
if (m_right < m_left)
|
||
std::swap(m_left, m_right);
|
||
} // void swapBounds()
|
||
|
||
public:
|
||
|
||
BoundingRect() = default;
|
||
|
||
BoundingRect(double left, double bottom, double right, double top)
|
||
: m_top(top), m_bottom(bottom), m_left(left), m_right(right)
|
||
{
|
||
swapBounds();
|
||
} // BoundingRect( double left, double bottom, double right, double top)
|
||
|
||
BoundingRect(const BoundingRect &) = default;
|
||
BoundingRect &operator =(const BoundingRect &) = default;
|
||
|
||
~BoundingRect() = default;
|
||
|
||
inline void operator()(double left, double bottom, double right, double top)
|
||
{
|
||
m_top = top;
|
||
m_bottom = bottom;
|
||
m_left = left;
|
||
m_right = right;
|
||
|
||
swapBounds();
|
||
} // inline void operator() (double left, double bottom, double right, double top)
|
||
|
||
double top() const { return m_top; }
|
||
double bottom() const { return m_bottom; }
|
||
double left() const { return m_left; }
|
||
double right() const { return m_right; }
|
||
|
||
bool contains(const QPointF &p) const { return ((p.x() >= m_left) && (p.x() <= m_right) && (p.y() >= m_bottom) && (p.y() <= m_top)); }
|
||
|
||
QRectF toQRectF() const
|
||
{
|
||
return QRectF(QPointF(m_left, m_top), QPointF(m_right, m_bottom));
|
||
} // toQRectF
|
||
}; // BoundingRect
|
||
|
||
|
||
|
||
struct ElevationMap::Internal
|
||
{
|
||
BoundingRect _location; // Lage der Karte in Lat/Lon
|
||
|
||
struct PointDouble
|
||
{
|
||
void setX(double x) { _x = x;}
|
||
void setY(double y) { _y = y;}
|
||
double x() { return _x;}
|
||
double y() { return _y;}
|
||
double _x{0.0};
|
||
double _y{0.0};
|
||
};
|
||
PointDouble _elexelSize; // LatLon Size of 1 elevation sample, redundat
|
||
|
||
int _totalColumns = 0;
|
||
int _totalRows = 0;
|
||
|
||
double _toMeterFactor = 1.0;
|
||
std::vector<int> _elevationField;
|
||
|
||
template<class BufferType> bool copyBufferData(GDALRasterBandH hBand, int xBlkSize, int yBlkSize)
|
||
{
|
||
// NLDEVCTRL-804: Konsistenz d. Z-Werte zwischen 2D/3D (also ElevationMap vs. ESRI::BaseSurface)
|
||
// This ElevationMap was using an approach similar to following.
|
||
// Now it´s more or less identical and the difference using 'Altmark_dt1.tif'
|
||
// are gone. AHTC is also fine :-)
|
||
/*
|
||
CPLErr GetHistogram( GDALRasterBand *poBand, GUIntBig *panHistogram )
|
||
{
|
||
memset( panHistogram, 0, sizeof(GUIntBig) * 256 );
|
||
|
||
CPLAssert( poBand->GetRasterDataType() == GDT_Byte );
|
||
|
||
int nXBlockSize, nYBlockSize;
|
||
|
||
poBand->GetBlockSize( &nXBlockSize, &nYBlockSize );
|
||
int nXBlocks = (poBand->GetXSize() + nXBlockSize - 1) / nXBlockSize;
|
||
int nYBlocks = (poBand->GetYSize() + nYBlockSize - 1) / nYBlockSize;
|
||
|
||
GByte *pabyData = (GByte *) CPLMalloc(nXBlockSize * nYBlockSize);
|
||
|
||
for( int iYBlock = 0; iYBlock < nYBlocks; iYBlock++ )
|
||
{
|
||
for( int iXBlock = 0; iXBlock < nXBlocks; iXBlock++ )
|
||
{
|
||
int nXValid, nYValid;
|
||
|
||
poBand->ReadBlock( iXBlock, iYBlock, pabyData );
|
||
|
||
// Compute the portion of the block that is valid
|
||
// for partial edge blocks.
|
||
poBand->GetActualBlockSize(iXBlock, iYBlock, &nXValid, &nYValid)
|
||
|
||
// Collect the histogram counts.
|
||
for( int iY = 0; iY < nYValid; iY++ )
|
||
{
|
||
for( int iX = 0; iX < nXValid; iX++ )
|
||
{
|
||
panHistogram[pabyData[iX + iY * nXBlockSize]] += 1;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
Source https://gdal.org/api/gdalrasterband_cpp.html
|
||
*/
|
||
|
||
const size_t bufSize = xBlkSize * yBlkSize;
|
||
std::unique_ptr<BufferType[]> bufferPtr(new BufferType[bufSize]);
|
||
|
||
int nXBlocks = (GDALGetRasterBandXSize(hBand) + xBlkSize - 1) / xBlkSize;
|
||
int nYBlocks = (GDALGetRasterBandYSize(hBand) + yBlkSize - 1) / yBlkSize;
|
||
|
||
for (int iYBlock = 0; iYBlock < nYBlocks; ++iYBlock)
|
||
{
|
||
for (int iXBlock = 0; iXBlock < nXBlocks; ++iXBlock)
|
||
{
|
||
|
||
CPLErr err = GDALReadBlock( hBand, iXBlock, iYBlock, bufferPtr.get() ); // Gdal liest spaltenweise !! x <-> y
|
||
if (err != CE_None)
|
||
{
|
||
return false;
|
||
}
|
||
|
||
int nXValid, nYValid;
|
||
|
||
// Unfortunately GDALGetActualBlockSize is not available at Suse 42.3 (gdal 2.1.3)
|
||
// so check bounds using _totalColumns, _totalRows:
|
||
#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(2, 2, 0)
|
||
# define HAS_ACTUAL_BLOCKSIZE
|
||
#endif
|
||
|
||
#ifdef HAS_ACTUAL_BLOCKSIZE
|
||
GDALGetActualBlockSize(hBand, iXBlock, iYBlock, &nXValid, &nYValid);
|
||
#else
|
||
nXValid = xBlkSize;
|
||
nYValid = yBlkSize;
|
||
#endif // HAS_ACTUAL_BLOCKSIZE
|
||
|
||
|
||
const int curYOffset = (iYBlock * yBlkSize);
|
||
const int curXOffset = (iXBlock * xBlkSize);
|
||
|
||
for (int iY = 0; iY < nYValid; ++iY)
|
||
{
|
||
const int row = iY + curYOffset; // Line
|
||
|
||
#ifndef HAS_ACTUAL_BLOCKSIZE
|
||
if (row >= _totalRows)
|
||
continue;
|
||
#endif // HAS_ACTUAL_BLOCKSIZE
|
||
|
||
for (int iX = 0; iX < nXValid; ++iX)
|
||
{
|
||
// identical to gdallocationinfo
|
||
const int column = iX + curXOffset; // Pixel
|
||
|
||
#ifndef HAS_ACTUAL_BLOCKSIZE
|
||
if (column >= _totalColumns)
|
||
continue;
|
||
#endif // HAS_ACTUAL_BLOCKSIZE
|
||
|
||
int bi = iX + (iY * xBlkSize);
|
||
const int z = static_cast<int>( static_cast<double>(bufferPtr[bi]) / _toMeterFactor );
|
||
|
||
_elevationField[column + row * _totalColumns] = z;
|
||
} // for
|
||
} // for
|
||
|
||
} // for
|
||
} // for
|
||
|
||
return true;
|
||
|
||
} // template<class BufferType> void copyBufferData(BufferType * bufferPtr, int xBlkSize, int yBlkSize)
|
||
|
||
}; // struct ElevationMap::Internal
|
||
|
||
|
||
ElevationMap::ElevationMap(const QFileInfo &fileInfo)
|
||
: m_internal(new Internal())
|
||
, m_aproximation(true)
|
||
, m_valid(false)
|
||
, m_logger(getLogger("ElevationMap"))
|
||
{
|
||
GDALAllRegister();
|
||
|
||
m_fileName = fileInfo.absoluteFilePath();
|
||
|
||
INFO_LOG(m_logger) << " ctor ElevationMap: map of file '" << m_fileName << "'"
|
||
<< ", this=" << this;
|
||
|
||
if (!fileInfo.exists())
|
||
{
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " elevation data file \'" << m_fileName << "\' doesn't exists.";
|
||
return;
|
||
}
|
||
|
||
if (!fileInfo.isReadable())
|
||
{
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " elevation data file \'" << m_fileName << "\' isn't readable.";
|
||
return;
|
||
}
|
||
|
||
auto cb = [this]()
|
||
{
|
||
return readFromFile();
|
||
};
|
||
auto cbFinished = [this]()
|
||
{
|
||
NOTICE_LOG(m_logger) << " Finished reading elevation map of file '" << m_fileName << "'"
|
||
<< ", this=" << this;
|
||
};
|
||
QObject::connect(&m_gdalReadWatcher, &QFutureWatcher<bool>::finished,
|
||
&m_gdalReadWatcher, cbFinished);
|
||
m_gdalReadWatcher.setFuture(QtConcurrent::run(cb));
|
||
// sync init, takes ~800ms during sitmap loading 72MB file
|
||
// QElapsedTimer t;
|
||
// t.start();
|
||
// m_gdalReadWatcher.waitForFinished();
|
||
// NOTICE_LOG(m_logger) << "Time loading elevation map: " << t.elapsed();
|
||
} // ElevationMap::ElevationMap( const QFile &file )
|
||
|
||
|
||
ElevationMap::~ElevationMap()
|
||
{
|
||
TRACE_LOG(m_logger) << LOGFUNCTION << " START";
|
||
m_gdalReadWatcher.cancel();
|
||
m_gdalReadWatcher.waitForFinished();
|
||
|
||
delete m_internal;
|
||
m_internal = nullptr;
|
||
TRACE_LOG(m_logger) << LOGFUNCTION << " END";
|
||
} // ElevationMap::~ElevationMap
|
||
|
||
bool ElevationMap::waitForInitialization()
|
||
{
|
||
m_gdalReadWatcher.waitForFinished();
|
||
return m_valid;
|
||
}
|
||
|
||
bool ElevationMap::isValid() const
|
||
{
|
||
// N. Thoden, 27.04.2020
|
||
// isFinished seems to be unreliable:
|
||
// if mainThread request ::isValid() diretly after ::waitForInitialization()
|
||
// you get false (by 'if (!m_gdalReadWatcher.isFinished())'
|
||
// but m_valid is true. So skip the additional check.
|
||
// if (!m_gdalReadWatcher.isFinished())
|
||
// return false;
|
||
|
||
return m_valid;
|
||
}
|
||
|
||
QString ElevationMap::filename() const
|
||
{
|
||
return m_fileName;
|
||
}
|
||
|
||
bool ElevationMap::readFromFile()
|
||
{
|
||
m_valid = false;
|
||
|
||
GDALDatasetH hDataset = GDALOpen(m_fileName.toLocal8Bit().constData(), GA_ReadOnly);
|
||
|
||
if (nullptr == hDataset)
|
||
{
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " GDALOpen failed, may be calling of GDALAllRegister is missing";
|
||
return false;
|
||
}
|
||
|
||
// Detect basic parameters
|
||
//
|
||
GDALDriverH hDriver = GDALGetDatasetDriver(hDataset);
|
||
int nRasterXSize = GDALGetRasterXSize(hDataset);
|
||
int nRasterYSize = GDALGetRasterYSize(hDataset);
|
||
double adfGeoTransform[6];
|
||
CPLErr err = GDALGetGeoTransform(hDataset, adfGeoTransform); // XY Spatial Referenz WGS84 (4326) expected!!!
|
||
|
||
if (err != CE_None)
|
||
{
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " GDAL failed to get geo raster transform.";
|
||
GDALClose(hDataset);
|
||
return false;
|
||
}
|
||
|
||
|
||
if (nullptr == hDriver) { GDALClose(hDataset); return false; }
|
||
if (nRasterXSize <= 0) { GDALClose(hDataset); return false; }
|
||
if (nRasterYSize <= 0) { GDALClose(hDataset); return false; }
|
||
|
||
double pixWidth = adfGeoTransform[1];
|
||
double pixHeight = adfGeoTransform[5];
|
||
|
||
auto toRAD = [](double x) -> double { return x * MathConst::pi / 180.0; };
|
||
|
||
|
||
m_internal->_location(toRAD(adfGeoTransform[0]) // adfGeoTransform is in degrees
|
||
, toRAD(adfGeoTransform[3])
|
||
, toRAD(adfGeoTransform[0] + double(nRasterXSize) * pixWidth + double(nRasterYSize) * adfGeoTransform[2])
|
||
, toRAD(adfGeoTransform[3] + double(nRasterYSize) * pixHeight + double(nRasterXSize) * adfGeoTransform[4]));
|
||
|
||
m_internal->_elexelSize.setX(fabs(toRAD(pixWidth)));
|
||
m_internal->_elexelSize.setY(fabs(toRAD(pixHeight)));
|
||
|
||
GDALRasterBandH hBand = GDALGetRasterBand(hDataset, 1);
|
||
|
||
if (nullptr == hBand)
|
||
{
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " GDAL failed to read raster band.";
|
||
GDALClose(hDataset);
|
||
return false;
|
||
}
|
||
|
||
int xBlkSize;
|
||
int yBlkSize;
|
||
GDALGetBlockSize(hBand, &xBlkSize, &yBlkSize);
|
||
|
||
if (xBlkSize <0 || yBlkSize < 0)
|
||
{
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " no buffer, next operations will possible fail....";
|
||
return false;
|
||
}
|
||
|
||
m_internal->_totalColumns = nRasterXSize;
|
||
m_internal->_totalRows = nRasterYSize;
|
||
m_internal->_elevationField.resize(nRasterXSize * nRasterYSize);
|
||
|
||
const GDALDataType type = GDALGetRasterDataType(hBand);
|
||
switch (type)
|
||
{
|
||
case GDT_UInt16:
|
||
m_internal->_toMeterFactor = 1.0;
|
||
m_valid = m_internal->copyBufferData<unsigned short int> ( hBand, xBlkSize, yBlkSize );
|
||
break;
|
||
case GDT_Int16:
|
||
m_internal->_toMeterFactor = 1.0;
|
||
m_valid = m_internal->copyBufferData<short int> ( hBand, xBlkSize, yBlkSize );
|
||
break;
|
||
case GDT_Float32:
|
||
m_internal->_toMeterFactor = 0.01;
|
||
m_valid = m_internal->copyBufferData<float> ( hBand, xBlkSize, yBlkSize );
|
||
break;
|
||
case GDT_Float64:
|
||
m_internal->_toMeterFactor = 0.01;
|
||
m_valid = m_internal->copyBufferData<double> ( hBand, xBlkSize, yBlkSize );
|
||
break;
|
||
default:
|
||
{
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " data type isn't supported type or a floating point type"
|
||
<< " (type=" << type << ")";
|
||
GDALClose(hDataset);
|
||
return false;
|
||
}
|
||
} // switch
|
||
|
||
if (!m_valid)
|
||
{
|
||
CPLErrorNum errNum = CPLGetLastErrorNo();
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << " GDAL failed to read data block (error no " << int(errNum) << ")";
|
||
}
|
||
|
||
GDALClose(hDataset);
|
||
|
||
INFO_LOG(m_logger) << LOGFUNCTION << " loaded, filename =" << m_fileName << "'"
|
||
<< ", m_valid=" << m_valid;
|
||
return m_valid;
|
||
|
||
} // bool ElevationMap::readFromFile(const QString & fileName)
|
||
|
||
|
||
|
||
|
||
QRectF ElevationMap::extent() const
|
||
{
|
||
if (!isValid())
|
||
return QRectF();
|
||
|
||
return m_internal->_location.toQRectF();
|
||
}
|
||
|
||
QPair<double, double> ElevationMap::elexelSize() const
|
||
{
|
||
return {m_internal->_elexelSize._x, m_internal->_elexelSize._y};
|
||
}
|
||
|
||
|
||
int ElevationMap::getElevation(int column, int row) const
|
||
{
|
||
row = m_internal->_totalRows - row;
|
||
const size_t index = column + row * m_internal->_totalColumns;
|
||
|
||
if (m_internal->_elevationField.empty() || m_internal->_elevationField.size() < index)
|
||
return std::numeric_limits<int>::min();
|
||
|
||
return m_internal->_elevationField.at(index);
|
||
} // int ElevationMap::getElevation( int column, int row)
|
||
|
||
|
||
double ElevationMap::elevation(double lon, double lat) const
|
||
{
|
||
if (!isValid())
|
||
return 0;
|
||
|
||
return getElevation(lon, lat) * m_internal->_toMeterFactor;
|
||
} // double ElevationMap::elevation( double lon, double lat)
|
||
|
||
|
||
double ElevationMap:: getElevation(double lon, double lat) const
|
||
{
|
||
if (!m_internal->_location.contains(QPointF(lon, lat)))
|
||
{
|
||
return std::numeric_limits<int>::min(); // ausserhalb der Höhenkarte !!
|
||
}
|
||
double x = (lon - m_internal->_location.left()) / m_internal->_elexelSize.x();
|
||
int column = floor(x);
|
||
|
||
double y = (lat - m_internal->_location.bottom()) / m_internal->_elexelSize.y();
|
||
int row = floor(y);
|
||
|
||
int elevation = getElevation(column, row);
|
||
if (elevation == std::numeric_limits<int>::min())
|
||
return elevation;
|
||
|
||
if (!m_aproximation)
|
||
return elevation;
|
||
|
||
/*
|
||
** Approximation: e0, e1, e2, e3 Hoehenwerte der 4 benachbarten Stützpunkte
|
||
**
|
||
** + e3 + e2
|
||
** y
|
||
** ^ + er er = e0 + x*(e1-e0) + y*(e3-e0) + x*y(e2-e3+e0-e1)
|
||
** |
|
||
** |
|
||
** |
|
||
** | + e0 + e1
|
||
** -----> x
|
||
**
|
||
** Besonderheiten: die Stützpunkte liegen definitiongemaess jeweils in der Mitte des Rasters,
|
||
** so dass je nach Groesse von x, y andere Nachbarpunkte zu berücksichtigen sind, und x, sowie y entsprechend
|
||
** angepasst werden muesen.
|
||
** Bei Randrasterpunkten kann keine Interpolation erfolgen, so dass hier entsprechend extrapoliert wird.
|
||
*/
|
||
|
||
x = (x - (double)column); // Residuum normiert auf Stützpunktabstand
|
||
y = (y - (double)row); // Residuum normiert auf Stützpunktabstand
|
||
|
||
x -= 0.5; // Mitte berücksichtigen
|
||
y -= 0.5; // Mitte berücksichtigen
|
||
|
||
if (x < 0.0) { column--; x += 1.0;} // linke Spalte zur Interpolation nutzen
|
||
|
||
if (y < 0.0) { row--; y += 1.0;} // untere Reihe zur Interpolation nutzen
|
||
|
||
// Randrasterpunkte berücksichtigen: extrapolieren statt interpolieren
|
||
if (column <= -1) { x = -x; column = 0;}
|
||
else if (column >= (m_internal->_totalColumns - 1)) { x += 1.0; column = m_internal->_totalColumns - 2;}
|
||
|
||
if (row == -1) { y = -y; row = 0;}
|
||
else if (row >= (m_internal->_totalRows - 1)) { y += 1.0; row = m_internal->_totalRows - 2;}
|
||
|
||
// Die Approximation:
|
||
double substitute = elevation;
|
||
double e0 = getOrSubstituteElevation(column, row, substitute);
|
||
double e1 = getOrSubstituteElevation(column + 1, row, substitute);
|
||
double e2 = getOrSubstituteElevation(column + 1, row + 1, substitute);
|
||
double e3 = getOrSubstituteElevation(column, row + 1, substitute);
|
||
|
||
double er = e0 + x * (e1 - e0) + y * (e3 - e0) + x * y * (e2 - e3 + e0 - e1);
|
||
|
||
return er; // floor(er+0.5);; runden herausgenommen 19.8.2008 JT
|
||
} // double ElevationMap:: getElevation( double lon, double lat)
|
||
|
||
|
||
double ElevationMap::getOrSubstituteElevation(int column, int row, double substitute) const
|
||
{
|
||
int elevation = getElevation(column, row);
|
||
|
||
if (elevation != std::numeric_limits<int>::min())
|
||
return elevation;
|
||
|
||
int substCount = 0;
|
||
double calcSubstitute = 0.0;
|
||
|
||
elevation = getElevation(column + 1, row);
|
||
if (elevation != std::numeric_limits<int>::min()) {substCount++; calcSubstitute += elevation; }
|
||
|
||
elevation = getElevation(column - 1, row);
|
||
if (elevation != std::numeric_limits<int>::min()) {substCount++; calcSubstitute += elevation; }
|
||
|
||
elevation = getElevation(column, row + 1);
|
||
if (elevation != std::numeric_limits<int>::min()) {substCount++; calcSubstitute += elevation; }
|
||
|
||
elevation = getElevation(column, row - 1);
|
||
if (elevation != std::numeric_limits<int>::min()) {substCount++; calcSubstitute += elevation; }
|
||
|
||
if (substCount)
|
||
substitute = calcSubstitute / substCount;
|
||
|
||
return substitute;
|
||
} // double ElevationMap::getOrSubstituteElevation( int column, int row, double substitute)
|
||
|
||
double ElevationMap::distance(const QPointF &start, const QPointF &end) const
|
||
{
|
||
if (!isValid())
|
||
return 0;
|
||
|
||
if (!m_internal->_location.contains(start) || !m_internal->_location.contains(end))
|
||
{
|
||
return -1.0; // ausserhalb{ der Höhenkarte !!
|
||
}
|
||
|
||
double start_z = getElevation(start.x(), start.y());
|
||
if (start_z <= std::numeric_limits<int>::min()) return -1.0;
|
||
|
||
double end_z = getElevation(end.x(), end.y());
|
||
if (end_z <= std::numeric_limits<int>::min()) return -1.0;
|
||
|
||
const QPointF &dist = distance_in_meter_WGS84(start, end);
|
||
|
||
double dz = end_z - start_z;
|
||
double d = sqrt(dist.x() * dist.x() + dist.y() * dist.y() + dz * dz) ;
|
||
|
||
return d;
|
||
} // double ElevationMap::distance( const QPointF &start, const QPointF &end )
|
||
|
||
|
||
double ElevationMap::profileDistance(const QPointF &start, const QPointF &end) const
|
||
{
|
||
if (!isValid())
|
||
return 0;
|
||
|
||
if (!m_internal->_location.contains(start) || !m_internal->_location.contains(end))
|
||
{
|
||
return -1.0; // ausserhalb der Höhenkarte !!
|
||
}
|
||
|
||
Vec3dVector elevProfile;
|
||
getSamplePoints(start, end, elevProfile);
|
||
|
||
for (Vec3dVector::iterator itr = elevProfile.begin(); itr != elevProfile.end(); ++itr)
|
||
itr->setZ(getElevation(itr->x(), itr->y()));
|
||
|
||
// Einzellaengen aufsummieren, Stützstellen mit undefinierten Hoehen werden ausgelassen!
|
||
Point3D curPos(*elevProfile.begin());
|
||
double sumLng = 0.0;
|
||
for (Vec3dVector::iterator itr = elevProfile.begin() + 1; itr != elevProfile.end(); ++itr)
|
||
{
|
||
if ((*itr).z() > std::numeric_limits<int>::min())
|
||
{
|
||
if (curPos.z() > std::numeric_limits<int>::min())
|
||
{
|
||
const QPointF &dist = distance_in_meter_WGS84(QPointF(curPos.x(), curPos.y()), QPointF(itr->x(), itr->y()));
|
||
double dz = itr->z() - curPos.z();
|
||
sumLng += sqrt(dist.x() * dist.x() + dist.y() * dist.y() + dz * dz) ;
|
||
}
|
||
else
|
||
{
|
||
double dx = itr->x() - curPos.x();
|
||
double dy = itr->y() - curPos.y();
|
||
sumLng += sqrt(dx * dx + dy * dy);
|
||
}
|
||
curPos = *itr;
|
||
}
|
||
}
|
||
return sumLng;
|
||
} // double ElevationMap::profileDistance( const QPointF &start, const QPointF &end)
|
||
|
||
|
||
unsigned int ElevationMap::elevationProfile(const QPointF &start, const QPointF &end, Vec3dVector &elevProfile) const
|
||
{
|
||
unsigned int startSize = elevProfile.size();
|
||
|
||
getSamplePoints(start, end, elevProfile);
|
||
for (Vec3dVector::iterator itr = elevProfile.begin() + startSize; itr != elevProfile.end(); ++itr)
|
||
itr->setZ(getElevation(itr->x(), itr->y()));
|
||
|
||
return (elevProfile.size() - startSize);
|
||
} // unsigned int ElevationMap::elevationProfile( const QPointF &start, const QPointF &end, Vec3dVector &elevProfile)
|
||
|
||
bool ElevationMap::contains(const QPointF &point) const
|
||
{
|
||
return m_internal->_location.contains(point);
|
||
}
|
||
|
||
bool ElevationMap::intersects(const QPolygonF &polygon) const
|
||
{
|
||
const BoundingRect &bbox = m_internal->_location;
|
||
|
||
// maybe QPointF uses float, but we can accept double to be reduced to float:
|
||
const QPolygonF &ourBorder = QPolygonF({QPointF(bbox.left() , bbox.bottom()),
|
||
QPointF(bbox.left() , bbox.top()),
|
||
QPointF(bbox.right(), bbox.top()),
|
||
QPointF(bbox.right(), bbox.bottom())});
|
||
|
||
#if QT_VERSION < QT_VERSION_CHECK(5,10,0)
|
||
return ourBorder.intersected(polygon).size() > 0;
|
||
#else
|
||
return ourBorder.intersects(polygon);
|
||
#endif
|
||
}
|
||
|
||
void ElevationMap::getSamplePoints(const QPointF &start, const QPointF &end, Vec3dVector &sPnts) const
|
||
{
|
||
sPnts.clear();
|
||
sPnts.push_back(Point3D::fromQPointF(start));
|
||
|
||
bool endIsClipped = false;
|
||
QPointF cStart(start);
|
||
QPointF cEnd(end);
|
||
|
||
if (!m_internal->_location.contains(start) || !m_internal->_location.contains(end))
|
||
{
|
||
int res = clipp2Map(start, end, cStart, cEnd);
|
||
|
||
if (res == 0) // no intersection
|
||
{
|
||
sPnts.push_back(Point3D::fromQPointF(end));
|
||
return;
|
||
}
|
||
if (res & 1) // start clipped
|
||
sPnts.push_back(Point3D::fromQPointF(cStart));
|
||
if (res & 2) // end clipped
|
||
endIsClipped = true;
|
||
}
|
||
|
||
QPointF delta(cEnd - cStart);
|
||
|
||
// Berechne alle Schnittpunkte des Vektors( cStart, cEnd) mit den x-Grenzen der Hoehenstuetzrechtecke.
|
||
// Die Stützpunkte liegen definitiongemaess jeweils in der Mitte des Raster-Rechtecks
|
||
int startCol = floor((cStart.x() - m_internal->_location.left()) / m_internal->_elexelSize.x());
|
||
if (startCol < 0) startCol = 0;
|
||
else if (startCol >= m_internal->_totalColumns) startCol = m_internal->_totalColumns - 1;
|
||
|
||
int endCol = floor((cEnd.x() - m_internal->_location.left()) / m_internal->_elexelSize.x());
|
||
if (endCol < 0) endCol = 0;
|
||
else if (endCol >= m_internal->_totalColumns) endCol = m_internal->_totalColumns - 1;
|
||
|
||
int incCol = 0;
|
||
if (endCol > startCol) incCol = 1;
|
||
else if (endCol < startCol) incCol = -1;
|
||
|
||
double m = delta.y() / delta.x();
|
||
Vec2dVector xCuts;
|
||
xCuts.reserve(std::abs(endCol - startCol));
|
||
for (int col = startCol; col != endCol; col += incCol)
|
||
{
|
||
double x = m_internal->_location.left() + ((col + 0.5) + (incCol * 0.5)) * m_internal->_elexelSize.x();
|
||
double y = cStart.y() + (x - cStart.x()) * m;
|
||
xCuts.push_back(QPointF(x, y));
|
||
}
|
||
|
||
// Berechne alle Schnittpunkte des Vektors( cStart, cEnd) mit den y-Grenzen der Hoehenstuetzrechtecke.
|
||
// Die Stützpunkte liegen definitiongemaess jeweils in der Mitte des Raster-Rechtecks
|
||
int startRow = floor((cStart.y() - m_internal->_location.bottom()) / m_internal->_elexelSize.y());
|
||
if (startRow < 0) startRow = 0;
|
||
else if (startRow >= m_internal->_totalRows) startRow = m_internal->_totalRows - 1;
|
||
|
||
int endRow = floor((end.y() - m_internal->_location.bottom()) / m_internal->_elexelSize.y());
|
||
if (endRow < 0) endRow = 0;
|
||
else if (endRow >= m_internal->_totalRows) endRow = m_internal->_totalRows - 1;
|
||
|
||
int incRow = 0;
|
||
if (endRow > startRow) incRow = 1;
|
||
else if (endRow < startRow) incRow = -1;
|
||
|
||
m = delta.x() / delta.y();
|
||
Vec2dVector yCuts;
|
||
yCuts.reserve(std::abs(endRow - startRow));
|
||
for (int row = startRow; row != endRow; row += incRow)
|
||
{
|
||
double y = m_internal->_location.bottom() + ((row + 0.5) + (incRow * 0.5)) * m_internal->_elexelSize.y();
|
||
double x = cStart.x() + (y - cStart.y()) * m;
|
||
yCuts.push_back(QPointF(x, y));
|
||
}
|
||
|
||
// xCuts und yCuts zusammenführen, unter Berücksichtigung der Reihenfolge
|
||
Vec2dVector mergeCuts;
|
||
Vec2dVector::iterator xCutItr = xCuts.begin();
|
||
Vec2dVector::iterator yCutItr = yCuts.begin();
|
||
|
||
while ((xCutItr != xCuts.end()) && (yCutItr != yCuts.end()))
|
||
{
|
||
if (incCol > 0)
|
||
{
|
||
if ((*xCutItr).x() < (*yCutItr).x()) { mergeCuts.push_back(*xCutItr); ++xCutItr; }
|
||
else { mergeCuts.push_back(*yCutItr); ++yCutItr; }
|
||
}
|
||
else
|
||
{
|
||
if ((*xCutItr).x() > (*yCutItr).x()) { mergeCuts.push_back(*xCutItr); ++xCutItr; }
|
||
else { mergeCuts.push_back(*yCutItr); ++yCutItr; }
|
||
}
|
||
}
|
||
while ((xCutItr != xCuts.end())) { mergeCuts.push_back(*xCutItr); ++xCutItr; }
|
||
while ((yCutItr != yCuts.end())) { mergeCuts.push_back(*yCutItr); ++yCutItr; }
|
||
|
||
// mergeCuts enthält sequenziell alle Schnittpunkte des Vektors mit den Hoehengridkanten
|
||
// Die 2 Schnittpunkte pro Hoehenwert werden gemittelt zu einem zusammengefasst
|
||
if (mergeCuts.size() > 1)
|
||
{
|
||
for (Vec2dVector::iterator itr = mergeCuts.begin() + 1; itr != mergeCuts.end(); ++itr)
|
||
sPnts.push_back(Point3D::fromQPointF((*(itr - 1) + *itr) * 0.5));
|
||
}
|
||
sPnts.push_back(Point3D::fromQPointF(cEnd));
|
||
|
||
// DEBUG_LOG(m_logger) << " xCuts.si " << xCuts.size()
|
||
// << ", yCuts.si " << yCuts.size()
|
||
// << ", mergeCuts.si " << mergeCuts.size();
|
||
|
||
|
||
if (endIsClipped) sPnts.push_back(Point3D::fromQPointF(end));
|
||
} // void ElevationMap::getSamplePoints( const QPointF &start, const QPointF &end, Vec3dVector & sPnts)
|
||
|
||
|
||
// start und/oder end liegt ausserhalb !! in Lat/Lon !!
|
||
int ElevationMap::clipp2Map(const QPointF &start, const QPointF &end, QPointF &cStart, QPointF &cEnd) const
|
||
{
|
||
Vec2dVector clippVec;
|
||
|
||
if (end.x() != start.x())
|
||
{
|
||
double x = m_internal->_location.right();
|
||
double m = (end.y() - start.y()) / (end.x() - start.x());
|
||
if (isInIntervall(x, start.x(), end.x()))
|
||
{
|
||
double yc = start.y() + (x - start.x()) * m;
|
||
if ((yc <= m_internal->_location.top()) && (yc >= m_internal->_location.bottom()))
|
||
clippVec.push_back(QPointF(x, yc));
|
||
}
|
||
x = m_internal->_location.left();
|
||
if (isInIntervall(x, start.x(), end.x()))
|
||
{
|
||
double yc = start.y() + (x - start.x()) * (end.y() - start.y()) / (end.x() - start.x());
|
||
if ((yc <= m_internal->_location.top()) && (yc >= m_internal->_location.bottom()))
|
||
clippVec.push_back(QPointF(x, yc));
|
||
}
|
||
}
|
||
|
||
if (end.y() != start.y())
|
||
{
|
||
double y = m_internal->_location.top();
|
||
double m = (end.x() - start.x()) / (end.y() - start.y());
|
||
if (isInIntervall(y, start.y(), end.y()))
|
||
{
|
||
double xc = start.x() + (y - start.y()) * m;
|
||
if ((xc < m_internal->_location.right()) && (xc > m_internal->_location.left()))
|
||
clippVec.push_back(QPointF(xc, y));
|
||
}
|
||
y = m_internal->_location.bottom();
|
||
if (isInIntervall(y, start.y(), end.y()))
|
||
{
|
||
double xc = start.x() + (y - start.y()) * m;
|
||
if ((xc < m_internal->_location.right()) && (xc > m_internal->_location.left()))
|
||
clippVec.push_back(QPointF(xc, y));
|
||
}
|
||
}
|
||
int size = (int)clippVec.size();
|
||
if (size == 0)
|
||
{
|
||
return 0; // komlett ausserhalb ( oder komplett innerhalb)
|
||
}
|
||
|
||
if (size == 1)
|
||
{
|
||
if (m_internal->_location.contains(start))
|
||
{
|
||
cEnd = clippVec[0]; return 2;
|
||
}
|
||
else
|
||
{
|
||
cStart = clippVec[0]; return 1;
|
||
}
|
||
}
|
||
|
||
if (size == 2)
|
||
{
|
||
if (squareLength(clippVec[0] - start) < squareLength(clippVec[1] - start))
|
||
{
|
||
cStart = clippVec[0];
|
||
cEnd = clippVec[1];
|
||
}
|
||
else
|
||
{
|
||
cStart = clippVec[1];
|
||
cEnd = clippVec[0];
|
||
}
|
||
return 3;
|
||
}
|
||
|
||
ERROR_LOG(m_logger) << LOGFUNCTION << "Internal Clipp Error!";
|
||
|
||
return 0;
|
||
} // int ElevationMap::clipp2Map( const QPointF &start, const QPointF &end, QPointF &cStart, QPointF &cEnd)
|
||
|
||
|
||
// Line of Sight Profile: berücksichtigt Erdkrümmung und WGS84 Erdellipsoid
|
||
void ElevationMap::LOSProfile(const QPointF &start, const QPointF &end, LOSVector &losProfile)
|
||
{
|
||
m_gdalReadWatcher.waitForFinished();
|
||
|
||
if (!m_internal->_location.contains(start)) return ; // ausserhalb der Höhenkarte !!
|
||
|
||
bool clipped = false;
|
||
double totalLng;
|
||
QPointF startG = start;
|
||
QPointF endG = end;
|
||
|
||
if (!m_internal->_location.contains(end)) // Clippen
|
||
{
|
||
clipped = true;
|
||
totalLng = length(endG - startG);
|
||
if (endG.x() < m_internal->_location.left())
|
||
{
|
||
double m = (endG.y() - startG.y()) / (endG.x() - startG.x());
|
||
double x = m_internal->_location.left();
|
||
endG.setX(x);
|
||
endG.setY(startG.y() + m * (x - startG.x()));
|
||
}
|
||
else if (endG.x() > m_internal->_location.right())
|
||
{
|
||
double m = (endG.y() - startG.y()) / (endG.x() - startG.x());
|
||
double x = m_internal->_location.right();
|
||
endG.setX(x);
|
||
endG.setY(startG.y() + m * (x - startG.x()));
|
||
}
|
||
if (endG.y() < m_internal->_location.bottom())
|
||
{
|
||
double m = (endG.x() - startG.x()) / (endG.y() - startG.y());
|
||
double y = m_internal->_location.bottom();
|
||
endG.setY(y);
|
||
endG.setX(startG.x() + m * (y - startG.y()));
|
||
}
|
||
else if (endG.y() > m_internal->_location.top())
|
||
{
|
||
double m = (endG.x() - startG.x()) / (endG.y() - startG.y());
|
||
double y = m_internal->_location.top();
|
||
endG.setY(y);
|
||
endG.setX(startG.x() + m * (y - startG.y()));
|
||
}
|
||
}
|
||
|
||
// Korrektursteigung aufgrund WGS84 Erdellipsoid, Lineare Approximation, Fehler bis 200 km: < 1.0m
|
||
double startRadius = radiusWGS84(startG.y());
|
||
double corrM = (radiusWGS84(endG.y()) - startRadius) / length(endG - startG);
|
||
|
||
Vec3dVector elevProfile;
|
||
getSamplePoints(start, end, elevProfile);
|
||
for (Vec3dVector::iterator itr = elevProfile.begin(); itr != elevProfile.end(); ++itr)
|
||
{
|
||
double dist = length(distance_in_meter_WGS84(start, itr->toPointF()));
|
||
double height = elevation(itr->x(), itr->y());
|
||
height += corrM * dist; // WGS84 Erdellipsoid Korrektur
|
||
height -= (startRadius * (1 - cos(atan(dist / startRadius)))); // Einbeziehen der Erdkruemmung
|
||
// alternativ als Naeherung: height -= ( dist * dist) / (startRadius * 2.0);
|
||
losProfile.push_back(LOSentry(dist, height));
|
||
}
|
||
|
||
if (clipped)
|
||
losProfile.push_back(LOSentry(totalLng, -NAN));
|
||
|
||
qCritical() << " losProfile.isze() = " << losProfile.size();
|
||
|
||
} // void ElevationMap::LOSProfile( const QPointF &start, const QPointF &end, LOSVector &losProfile)
|
||
|
||
|
||
|
||
} // SituationMap
|
||
|
||
|
||
|
||
/* -----------------------------------------------------------------------------------------------
|
||
** CVS-Id/Release ...: $Header$
|
||
** -----------------------------------------------------------------------------------------------
|
||
**
|
||
** Historie:
|
||
**
|
||
** $Log: not supported by cvs2svn $
|
||
** Revision 1.17.1.1 2008/06/16 06:05:37 joerg
|
||
** Transparenz und Bias einheitlich für alle Objekte implementiert
|
||
**
|
||
** Revision 1.17 2006/03/24 09:05:01 klaus
|
||
** Elevationprofile: Begrenzung auf Kartengröße aufgehoben ( Clipping)
|
||
**
|
||
** Revision 1.16 2005/11/09 13:44:36 klaus
|
||
** Umsetzung auf Linux;
|
||
**
|
||
** Revision 1.15 2005/04/15 13:05:45 joerg
|
||
** Resource Control
|
||
**
|
||
** Revision 1.14 2005/03/04 14:29:26 klaus
|
||
** Elevationmap mit CachedTileList versehen und automatischem Entladen von Tiles
|
||
**
|
||
** Revision 1.13 2005/03/03 12:37:53 klaus
|
||
** dynamischen Laden und Entladen der ElevationTiles für Reliefmap
|
||
**
|
||
** Revision 1.12 2005/03/01 12:36:59 klaus
|
||
** Reliefmap eingeführt, aus Elevationmap entfernt
|
||
**
|
||
** Revision 1.11 2005/02/25 15:32:56 klaus
|
||
** Ergänzung um Reliefkarte
|
||
**
|
||
** Revision 1.10 2005/02/10 14:15:19 klaus
|
||
** Bias in Basisklasse MapObject überführt
|
||
**
|
||
** Revision 1.9 2005/02/10 13:50:33 klaus
|
||
** Bias für SightFans eingeführt
|
||
**
|
||
** Revision 1.8 2005/02/09 16:45:26 klaus
|
||
** 1. Version der SightFans
|
||
**
|
||
** Revision 1.7 2005/01/26 12:26:23 klaus
|
||
** todos bearbeitet, Fehlerrückmeldnugen
|
||
**
|
||
** Revision 1.6 2004/12/14 10:29:48 klaus
|
||
** 7. Erweiterung für Höhenkarten
|
||
**
|
||
** Revision 1.5 2004/12/07 11:08:28 klaus
|
||
** 5. Erweiterung für Höhenkarten
|
||
**
|
||
** Revision 1.4 2004/12/03 16:25:14 klaus
|
||
** 4. Erweiterung für Höhenkarten
|
||
**
|
||
** Revision 1.3 2004/11/30 13:18:13 klaus
|
||
** 3. Erweiterung für Höhenkarten
|
||
**
|
||
** Revision 1.2 2004/11/23 14:36:40 klaus
|
||
** 2. Erweiterung für Höhenkarten
|
||
**
|
||
**
|
||
** ----------------------------------------------------------------------------------------------- */
|
||
|