Commit 78b8313d authored by Dennis Nienhüser's avatar Dennis Nienhüser Committed by Friedrich W. H. Kossebau
Browse files

Use Douglas-Peucker algorithm for line/area simplification

parent 10611fac
......@@ -25,39 +25,34 @@ namespace Marble {
NodeReducer::NodeReducer(GeoDataDocument* document, int zoomLevel) :
BaseFilter(document),
m_resolution(resolutionForLevel(zoomLevel)),
m_removedNodes(0),
m_remainingNodes(0)
{
/*
* @todo FIXME distance based reduction is too simple for polygons, needs sth like Douglas-Peucker
*/
foreach (GeoDataPlacemark* placemark, placemarks()) {
GeoDataGeometry const * const geometry = placemark->geometry();
if(geometry->nodeType() == GeoDataTypes::GeoDataLineStringType) {
GeoDataLineString const * prevLine = static_cast<GeoDataLineString const *>(geometry);
GeoDataLineString* reducedLine = new GeoDataLineString;
reduce(prevLine, placemark->osmData(), reducedLine);
reduce(*prevLine, placemark, reducedLine, zoomLevel);
placemark->setGeometry(reducedLine);
} else if (zoomLevel < 17) {
if(geometry->nodeType() == GeoDataTypes::GeoDataLinearRingType) {
GeoDataLinearRing const * prevRing = static_cast<GeoDataLinearRing const *>(geometry);
GeoDataLinearRing* reducedRing = new GeoDataLinearRing;
reduce(prevRing, placemark->osmData(), reducedRing);
reduce(*prevRing, placemark, reducedRing, zoomLevel);
placemark->setGeometry(reducedRing);
} else if(geometry->nodeType() == GeoDataTypes::GeoDataPolygonType) {
GeoDataPolygon* reducedPolygon = new GeoDataPolygon;
GeoDataPolygon const * prevPolygon = static_cast<GeoDataPolygon const *>(geometry);
GeoDataLinearRing const * prevRing = &(prevPolygon->outerBoundary());
GeoDataLinearRing reducedRing;
reduce(prevRing, placemark->osmData(), &reducedRing);
reduce(*prevRing, placemark, &reducedRing, zoomLevel);
reducedPolygon->setOuterBoundary(reducedRing);
QVector<GeoDataLinearRing> const & innerBoundaries = prevPolygon->innerBoundaries();
for(int i = 0; i < innerBoundaries.size(); i++) {
prevRing = &innerBoundaries[i];
GeoDataLinearRing reducedInnerRing;
reduce(prevRing, placemark->osmData(), &reducedInnerRing);
reduce(*prevRing, placemark, &reducedInnerRing, zoomLevel);
reducedPolygon->appendInnerBoundary(reducedInnerRing);
}
placemark->setGeometry(reducedPolygon);
......@@ -71,69 +66,51 @@ qint64 NodeReducer::remainingNodes() const
return m_remainingNodes;
}
qint64 NodeReducer::removedNodes() const
qreal NodeReducer::epsilonForString(int detailLevel) const
{
return m_removedNodes;
int const factor = 1 << (qAbs(detailLevel-12));
return detailLevel < 12 ? 60.0 * factor : 60.0 / factor;
}
qreal NodeReducer::resolutionForLevel(int level) {
switch (level) {
case 0:
return 0.0655360;
break;
case 1:
return 0.0327680;
break;
case 2:
return 0.0163840;
break;
case 3:
return 0.0081920;
break;
case 4:
return 0.0040960;
break;
case 5:
return 0.0020480;
break;
case 6:
return 0.0010240;
break;
case 7:
return 0.0005120;
break;
case 8:
return 0.0002560;
break;
case 9:
return 0.0001280;
break;
case 10:
return 0.0000640;
break;
case 11:
return 0.0000320;
break;
case 12:
return 0.0000160;
break;
case 13:
return 0.0000080;
break;
case 14:
return 0.0000040;
break;
case 15:
return 0.0000020;
break;
case 16:
return 0.0000010;
break;
default:
case 17:
return 0.0000005;
break;
qreal NodeReducer::epsilonForArea(int detailLevel) const
{
int const factor = 1 << (qAbs(detailLevel-12));
return detailLevel < 12 ? 90.0 * factor : 90.0 / factor;
}
qreal NodeReducer::perpendicularDistance(const GeoDataCoordinates &a, const GeoDataCoordinates &b, const GeoDataCoordinates &c) const
{
qreal ret;
qreal const y0 = a.latitude();
qreal const x0 = a.longitude();
qreal const y1 = b.latitude();
qreal const x1 = b.longitude();
qreal const y2 = c.latitude();
qreal const x2 = c.longitude();
qreal const y01 = x0 - x1;
qreal const x01 = y0 - y1;
qreal const y10 = x1 - x0;
qreal const x10 = y1 - y0;
qreal const y21 = x2 - x1;
qreal const x21 = y2 - y1;
qreal const len = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);
qreal const t = (x01 * x21 + y01 * y21) / len;
if ( t < 0.0 ) {
ret = EARTH_RADIUS * distanceSphere(a, b);
} else if ( t > 1.0 ) {
ret = EARTH_RADIUS * distanceSphere(a, c);
} else {
qreal const nom = qAbs( x21 * y10 - x10 * y21 );
qreal const den = sqrt( x21 * x21 + y21 * y21 );
ret = EARTH_RADIUS * nom / den;
}
return ret;
}
qint64 NodeReducer::removedNodes() const
{
return m_removedNodes;
}
}
......@@ -14,6 +14,7 @@
#include "BaseFilter.h"
#include "MarbleMath.h"
#include "OsmPlacemarkData.h"
#include "VectorClipper.h"
namespace Marble {
......@@ -25,39 +26,78 @@ public:
qint64 remainingNodes() const;
private:
qreal epsilonForString(int detailLevel) const;
qreal epsilonForArea(int detailLevel) const;
qreal perpendicularDistance(const GeoDataCoordinates &a, const GeoDataCoordinates &b, const GeoDataCoordinates &c) const;
template<class T>
void reduce(T const * lineString, const OsmPlacemarkData &osmData, T* reducedLine)
void reduce(T const & lineString, const GeoDataPlacemark* placemark, T* reducedLine, int tileLevel)
{
qint64 const prevSize = lineString->size();
if (prevSize < 2) {
m_remainingNodes += prevSize;
return;
}
auto iter = lineString->begin();
GeoDataCoordinates currentCoords = *iter;
reducedLine->append(*iter);
++iter;
for (auto const end = lineString->end() - 1; iter != end; ++iter) {
if (!osmData.nodeReference(currentCoords).isEmpty()) {
continue; // do not remove nodes with tags
}
if (distanceSphere( currentCoords, *iter ) >= m_resolution) {
currentCoords = *iter;
reducedLine->append(*iter);
}
}
reducedLine->append(*iter);
bool const isArea = lineString.isClosed() && VectorClipper::canBeArea(placemark->visualCategory());
qreal const epsilon = isArea ? epsilonForArea(tileLevel) : epsilonForString(tileLevel);
*reducedLine = douglasPeucker(lineString, placemark->osmData(), epsilon);
qint64 prevSize = lineString.size();
qint64 reducedSize = reducedLine->size();
m_removedNodes += (prevSize - reducedSize);
m_remainingNodes += reducedSize;
//qDebug()<<"Nodes reduced "<<(prevSize - reducedSize)<<endl;
}
static qreal resolutionForLevel(int level);
template<class T>
T extract(T const & lineString, int start, int end) const
{
T result;
for (int i=start; i<=end; ++i) {
result << lineString[i];
}
return result;
}
template<class T>
T merge(T const & a, T const &b) const
{
T result = a;
for (int i=1, n=b.size(); i<n; ++i) {
result << b[i];
}
return result;
}
template<class T>
T douglasPeucker(T const & lineString, const OsmPlacemarkData &osmData, qreal epsilon) const
{
if (lineString.size() < 3) {
return lineString;
}
// @todo Keep nodes with tags
// if (!osmData.nodeReference(currentCoords).isEmpty()) {
// continue; // do not remove nodes with tags
// }
double maxDistance = 0.0;
int index = 0;
int const end = lineString.size()-1;
for (int i = 1; i<end; ++i) {
double const distance = perpendicularDistance(lineString[i], lineString[0], lineString[end]);
if (distance > maxDistance) {
index = i;
maxDistance = distance;
}
}
if (maxDistance >= epsilon) {
T const left = douglasPeucker(extract(lineString, 0, index), osmData, epsilon);
T const right = douglasPeucker(extract(lineString, index, end), osmData, epsilon);
return merge(left, right);
}
T result;
result << lineString[0];
result << lineString[end];
return result;
}
qreal m_resolution;
qint64 m_removedNodes;
qint64 m_remainingNodes;
};
......
......@@ -289,7 +289,7 @@ ClipperLib::Path VectorClipper::clipPath(const GeoDataLatLonBox &box) const
return path;
}
bool VectorClipper::canBeArea(GeoDataPlacemark::GeoDataVisualCategory visualCategory) const
bool VectorClipper::canBeArea(GeoDataPlacemark::GeoDataVisualCategory visualCategory)
{
if (visualCategory >= GeoDataPlacemark::HighwaySteps && visualCategory <= GeoDataPlacemark::HighwayMotorway) {
return false;
......
......@@ -34,12 +34,12 @@ public:
GeoDataDocument* clipTo(const GeoDataLatLonBox &box, bool filterSmallAreas);
GeoDataDocument* clipTo(unsigned int zoomLevel, unsigned int tileX, unsigned int tileY);
static bool canBeArea(GeoDataPlacemark::GeoDataVisualCategory visualCategory);
private:
GeoDataDocument* clipToBaseClipper(const GeoDataLatLonBox &box);
QVector<GeoDataPlacemark*> potentialIntersections(const GeoDataLatLonBox &box) const;
ClipperLib::Path clipPath(const GeoDataLatLonBox &box) const;
bool canBeArea(GeoDataPlacemark::GeoDataVisualCategory visualCategory) const;
qreal area(const GeoDataLinearRing &ring);
template<class T>
......
......@@ -273,9 +273,9 @@ int main(int argc, char *argv[])
originalWays = concatenator.originalWays();
mergedWays = concatenator.mergedWays();
}
NodeReducer nodeReducer(tile1.data(), zoomLevel);
GeoDocPtr tile2 = GeoDocPtr(loader.clip(zoomLevel, tileId.x(), tileId.y()));
GeoDocPtr combined = GeoDocPtr(mergeDocuments(tile1.data(), tile2.data()));
NodeReducer nodeReducer(combined.data(), zoomLevel);
if (boundaryTiles.contains(iter.key())) {
writeBoundaryTile(tile1.data(), region, parser, tileId.x(), tileId.y(), zoomLevel);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment