diff --git a/src/lib/marble/geodata/data/GeoDataCoordinates.cpp b/src/lib/marble/geodata/data/GeoDataCoordinates.cpp index 854ed2965..64989f9ef 100644 --- a/src/lib/marble/geodata/data/GeoDataCoordinates.cpp +++ b/src/lib/marble/geodata/data/GeoDataCoordinates.cpp @@ -1,1184 +1,1197 @@ // // This file is part of the Marble Virtual Globe. // // This program is free software licensed under the GNU LGPL. You can // find a copy of this license in LICENSE.txt in the top directory of // the source code. // // Copyright 2004-2007 Torsten Rahn // Copyright 2007-2008 Inge Wallin // Copyright 2008 Patrick Spendrin // Copyright 2011 Friedrich W. H. Kossebau // Copyright 2011 Bernhard Beschow // Copyright 2015 Alejandro Garcia Montoro // #include "GeoDataCoordinates.h" #include "GeoDataCoordinates_p.h" #include "LonLatParser_p.h" #include #include #include #include "MarbleGlobal.h" #include "MarbleDebug.h" #include "MarbleMath.h" #include "Quaternion.h" namespace Marble { const qreal GeoDataCoordinatesPrivate::sm_semiMajorAxis = 6378137.0; const qreal GeoDataCoordinatesPrivate::sm_semiMinorAxis = 6356752.314; const qreal GeoDataCoordinatesPrivate::sm_eccentricitySquared = 6.69437999013e-03; const qreal GeoDataCoordinatesPrivate::sm_utmScaleFactor = 0.9996; GeoDataCoordinates::Notation GeoDataCoordinates::s_notation = GeoDataCoordinates::DMS; const GeoDataCoordinates GeoDataCoordinates::null = GeoDataCoordinates( 0, 0, 0 ); // don't use default constructor! GeoDataCoordinates::GeoDataCoordinates( qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit, int _detail ) : d( new GeoDataCoordinatesPrivate( _lon, _lat, _alt, unit, _detail ) ) { d->ref.ref(); } /* simply copy the d pointer * it will be replaced in the detach function instead */ GeoDataCoordinates::GeoDataCoordinates( const GeoDataCoordinates& other ) : d( other.d ) { d->ref.ref(); } /* simply copy null's d pointer * it will be replaced in the detach function */ GeoDataCoordinates::GeoDataCoordinates() : d( null.d ) { d->ref.ref(); } /* * only delete the private d pointer if the number of references is 0 * remember that all copies share the same d pointer! */ GeoDataCoordinates::~GeoDataCoordinates() { delete d->m_q; d->m_q = nullptr; if (!d->ref.deref()) delete d; #ifdef DEBUG_GEODATA // mDebug() << "delete coordinates"; #endif } bool GeoDataCoordinates::isValid() const { return d != null.d; } /* * if only one copy exists, return * else make a new private d pointer object and assign the values of the current * one to it * at the end, if the number of references thus reaches 0 delete it * this state shouldn't happen, but if it does, we have to clean up behind us. */ void GeoDataCoordinates::detach() { delete d->m_q; d->m_q = nullptr; if(d->ref.load() == 1) { return; } GeoDataCoordinatesPrivate *new_d = new GeoDataCoordinatesPrivate( *d ); if (!d->ref.deref()) { delete d; } d = new_d; d->ref.ref(); } /* * call detach() at the start of all non-static, non-const functions */ void GeoDataCoordinates::set( qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit ) { detach(); d->m_altitude = _alt; switch( unit ){ default: case Radian: d->m_lon = _lon; d->m_lat = _lat; break; case Degree: d->m_lon = _lon * DEG2RAD; d->m_lat = _lat * DEG2RAD; break; } } /* * call detach() at the start of all non-static, non-const functions */ void GeoDataCoordinates::setLongitude( qreal _lon, GeoDataCoordinates::Unit unit ) { detach(); switch( unit ){ default: case Radian: d->m_lon = _lon; break; case Degree: d->m_lon = _lon * DEG2RAD; break; } } /* * call detach() at the start of all non-static, non-const functions */ void GeoDataCoordinates::setLatitude( qreal _lat, GeoDataCoordinates::Unit unit ) { detach(); switch( unit ){ case Radian: d->m_lat = _lat; break; case Degree: d->m_lat = _lat * DEG2RAD; break; } } void GeoDataCoordinates::geoCoordinates( qreal& lon, qreal& lat, GeoDataCoordinates::Unit unit ) const { switch ( unit ) { default: case Radian: lon = d->m_lon; lat = d->m_lat; break; case Degree: lon = d->m_lon * RAD2DEG; lat = d->m_lat * RAD2DEG; break; } } void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat) const { lon = d->m_lon; lat = d->m_lat; } void GeoDataCoordinates::geoCoordinates( qreal& lon, qreal& lat, qreal& alt, GeoDataCoordinates::Unit unit ) const { geoCoordinates( lon, lat, unit ); alt = d->m_altitude; } void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat, qreal &alt) const { lon = d->m_lon; lat = d->m_lat; alt = d->m_altitude; } qreal GeoDataCoordinates::longitude( GeoDataCoordinates::Unit unit ) const { switch ( unit ) { default: case Radian: return d->m_lon; case Degree: return d->m_lon * RAD2DEG; } } qreal GeoDataCoordinates::longitude() const { return d->m_lon; } qreal GeoDataCoordinates::latitude( GeoDataCoordinates::Unit unit ) const { switch ( unit ) { default: case Radian: return d->m_lat; case Degree: return d->m_lat * RAD2DEG; } } qreal GeoDataCoordinates::latitude() const { return d->m_lat; } //static GeoDataCoordinates::Notation GeoDataCoordinates::defaultNotation() { return s_notation; } //static void GeoDataCoordinates::setDefaultNotation( GeoDataCoordinates::Notation notation ) { s_notation = notation; } //static qreal GeoDataCoordinates::normalizeLon( qreal lon, GeoDataCoordinates::Unit unit ) { qreal halfCircle; if ( unit == GeoDataCoordinates::Radian ) { halfCircle = M_PI; } else { halfCircle = 180; } if ( lon > halfCircle ) { int cycles = (int)( ( lon + halfCircle ) / ( 2 * halfCircle ) ); return lon - ( cycles * 2 * halfCircle ); } if ( lon < -halfCircle ) { int cycles = (int)( ( lon - halfCircle ) / ( 2 * halfCircle ) ); return lon - ( cycles * 2 * halfCircle ); } return lon; } //static qreal GeoDataCoordinates::normalizeLat( qreal lat, GeoDataCoordinates::Unit unit ) { qreal halfCircle; if ( unit == GeoDataCoordinates::Radian ) { halfCircle = M_PI; } else { halfCircle = 180; } if ( lat > ( halfCircle / 2.0 ) ) { int cycles = (int)( ( lat + halfCircle ) / ( 2 * halfCircle ) ); qreal temp; if( cycles == 0 ) { // pi/2 < lat < pi temp = halfCircle - lat; } else { temp = lat - ( cycles * 2 * halfCircle ); } if ( temp > ( halfCircle / 2.0 ) ) { return ( halfCircle - temp ); } if ( temp < ( -halfCircle / 2.0 ) ) { return ( -halfCircle - temp ); } return temp; } if ( lat < ( -halfCircle / 2.0 ) ) { int cycles = (int)( ( lat - halfCircle ) / ( 2 * halfCircle ) ); qreal temp; if( cycles == 0 ) { temp = -halfCircle - lat; } else { temp = lat - ( cycles * 2 * halfCircle ); } if ( temp > ( +halfCircle / 2.0 ) ) { return ( +halfCircle - temp ); } if ( temp < ( -halfCircle / 2.0 ) ) { return ( -halfCircle - temp ); } return temp; } return lat; } //static void GeoDataCoordinates::normalizeLonLat( qreal &lon, qreal &lat, GeoDataCoordinates::Unit unit ) { qreal halfCircle; if ( unit == GeoDataCoordinates::Radian ) { halfCircle = M_PI; } else { halfCircle = 180; } if ( lon > +halfCircle ) { int cycles = (int)( ( lon + halfCircle ) / ( 2 * halfCircle ) ); lon = lon - ( cycles * 2 * halfCircle ); } if ( lon < -halfCircle ) { int cycles = (int)( ( lon - halfCircle ) / ( 2 * halfCircle ) ); lon = lon - ( cycles * 2 * halfCircle ); } if ( lat > ( +halfCircle / 2.0 ) ) { int cycles = (int)( ( lat + halfCircle ) / ( 2 * halfCircle ) ); qreal temp; if( cycles == 0 ) { // pi/2 < lat < pi temp = halfCircle - lat; } else { temp = lat - ( cycles * 2 * halfCircle ); } if ( temp > ( +halfCircle / 2.0 ) ) { lat = +halfCircle - temp; } if ( temp < ( -halfCircle / 2.0 ) ) { lat = -halfCircle - temp; } lat = temp; if( lon > 0 ) { lon = -halfCircle + lon; } else { lon = halfCircle + lon; } } if ( lat < ( -halfCircle / 2.0 ) ) { int cycles = (int)( ( lat - halfCircle ) / ( 2 * halfCircle ) ); qreal temp; if( cycles == 0 ) { temp = -halfCircle - lat; } else { temp = lat - ( cycles * 2 * halfCircle ); } if ( temp > ( +halfCircle / 2.0 ) ) { lat = +halfCircle - temp; } if ( temp < ( -halfCircle / 2.0 ) ) { lat = -halfCircle - temp; } lat = temp; if( lon > 0 ) { lon = -halfCircle + lon; } else { lon = halfCircle + lon; } } return; } GeoDataCoordinates GeoDataCoordinates::fromString( const QString& string, bool& successful ) { LonLatParser parser; successful = parser.parse(string); if (successful) { return GeoDataCoordinates( parser.lon(), parser.lat(), 0, GeoDataCoordinates::Degree ); } else { return GeoDataCoordinates(); } } QString GeoDataCoordinates::toString() const { return GeoDataCoordinates::toString( s_notation ); } QString GeoDataCoordinates::toString( GeoDataCoordinates::Notation notation, int precision ) const { QString coordString; if( notation == GeoDataCoordinates::UTM ){ int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat); // Handle lack of UTM zone number in the poles const QString zoneString = (zoneNumber > 0) ? QString::number(zoneNumber) : QString(); QString bandString = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat); QString eastingString = QString::number(GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat), 'f', 2); QString northingString = QString::number(GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat), 'f', 2); return QString("%1%2 %3 m E, %4 m N").arg(zoneString, bandString, eastingString, northingString); } else{ coordString = lonToString( d->m_lon, notation, Radian, precision ) + QLatin1String(", ") + latToString( d->m_lat, notation, Radian, precision ); } return coordString; } QString GeoDataCoordinates::lonToString( qreal lon, GeoDataCoordinates::Notation notation, GeoDataCoordinates::Unit unit, int precision, char format ) { if( notation == GeoDataCoordinates::UTM ){ /** * @FIXME: UTM needs lon + lat to know zone number and easting * By now, this code returns the zone+easting of the point * (lon, equator), but this can differ a lot at different locations * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536 */ qreal lonRad = ( unit == Radian ) ? lon : lon * DEG2RAD; int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(lonRad, 0); // Handle lack of UTM zone number in the poles QString result = (zoneNumber > 0) ? QString::number(zoneNumber) : QString(); if(precision > 0){ QString eastingString = QString::number( GeoDataCoordinatesPrivate::lonLatToEasting(lonRad, 0), 'f', 2 ); result += QString(" %1 m E").arg(eastingString); } return result; } QString weString = ( lon < 0 ) ? tr("W") : tr("E"); QString lonString; qreal lonDegF = ( unit == Degree ) ? fabs( lon ) : fabs( (qreal)(lon) * RAD2DEG ); // Take care of -1 case precision = ( precision < 0 ) ? 5 : precision; if ( notation == DMS || notation == DM ) { int lonDeg = (int) lonDegF; qreal lonMinF = 60 * (lonDegF - lonDeg); int lonMin = (int) lonMinF; qreal lonSecF = 60 * (lonMinF - lonMin); int lonSec = (int) lonSecF; // Adjustment for fuzziness (like 49.999999999999999999999) if ( precision == 0 ) { lonDeg = qRound( lonDegF ); } else if ( precision <= 2 ) { lonMin = qRound( lonMinF ); } else if ( precision <= 4 && notation == DMS ) { lonSec = qRound( lonSecF ); } else { if ( notation == DMS ) { lonSec = lonSecF = qRound( lonSecF * qPow( 10, precision - 4 ) ) / qPow( 10, precision - 4 ); } else { lonMin = lonMinF = qRound( lonMinF * qPow( 10, precision - 2 ) ) / qPow( 10, precision - 2 ); } } if (lonSec > 59 && notation == DMS ) { lonSec = lonSecF = 0; lonMin = lonMinF = lonMinF + 1; } if (lonMin > 59) { lonMin = lonMinF = 0; lonDeg = lonDegF = lonDegF + 1; } // Evaluate the string lonString = QString::fromUtf8("%1\xc2\xb0").arg(lonDeg, 3, 10, QLatin1Char(' ')); if ( precision == 0 ) { return lonString + weString; } if ( notation == DMS || precision < 3 ) { lonString += QString(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0')); } if ( precision < 3 ) { return lonString + weString; } if ( notation == DMS ) { // Includes -1 case! if ( precision < 5 ) { lonString += QString(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0')); return lonString + weString; } lonString += QString(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0')); } else { lonString += QString(" %L3'").arg(lonMinF, precision + 1, 'f', precision - 2, QLatin1Char('0')); } } else if ( notation == GeoDataCoordinates::Decimal ) { lonString = QString::fromUtf8("%L1\xc2\xb0").arg(lonDegF, 4 + precision, format, precision, QLatin1Char(' ')); } else if ( notation == GeoDataCoordinates::Astro ) { if (lon < 0) { lon += ( unit == Degree ) ? 360 : 2 * M_PI; } qreal lonHourF = ( unit == Degree ) ? fabs( lon/15.0 ) : fabs( (qreal)(lon/15.0) * RAD2DEG ); int lonHour = (int) lonHourF; qreal lonMinF = 60 * (lonHourF - lonHour); int lonMin = (int) lonMinF; qreal lonSecF = 60 * (lonMinF - lonMin); int lonSec = (int) lonSecF; // Adjustment for fuzziness (like 49.999999999999999999999) if ( precision == 0 ) { lonHour = qRound( lonHourF ); } else if ( precision <= 2 ) { lonMin = qRound( lonMinF ); } else if ( precision <= 4 ) { lonSec = qRound( lonSecF ); } else { lonSec = lonSecF = qRound( lonSecF * qPow( 10, precision - 4 ) ) / qPow( 10, precision - 4 ); } if (lonSec > 59 ) { lonSec = lonSecF = 0; lonMin = lonMinF = lonMinF + 1; } if (lonMin > 59) { lonMin = lonMinF = 0; lonHour = lonHourF = lonHourF + 1; } // Evaluate the string lonString = QString::fromUtf8("%1h").arg(lonHour, 3, 10, QLatin1Char(' ')); if ( precision == 0 ) { return lonString; } lonString += QString(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0')); if ( precision < 3 ) { return lonString; } // Includes -1 case! if ( precision < 5 ) { lonString += QString(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0')); return lonString; } lonString += QString(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0')); return lonString; } return lonString + weString; } QString GeoDataCoordinates::lonToString() const { return GeoDataCoordinates::lonToString( d->m_lon , s_notation ); } QString GeoDataCoordinates::latToString( qreal lat, GeoDataCoordinates::Notation notation, GeoDataCoordinates::Unit unit, int precision, char format ) { if( notation == GeoDataCoordinates::UTM ){ /** * @FIXME: UTM needs lon + lat to know latitude band and northing * By now, this code returns the band+northing of the point * (meridian, lat), but this can differ a lot at different locations * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536 */ qreal latRad = ( unit == Radian ) ? lat : lat * DEG2RAD; QString result = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(0, latRad); if ( precision > 0 ){ QString northingString = QString::number( GeoDataCoordinatesPrivate::lonLatToNorthing(0, latRad), 'f', 2 ); result += QString(" %1 m N").arg(northingString); } return result; } QString pmString; QString nsString; if (notation == GeoDataCoordinates::Astro){ pmString = ( lat > 0 ) ? "+" : "-"; } else { nsString = ( lat > 0 ) ? tr("N") : tr("S"); } QString latString; qreal latDegF = ( unit == Degree ) ? fabs( lat ) : fabs( (qreal)(lat) * RAD2DEG ); // Take care of -1 case precision = ( precision < 0 ) ? 5 : precision; if ( notation == DMS || notation == DM || notation == Astro) { int latDeg = (int) latDegF; qreal latMinF = 60 * (latDegF - latDeg); int latMin = (int) latMinF; qreal latSecF = 60 * (latMinF - latMin); int latSec = (int) latSecF; // Adjustment for fuzziness (like 49.999999999999999999999) if ( precision == 0 ) { latDeg = qRound( latDegF ); } else if ( precision <= 2 ) { latMin = qRound( latMinF ); } else if ( precision <= 4 && notation == DMS ) { latSec = qRound( latSecF ); } else { if ( notation == DMS || notation == Astro ) { latSec = latSecF = qRound( latSecF * qPow( 10, precision - 4 ) ) / qPow( 10, precision - 4 ); } else { latMin = latMinF = qRound( latMinF * qPow( 10, precision - 2 ) ) / qPow( 10, precision - 2 ); } } if (latSec > 59 && ( notation == DMS || notation == Astro )) { latSecF = 0; latSec = latSecF; latMin = latMin + 1; } if (latMin > 59) { latMinF = 0; latMin = latMinF; latDeg = latDeg + 1; } // Evaluate the string latString = QString::fromUtf8("%1\xc2\xb0").arg(latDeg, 3, 10, QLatin1Char(' ')); if ( precision == 0 ) { return pmString + latString + nsString; } if ( notation == DMS || notation == Astro || precision < 3 ) { latString += QString(" %2\'").arg(latMin, 2, 10, QLatin1Char('0')); } if ( precision < 3 ) { return pmString + latString + nsString; } if ( notation == DMS || notation == Astro ) { // Includes -1 case! if ( precision < 5 ) { latString += QString(" %3\"").arg(latSec, 2, 'f', 0, QLatin1Char('0')); return latString + nsString; } latString += QString(" %L3\"").arg(latSecF, precision - 1, 'f', precision - 4, QLatin1Char('0')); } else { latString += QString(" %L3'").arg(latMinF, precision + 1, 'f', precision - 2, QLatin1Char('0')); } } else // notation = GeoDataCoordinates::Decimal { latString = QString::fromUtf8("%L1\xc2\xb0").arg(latDegF, 4 + precision, format, precision, QLatin1Char(' ')); } return pmString + latString + nsString; } QString GeoDataCoordinates::latToString() const { return GeoDataCoordinates::latToString( d->m_lat, s_notation ); } bool GeoDataCoordinates::operator==( const GeoDataCoordinates &rhs ) const { return *d == *rhs.d; } bool GeoDataCoordinates::operator!=( const GeoDataCoordinates &rhs ) const { return *d != *rhs.d; } void GeoDataCoordinates::setAltitude( const qreal altitude ) { detach(); d->m_altitude = altitude; } qreal GeoDataCoordinates::altitude() const { return d->m_altitude; } int GeoDataCoordinates::utmZone() const{ return GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat); } qreal GeoDataCoordinates::utmEasting() const{ return GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat); } QString GeoDataCoordinates::utmLatitudeBand() const{ return GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat); } qreal GeoDataCoordinates::utmNorthing() const{ return GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat); } quint8 GeoDataCoordinates::detail() const { return d->m_detail; } void GeoDataCoordinates::setDetail(quint8 detail) { detach(); d->m_detail = detail; } GeoDataCoordinates GeoDataCoordinates::rotateAround( const GeoDataCoordinates &axis, qreal angle, Unit unit ) const { const Quaternion quatAxis = Quaternion::fromEuler( -axis.latitude() , axis.longitude(), 0 ); const Quaternion rotationAmount = Quaternion::fromEuler( 0, 0, unit == Radian ? angle : angle * DEG2RAD ); const Quaternion resultAxis = quatAxis * rotationAmount * quatAxis.inverse(); return rotateAround(resultAxis); } GeoDataCoordinates GeoDataCoordinates::rotateAround(const Quaternion &rotAxis) const { Quaternion rotatedQuat = quaternion(); rotatedQuat.rotateAroundAxis(rotAxis); qreal rotatedLon, rotatedLat; rotatedQuat.getSpherical(rotatedLon, rotatedLat); return GeoDataCoordinates(rotatedLon, rotatedLat, altitude()); } qreal GeoDataCoordinates::bearing( const GeoDataCoordinates &other, Unit unit, BearingType type ) const { if ( type == FinalBearing ) { double const offset = unit == Degree ? 180.0 : M_PI; return offset + other.bearing( *this, unit, InitialBearing ); } qreal const delta = other.d->m_lon - d->m_lon; double const bearing = atan2( sin ( delta ) * cos ( other.d->m_lat ), cos( d->m_lat ) * sin( other.d->m_lat ) - sin( d->m_lat ) * cos( other.d->m_lat ) * cos ( delta ) ); return unit == Radian ? bearing : bearing * RAD2DEG; } GeoDataCoordinates GeoDataCoordinates::moveByBearing( qreal bearing, qreal distance ) const { qreal newLat = asin( sin(d->m_lat) * cos(distance) + cos(d->m_lat) * sin(distance) * cos(bearing) ); qreal newLon = d->m_lon + atan2( sin(bearing) * sin(distance) * cos(d->m_lat), cos(distance) - sin(d->m_lat) * sin(newLat) ); return GeoDataCoordinates( newLon, newLat ); } const Quaternion& GeoDataCoordinates::quaternion() const { if (d->m_q == nullptr) { d->m_q = new Quaternion(Quaternion::fromSpherical( d->m_lon , d->m_lat )); } return *d->m_q; } GeoDataCoordinates GeoDataCoordinates::interpolate( const GeoDataCoordinates &target, double t_ ) const { double const t = qBound( 0.0, t_, 1.0 ); Quaternion const quat = Quaternion::slerp( quaternion(), target.quaternion(), t ); qreal lon, lat; quat.getSpherical( lon, lat ); double const alt = (1.0-t) * d->m_altitude + t * target.d->m_altitude; return GeoDataCoordinates( lon, lat, alt ); } +GeoDataCoordinates GeoDataCoordinates::nlerp(const GeoDataCoordinates &target, double t) const +{ + qreal lon = 0.0; + qreal lat = 0.0; + + const Quaternion itpos = Quaternion::nlerp(quaternion(), target.quaternion(), t); + itpos.getSpherical(lon, lat); + + const qreal altitude = 0.5 * (d->m_altitude + target.altitude()); + + return GeoDataCoordinates(lon, lat, altitude); +} + GeoDataCoordinates GeoDataCoordinates::interpolate( const GeoDataCoordinates &before, const GeoDataCoordinates &target, const GeoDataCoordinates &after, double t_ ) const { double const t = qBound( 0.0, t_, 1.0 ); Quaternion const b1 = GeoDataCoordinatesPrivate::basePoint( before.quaternion(), quaternion(), target.quaternion() ); Quaternion const a2 = GeoDataCoordinatesPrivate::basePoint( quaternion(), target.quaternion(), after.quaternion() ); Quaternion const a = Quaternion::slerp( quaternion(), target.quaternion(), t ); Quaternion const b = Quaternion::slerp( b1, a2, t ); Quaternion c = Quaternion::slerp( a, b, 2 * t * (1.0-t) ); qreal lon, lat; c.getSpherical( lon, lat ); // @todo spline interpolation of altitude? double const alt = (1.0-t) * d->m_altitude + t * target.d->m_altitude; return GeoDataCoordinates( lon, lat, alt ); } bool GeoDataCoordinates::isPole( Pole pole ) const { // Evaluate the most likely case first: // The case where we haven't hit the pole and where our latitude is normalized // to the range of 90 deg S ... 90 deg N if ( fabs( (qreal) 2.0 * d->m_lat ) < M_PI ) { return false; } else { if ( fabs( (qreal) 2.0 * d->m_lat ) == M_PI ) { // Ok, we have hit a pole. Now let's check whether it's the one we've asked for: if ( pole == AnyPole ){ return true; } else { if ( pole == NorthPole && 2.0 * d->m_lat == +M_PI ) { return true; } if ( pole == SouthPole && 2.0 * d->m_lat == -M_PI ) { return true; } return false; } } // else { // FIXME: Should we just normalize latitude and longitude and be done? // While this might work well for persistent data it would create some // possible overhead for temporary data, so this needs careful thinking. mDebug() << "GeoDataCoordinates not normalized!"; // Only as a last resort we cover the unlikely case where // the latitude is not normalized to the range of // 90 deg S ... 90 deg N if ( fabs( (qreal) 2.0 * normalizeLat( d->m_lat ) ) < M_PI ) { return false; } else { // Ok, we have hit a pole. Now let's check whether it's the one we've asked for: if ( pole == AnyPole ){ return true; } else { if ( pole == NorthPole && 2.0 * d->m_lat == +M_PI ) { return true; } if ( pole == SouthPole && 2.0 * d->m_lat == -M_PI ) { return true; } return false; } } } } } qreal GeoDataCoordinates::sphericalDistanceTo(const GeoDataCoordinates &other) const { qreal lon2, lat2; other.geoCoordinates( lon2, lat2 ); // FIXME: Take the altitude into account! return distanceSphere(d->m_lon, d->m_lat, lon2, lat2); } GeoDataCoordinates& GeoDataCoordinates::operator=( const GeoDataCoordinates &other ) { qAtomicAssign(d, other.d); return *this; } void GeoDataCoordinates::pack( QDataStream& stream ) const { stream << d->m_lon; stream << d->m_lat; stream << d->m_altitude; } void GeoDataCoordinates::unpack( QDataStream& stream ) { // call detach even though it shouldn't be needed - one never knows detach(); stream >> d->m_lon; stream >> d->m_lat; stream >> d->m_altitude; } Quaternion GeoDataCoordinatesPrivate::basePoint( const Quaternion &q1, const Quaternion &q2, const Quaternion &q3 ) { Quaternion const a = (q2.inverse() * q3).log(); Quaternion const b = (q2.inverse() * q1).log(); return q2 * ((a+b)*-0.25).exp(); } qreal GeoDataCoordinatesPrivate::arcLengthOfMeridian( qreal phi ) { // Precalculate n qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis) / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis); // Precalculate alpha qreal const alpha = ( (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0) * (1.0 + (qPow (n, 2.0) / 4.0) + (qPow (n, 4.0) / 64.0) ); // Precalculate beta qreal const beta = (-3.0 * n / 2.0) + (9.0 * qPow (n, 3.0) / 16.0) + (-3.0 * qPow (n, 5.0) / 32.0); // Precalculate gamma qreal const gamma = (15.0 * qPow (n, 2.0) / 16.0) + (-15.0 * qPow (n, 4.0) / 32.0); // Precalculate delta qreal const delta = (-35.0 * qPow (n, 3.0) / 48.0) + (105.0 * qPow (n, 5.0) / 256.0); // Precalculate epsilon qreal const epsilon = (315.0 * qPow (n, 4.0) / 512.0); // Now calculate the sum of the series and return qreal const result = alpha * (phi + (beta * qSin (2.0 * phi)) + (gamma * qSin (4.0 * phi)) + (delta * qSin (6.0 * phi)) + (epsilon * qSin (8.0 * phi))); return result; } qreal GeoDataCoordinatesPrivate::centralMeridianUTM( qreal zone ) { return DEG2RAD*(-183.0 + (zone * 6.0)); } qreal GeoDataCoordinatesPrivate::footpointLatitude( qreal northing ) { // Precalculate n (Eq. 10.18) qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis) / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis); // Precalculate alpha (Eq. 10.22) // (Same as alpha in Eq. 10.17) qreal const alpha = ((GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0) * (1 + (qPow (n, 2.0) / 4) + (qPow (n, 4.0) / 64)); // Precalculate y (Eq. 10.23) qreal const y = northing / alpha; // Precalculate beta (Eq. 10.22) qreal const beta = (3.0 * n / 2.0) + (-27.0 * qPow (n, 3.0) / 32.0) + (269.0 * qPow (n, 5.0) / 512.0); // Precalculate gamma (Eq. 10.22) qreal const gamma = (21.0 * qPow (n, 2.0) / 16.0) + (-55.0 * qPow (n, 4.0) / 32.0); // Precalculate delta (Eq. 10.22) qreal const delta = (151.0 * qPow (n, 3.0) / 96.0) + (-417.0 * qPow (n, 5.0) / 128.0); // Precalculate epsilon (Eq. 10.22) qreal const epsilon = (1097.0 * qPow (n, 4.0) / 512.0); // Now calculate the sum of the series (Eq. 10.21) qreal const result = y + (beta * qSin (2.0 * y)) + (gamma * qSin (4.0 * y)) + (delta * qSin (6.0 * y)) + (epsilon * qSin (8.0 * y)); return result; } QPointF GeoDataCoordinatesPrivate::mapLonLatToXY( qreal lambda, qreal phi, qreal lambda0 ) { // Equation (10.15) // Precalculate second numerical eccentricity const qreal ep2 = (qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) - qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0)) / qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0); // Precalculate the square of nu, just an auxilar quantity const qreal nu2 = ep2 * qPow(qCos(phi), 2.0); // Precalculate the radius of curvature in prime vertical const qreal N = qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) / (GeoDataCoordinatesPrivate::sm_semiMinorAxis * qSqrt(1 + nu2)); // Precalculate the tangent of phi and its square const qreal t = qTan(phi); const qreal t2 = t * t; // Precalculate longitude difference const qreal l = lambda - lambda0; /* * Precalculate coefficients for l**n in the equations below * so a normal human being can read the expressions for easting * and northing * -- l**1 and l**2 have coefficients of 1.0 * * The actual used coefficients starts at coef[1], just to * follow the meaningful nomenclature in equation 10.15 * (coef[n] corresponds to qPow(l,n) factor) */ const qreal coef1 = 1; const qreal coef2 = 1; const qreal coef3 = 1.0 - t2 + nu2; const qreal coef4 = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2); const qreal coef5 = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - 58.0 * t2 * nu2; const qreal coef6 = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - 330.0 * t2 * nu2; const qreal coef7 = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2); const qreal coef8 = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2); // Calculate easting (x) const qreal easting = N * qCos(phi) * coef1 * l + (N / 6.0 * qPow(qCos(phi), 3.0) * coef3 * qPow(l, 3.0)) + (N / 120.0 * qPow(qCos(phi), 5.0) * coef5 * qPow(l, 5.0)) + (N / 5040.0 * qPow(qCos(phi), 7.0) * coef7 * qPow(l, 7.0)); // Calculate northing (y) const qreal northing = arcLengthOfMeridian (phi) + (t / 2.0 * N * qPow(qCos(phi), 2.0) * coef2 * qPow(l, 2.0)) + (t / 24.0 * N * qPow(qCos(phi), 4.0) * coef4 * qPow(l, 4.0)) + (t / 720.0 * N * qPow(qCos(phi), 6.0) * coef6 * qPow(l, 6.0)) + (t / 40320.0 * N * qPow(qCos(phi), 8.0) * coef8 * qPow(l, 8.0)); return QPointF(easting, northing); } int GeoDataCoordinatesPrivate::lonLatToZone( qreal lon, qreal lat ){ // Converts lon and lat to degrees qreal lonDeg = lon * RAD2DEG; qreal latDeg = lat * RAD2DEG; /* Round the value of the longitude when the distance to the nearest integer * is less than 0.0000001. This avoids fuzzy values such as -114.0000000001, which * can produce a misbehaviour when calculating the zone associated at the borders * of the zone intervals (for example, the interval [-114, -108[ is associated with * zone number 12; if the following rounding is not done, the value returned by * lonLatToZone(114,0) is 11 instead of 12, as the function actually receives * -114.0000000001, which is in the interval [-120,-114[, associated to zone 11 */ qreal precision = 0.0000001; if ( qAbs(lonDeg - qFloor(lonDeg)) < precision || qAbs(lonDeg - qCeil(lonDeg)) < precision ){ lonDeg = qRound(lonDeg); } // There is no numbering associated to the poles, special value 0 is returned. if ( latDeg < -80 || latDeg > 84 ) { return 0; } // Obtains the zone number handling all the so called "exceptions" // See problem: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions // See solution: http://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point // General int zoneNumber = static_cast( (lonDeg+180) / 6.0 ) + 1; // Southwest Norway if ( latDeg >= 56 && latDeg < 64 && lonDeg >= 3 && lonDeg < 12 ) { zoneNumber = 32; } // Svalbard if ( latDeg >= 72 && latDeg < 84 ) { if ( lonDeg >= 0 && lonDeg < 9 ) { zoneNumber = 31; } else if ( lonDeg >= 9 && lonDeg < 21 ) { zoneNumber = 33; } else if ( lonDeg >= 21 && lonDeg < 33 ) { zoneNumber = 35; } else if ( lonDeg >= 33 && lonDeg < 42 ) { zoneNumber = 37; } } return zoneNumber; } qreal GeoDataCoordinatesPrivate::lonLatToEasting( qreal lon, qreal lat ){ int zoneNumber = lonLatToZone( lon, lat ); if ( zoneNumber == 0 ){ qreal lonDeg = lon * RAD2DEG; zoneNumber = static_cast( (lonDeg+180) / 6.0 ) + 1; } QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY( lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber) ); // Adjust easting and northing for UTM system qreal easting = coordinates.x() * GeoDataCoordinatesPrivate::sm_utmScaleFactor + 500000.0; return easting; } QString GeoDataCoordinatesPrivate::lonLatToLatitudeBand( qreal lon, qreal lat ){ // Obtains the latitude bands handling all the so called "exceptions" // Converts lon and lat to degrees qreal lonDeg = lon * RAD2DEG; qreal latDeg = lat * RAD2DEG; // Regular latitude bands between 80 S and 80 N (that is, between 10 and 170 in the [0,180] interval) int bandLetterIndex = 24; //Avoids "may be used uninitialized" warning if ( latDeg < -80 ) { // South pole (A for zones 1-30, B for zones 31-60) bandLetterIndex = ( (lonDeg+180) < 6*31 ) ? 0 : 1; } else if ( latDeg >= -80 && latDeg <= 80 ) { // General (+2 because the general lettering starts in C) bandLetterIndex = static_cast( (latDeg+80.0) / 8.0 ) + 2; } else if ( latDeg >= 80 && latDeg < 84 ) { // Band X is extended 4 more degrees bandLetterIndex = 21; } else if ( latDeg >= 84 ) { // North pole (Y for zones 1-30, Z for zones 31-60) bandLetterIndex = ((lonDeg+180) < 6*31) ? 22 : 23; } return QString( "ABCDEFGHJKLMNPQRSTUVWXYZ?" ).at( bandLetterIndex ); } qreal GeoDataCoordinatesPrivate::lonLatToNorthing( qreal lon, qreal lat ){ int zoneNumber = lonLatToZone( lon, lat ); if ( zoneNumber == 0 ){ qreal lonDeg = lon * RAD2DEG; zoneNumber = static_cast( (lonDeg+180) / 6.0 ) + 1; } QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY( lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber) ); qreal northing = coordinates.y() * GeoDataCoordinatesPrivate::sm_utmScaleFactor; if ( northing < 0.0 ) { northing += 10000000.0; } return northing; } uint qHash(const GeoDataCoordinates &coordinates) { uint seed = ::qHash(coordinates.altitude()); seed = ::qHash(coordinates.latitude(), seed); return ::qHash(coordinates.longitude(), seed); } } diff --git a/src/lib/marble/geodata/data/GeoDataCoordinates.h b/src/lib/marble/geodata/data/GeoDataCoordinates.h index adce777b2..f8048e62d 100644 --- a/src/lib/marble/geodata/data/GeoDataCoordinates.h +++ b/src/lib/marble/geodata/data/GeoDataCoordinates.h @@ -1,438 +1,446 @@ // // This file is part of the Marble Virtual Globe. // // This program is free software licensed under the GNU LGPL. You can // find a copy of this license in LICENSE.txt in the top directory of // the source code. // // Copyright 2006-2007 Torsten Rahn // Copyright 2007-2008 Inge Wallin // Copyright 2008 Patrick Spendrin // Copyright 2015 Alejandro Garcia Montoro // #ifndef MARBLE_GEODATACOORDINATES_H #define MARBLE_GEODATACOORDINATES_H #include #include #include #include "geodata_export.h" #include "MarbleGlobal.h" class QString; namespace Marble { class GeoDataCoordinatesPrivate; class Quaternion; /** * @short A 3d point representation * * GeoDataCoordinates is the simple representation of a single three * dimensional point. It can be used all through out marble as the data type * for three dimensional objects. it comprises of a Quaternion for speed issues. * This class was introduced to reflect the difference between a simple 3d point * and the GeoDataGeometry object containing such a point. The latter is a * GeoDataPoint and is simply derived from GeoDataCoordinates. * @see GeoDataPoint */ class GEODATA_EXPORT GeoDataCoordinates { Q_DECLARE_TR_FUNCTIONS(GeoDataCoordinates) public: /** * @brief enum used constructor to specify the units used * * Internally we always use radian for mathematical convenience. * However the Marble's interfaces to the outside should default * to degrees. */ enum Unit{ Radian, Degree }; /** * @brief enum used to specify the notation / numerical system * * For degrees there exist two notations: * "Decimal" (base-10) and the "Sexagesimal DMS" (base-60) which is * traditionally used in cartography. Decimal notation * uses floating point numbers to specify parts of a degree. The * Sexagesimal DMS notation uses integer based * Degrees-(Arc)Minutes-(Arc)Seconds to describe parts of a degree. */ enum Notation{ Decimal, ///< "Decimal" notation (base-10) DMS, ///< "Sexagesimal DMS" notation (base-60) DM, ///< "Sexagesimal DM" notation (base-60) UTM, Astro /// < "RA and DEC" notation (used for astronomical sky coordinates) }; /** * @brief The BearingType enum specifies where to measure the bearing * along great circle arcs * * When traveling along a great circle arc defined by the two points * A and B, the bearing varies along the arc. The "InitialBearing" bearing * corresponds to the bearing value at A, the "FinalBearing" bearing to that * at B. */ enum BearingType { InitialBearing, FinalBearing }; // Type definitions using Vector = QVector; using PtrVector = QVector; GeoDataCoordinates( const GeoDataCoordinates& other ); /** * @brief constructs an invalid instance * * Constructs an invalid instance such that calling isValid() * on it will return @code false @endcode. */ GeoDataCoordinates(); /** * @brief create a geocoordinate from longitude and latitude * @param lon longitude * @param lat latitude * @param alt altitude in meters (default: 0) * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) * @param detail detail (default: 0) */ GeoDataCoordinates( qreal lon, qreal lat, qreal alt = 0, GeoDataCoordinates::Unit unit = GeoDataCoordinates::Radian, int detail = 0 ); virtual ~GeoDataCoordinates(); /** * @brief Returns @code true @endcode if the coordinate is valid, @code false @endcode otherwise. * @return whether the coordinate is valid * * A coordinate is valid, if at least one component has been set and the last * assignment was not an invalid GeoDataCoordinates object. */ bool isValid() const; /** * @brief (re)set the coordinates in a GeoDataCoordinates object * @param lon longitude * @param lat latitude * @param alt altitude in meters (default: 0) * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) */ void set( qreal lon, qreal lat, qreal alt = 0, GeoDataCoordinates::Unit unit = GeoDataCoordinates::Radian ); /** * @brief use this function to get the longitude and latitude with one * call - use the unit parameter to switch between Radian and DMS * @param lon longitude * @param lat latitude * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) */ void geoCoordinates(qreal& lon, qreal& lat, GeoDataCoordinates::Unit unit) const; void geoCoordinates(qreal& lon, qreal& lat) const; /** * @brief use this function to get the longitude, latitude and altitude * with one call - use the unit parameter to switch between Radian and DMS * @param lon longitude * @param lat latitude * @param alt altitude in meters * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) */ void geoCoordinates(qreal& lon, qreal& lat, qreal& alt, GeoDataCoordinates::Unit unit) const; void geoCoordinates(qreal& lon, qreal& lat, qreal& alt) const; /** * @brief set the longitude in a GeoDataCoordinates object * @param lon longitude * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) */ void setLongitude( qreal lon, GeoDataCoordinates::Unit unit = GeoDataCoordinates::Radian ); /** * @brief retrieves the longitude of the GeoDataCoordinates object * use the unit parameter to switch between Radian and DMS * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) * @return longitude */ qreal longitude(GeoDataCoordinates::Unit unit) const; qreal longitude() const; /** * @brief retrieves the latitude of the GeoDataCoordinates object * use the unit parameter to switch between Radian and DMS * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) * @return latitude */ qreal latitude( GeoDataCoordinates::Unit unit ) const; qreal latitude() const; /** * @brief set the longitude in a GeoDataCoordinates object * @param lat longitude * @param unit units that lon and lat get measured in * (default for Radian: north pole at pi/2, southpole at -pi/2) */ void setLatitude( qreal lat, GeoDataCoordinates::Unit unit = GeoDataCoordinates::Radian ); /** * @brief return the altitude of the Point in meters */ qreal altitude() const; /** * @brief set the altitude of the Point in meters * @param altitude altitude */ void setAltitude( const qreal altitude ); /** * @brief retrieves the UTM zone of the GeoDataCoordinates object. * If the point is located on one of the poles (latitude < 80S or * latitude > 84N) there is no UTM zone associated; in this case, * 0 is returned. * @return UTM zone. */ int utmZone() const; /** * @brief retrieves the UTM easting of the GeoDataCoordinates object, * in meters. * @return UTM easting */ qreal utmEasting() const; /** * @brief retrieves the UTM latitude band of the GeoDataCoordinates object * @return UTM latitude band */ QString utmLatitudeBand() const; /** * @brief retrieves the UTM northing of the GeoDataCoordinates object, * in meters * @return UTM northing */ qreal utmNorthing() const; /** * @brief return the detail flag * detail range: 0 for most important points, 5 for least important */ quint8 detail() const; /** * @brief set the detail flag * @param detail detail */ void setDetail(quint8 detail); /** * @brief Rotates one coordinate around another. * @param axis The coordinate that serves as a rotation axis * @param angle Rotation angle * @param unit Unit of the result * @return The coordinate rotated in anticlockwise direction */ GeoDataCoordinates rotateAround( const GeoDataCoordinates &axis, qreal angle, Unit unit = Radian ) const; GeoDataCoordinates rotateAround(const Quaternion &rotAxis) const; /** * @brief Returns the bearing (true bearing, the angle between the line defined * by this point and the other and the prime meridian) * @param other The second point that, together with this point, defines a line * @param unit Unit of the result * @param type Type of the bearing * @return The true bearing in the requested unit, not range normalized, * in clockwise direction, with the value 0 corresponding to north */ qreal bearing( const GeoDataCoordinates &other, Unit unit = Radian, BearingType type = InitialBearing ) const; /** * @brief Returns the coordinates of the resulting point after moving this point * according to the distance and bearing parameters * @param bearing the same as above * @param distance the distance on a unit sphere */ GeoDataCoordinates moveByBearing( qreal bearing, qreal distance ) const; /** * @brief return a Quaternion with the used coordinates */ const Quaternion &quaternion() const; /** * @brief slerp (spherical linear) interpolation between this coordinate and the given target coordinate * @param target Destination coordinate * @param t Fraction 0..1 to weight between this and target * @return Interpolated coordinate between this (t<=0.0) and target (t>=1.0) */ GeoDataCoordinates interpolate( const GeoDataCoordinates &target, double t ) const; + /** + * @brief nlerp (normalized linear interpolation) between this coordinates and the given target coordinates + * @param target Destination coordinates + * @param t Fraction 0..1 to weight between this and target + * @return Interpolated coordinate between this (t<=0.0) and target (t>=1.0) + */ + GeoDataCoordinates nlerp(const GeoDataCoordinates &target, double t) const; + /** * @brief squad (spherical and quadrangle) interpolation between b and c * @param before First base point * @param target Third base point (second interpolation point) * @param after Fourth base point * @param t Offset between b (t<=0) and c (t>=1) */ GeoDataCoordinates interpolate( const GeoDataCoordinates &before, const GeoDataCoordinates &target, const GeoDataCoordinates &after, double t ) const; /** * @brief return whether our coordinates represent a pole * This method can be used to check whether the coordinate equals one of * the poles. */ bool isPole( Pole = AnyPole ) const; /** * @brief This method calculates the shortest distance between two points on a sphere. * @brief See: http://en.wikipedia.org/wiki/Great-circle_distance */ qreal sphericalDistanceTo(const GeoDataCoordinates &other) const; /** * @brief return Notation of string representation */ static GeoDataCoordinates::Notation defaultNotation(); /** * @brief set the Notation of the string representation * @param notation Notation */ static void setDefaultNotation( GeoDataCoordinates::Notation notation ); /** * @brief normalize the longitude to always be -M_PI <= lon <= +M_PI (Radian). * @param lon longitude * @param unit unit of the result */ static qreal normalizeLon( qreal lon, GeoDataCoordinates::Unit = GeoDataCoordinates::Radian ); /** * @brief normalize latitude to always be in -M_PI / 2. <= lat <= +M_PI / 2 (Radian). * @param lat latitude * @param unit unit of the result */ static qreal normalizeLat( qreal lat, GeoDataCoordinates::Unit = GeoDataCoordinates::Radian ); /** * @brief normalize both longitude and latitude at the same time * This method normalizes both latitude and longitude, so that the * latitude and the longitude stay within the "usual" range. * NOTE: If the latitude exceeds M_PI/2 (+90.0 deg) or -M_PI/2 (-90.0 deg) * then this will be interpreted as a pole traversion where the point will * end up on the opposite side of the globe. Therefore the longitude will * change by M_PI (180 deg). * If you don't want this behaviour use both normalizeLat() and * normalizeLon() instead. * @param lon the longitude value * @param lat the latitude value * @param unit unit of the result */ static void normalizeLonLat( qreal &lon, qreal &lat, GeoDataCoordinates::Unit = GeoDataCoordinates::Radian ); /** * @brief try to parse the string into a coordinate pair * @param string the string * @param successful becomes true if the conversion succeeds * @return the geodatacoordinates */ static GeoDataCoordinates fromString( const QString &string, bool& successful ); /** * @brief return a string representation of the coordinate * this is a convenience function which uses the default notation */ QString toString() const; /** * @brief return a string with the notation given by notation * * @param notation set a notation different from the default one * @param precision set the number of digits below degrees. * The precision depends on the current notation: * For Decimal representation the precision is the number of * digits after the decimal point. * In DMS a precision of 1 or 2 shows the arc minutes; a precision * of 3 or 4 will show arc seconds. A precision beyond that will * increase the number of digits after the arc second decimal point. */ QString toString( GeoDataCoordinates::Notation notation, int precision = -1 ) const; static QString lonToString( qreal lon, GeoDataCoordinates::Notation notation, GeoDataCoordinates::Unit unit = Radian, int precision = -1, char format = 'f' ); /** * @brief return a string representation of longitude of the coordinate * convenience function that uses the default notation */ QString lonToString() const; static QString latToString( qreal lat, GeoDataCoordinates::Notation notation, GeoDataCoordinates::Unit unit = Radian, int precision = -1, char format = 'f' ); /** * @brief return a string representation of latitude of the coordinate * convenience function that uses the default notation */ QString latToString() const; bool operator==(const GeoDataCoordinates &other) const; bool operator!=(const GeoDataCoordinates &other) const; GeoDataCoordinates& operator=( const GeoDataCoordinates &other ); /** Serialize the contents of the feature to @p stream. */ void pack(QDataStream &stream) const; /** Unserialize the contents of the feature from @p stream. */ void unpack(QDataStream &stream); private: void detach(); GeoDataCoordinatesPrivate *d; static GeoDataCoordinates::Notation s_notation; static const GeoDataCoordinates null; }; GEODATA_EXPORT uint qHash(const GeoDataCoordinates& coordinates ); } Q_DECLARE_METATYPE( Marble::GeoDataCoordinates ) #endif diff --git a/src/lib/marble/geodata/data/GeoDataLineString.cpp b/src/lib/marble/geodata/data/GeoDataLineString.cpp index 1ae6aebef..5a2da6042 100644 --- a/src/lib/marble/geodata/data/GeoDataLineString.cpp +++ b/src/lib/marble/geodata/data/GeoDataLineString.cpp @@ -1,964 +1,954 @@ // // This file is part of the Marble Virtual Globe. // // This program is free software licensed under the GNU LGPL. You can // find a copy of this license in LICENSE.txt in the top directory of // the source code. // // Copyright 2008 Torsten Rahn // Copyright 2009 Patrick Spendrin // #include "GeoDataLineString.h" #include "GeoDataLineString_p.h" #include "GeoDataLinearRing.h" #include "GeoDataTypes.h" #include "Quaternion.h" #include "MarbleDebug.h" #include namespace Marble { GeoDataLineString::GeoDataLineString( TessellationFlags f ) : GeoDataGeometry( new GeoDataLineStringPrivate( f ) ) { // mDebug() << "1) GeoDataLineString created:" << p(); } GeoDataLineString::GeoDataLineString( GeoDataLineStringPrivate* priv ) : GeoDataGeometry( priv ) { // mDebug() << "2) GeoDataLineString created:" << p(); } GeoDataLineString::GeoDataLineString( const GeoDataGeometry & other ) : GeoDataGeometry( other ) { // mDebug() << "3) GeoDataLineString created:" << p(); } GeoDataLineString::~GeoDataLineString() { #ifdef DEBUG_GEODATA mDebug() << "delete Linestring"; #endif } const char *GeoDataLineString::nodeType() const { return GeoDataTypes::GeoDataLineStringType; } EnumGeometryId GeoDataLineString::geometryId() const { return GeoDataLineStringId; } GeoDataGeometry *GeoDataLineString::copy() const { return new GeoDataLineString(*this); } void GeoDataLineStringPrivate::interpolateDateLine( const GeoDataCoordinates & previousCoords, const GeoDataCoordinates & currentCoords, GeoDataCoordinates & previousAtDateLine, GeoDataCoordinates & currentAtDateLine, TessellationFlags f ) const { GeoDataCoordinates dateLineCoords; // mDebug() << Q_FUNC_INFO; if ( f.testFlag( RespectLatitudeCircle ) && previousCoords.latitude() == currentCoords.latitude() ) { dateLineCoords = currentCoords; } else { int recursionCounter = 0; dateLineCoords = findDateLine( previousCoords, currentCoords, recursionCounter ); } previousAtDateLine = dateLineCoords; currentAtDateLine = dateLineCoords; if ( previousCoords.longitude() < 0 ) { previousAtDateLine.setLongitude( -M_PI ); currentAtDateLine.setLongitude( +M_PI ); } else { previousAtDateLine.setLongitude( +M_PI ); currentAtDateLine.setLongitude( -M_PI ); } } GeoDataCoordinates GeoDataLineStringPrivate::findDateLine( const GeoDataCoordinates & previousCoords, const GeoDataCoordinates & currentCoords, int recursionCounter ) const { int currentSign = ( currentCoords.longitude() < 0.0 ) ? -1 : +1 ; int previousSign = ( previousCoords.longitude() < 0.0 ) ? -1 : +1 ; qreal longitudeDiff = fabs( previousSign * M_PI - previousCoords.longitude() ) + fabs( currentSign * M_PI - currentCoords.longitude() ); if ( longitudeDiff < 0.001 || recursionCounter == 100 ) { // mDebug() << "stopped at recursion" << recursionCounter << " and longitude difference " << longitudeDiff; return currentCoords; } ++recursionCounter; - qreal lon = 0.0; - qreal lat = 0.0; - - qreal altDiff = currentCoords.altitude() - previousCoords.altitude(); - - const Quaternion itpos = Quaternion::nlerp( previousCoords.quaternion(), currentCoords.quaternion(), 0.5 ); - itpos.getSpherical( lon, lat ); - - qreal altitude = previousCoords.altitude() + 0.5 * altDiff; - - GeoDataCoordinates interpolatedCoords( lon, lat, altitude ); + const GeoDataCoordinates interpolatedCoords = previousCoords.nlerp(currentCoords, 0.5); int interpolatedSign = ( interpolatedCoords.longitude() < 0.0 ) ? -1 : +1 ; /* mDebug() << "SRC" << previousCoords.toString(); mDebug() << "TAR" << currentCoords.toString(); mDebug() << "IPC" << interpolatedCoords.toString(); */ if ( interpolatedSign != currentSign ) { return findDateLine( interpolatedCoords, currentCoords, recursionCounter ); } return findDateLine( previousCoords, interpolatedCoords, recursionCounter ); } quint8 GeoDataLineStringPrivate::levelForResolution(qreal resolution) const { if (m_previousResolution == resolution) return m_level; m_previousResolution = resolution; if (resolution < 0.0000005) m_level = 17; else if (resolution < 0.0000010) m_level = 16; else if (resolution < 0.0000020) m_level = 15; else if (resolution < 0.0000040) m_level = 14; else if (resolution < 0.0000080) m_level = 13; else if (resolution < 0.0000160) m_level = 12; else if (resolution < 0.0000320) m_level = 11; else if (resolution < 0.0000640) m_level = 10; else if (resolution < 0.0001280) m_level = 9; else if (resolution < 0.0002560) m_level = 8; else if (resolution < 0.0005120) m_level = 7; else if (resolution < 0.0010240) m_level = 6; else if (resolution < 0.0020480) m_level = 5; else if (resolution < 0.0040960) m_level = 4; else if (resolution < 0.0081920) m_level = 3; else if (resolution < 0.0163840) m_level = 2; else m_level = 1; return m_level; } qreal GeoDataLineStringPrivate::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; } } void GeoDataLineStringPrivate::optimize (GeoDataLineString& lineString) const { QVector::iterator itCoords = lineString.begin(); QVector::const_iterator itEnd = lineString.constEnd(); if (lineString.size() < 2) return; // Calculate the least non-zero detail-level by checking the bounding box quint8 startLevel = levelForResolution( ( lineString.latLonAltBox().width() + lineString.latLonAltBox().height() ) / 2 ); quint8 currentLevel = startLevel; quint8 maxLevel = startLevel; GeoDataCoordinates currentCoords; lineString.first().setDetail(startLevel); // Iterate through the linestring to assign different detail levels to the nodes. // In general the first and last node should have the start level assigned as // a detail level. // Starting from the first node the algorithm picks those nodes which // have a distance from each other that is just above the resolution that is // associated with the start level (which we use as a "current level"). // Each of those nodes get the current level assigned as the detail level. // After iterating through the linestring we increment the current level value // and starting again with the first node we assign detail values in a similar way // to the remaining nodes which have no final detail level assigned yet. // We do as many iterations through the lineString as needed and bump up the // current level until all nodes have a non-zero detail level assigned. while ( currentLevel < 16 && currentLevel <= maxLevel + 1 ) { itCoords = lineString.begin(); currentCoords = *itCoords; ++itCoords; for( ; itCoords != itEnd; ++itCoords) { if (itCoords->detail() != 0 && itCoords->detail() < currentLevel) continue; if ( currentLevel == startLevel && (itCoords->longitude() == -M_PI || itCoords->longitude() == M_PI || itCoords->latitude() < -89 * DEG2RAD || itCoords->latitude() > 89 * DEG2RAD)) { itCoords->setDetail(startLevel); currentCoords = *itCoords; maxLevel = currentLevel; continue; } if (currentCoords.sphericalDistanceTo(*itCoords) < resolutionForLevel(currentLevel + 1)) { itCoords->setDetail(currentLevel + 1); } else { itCoords->setDetail(currentLevel); currentCoords = *itCoords; maxLevel = currentLevel; } } ++currentLevel; } lineString.last().setDetail(startLevel); } bool GeoDataLineString::isEmpty() const { Q_D(const GeoDataLineString); return d->m_vector.isEmpty(); } int GeoDataLineString::size() const { Q_D(const GeoDataLineString); return d->m_vector.size(); } GeoDataCoordinates& GeoDataLineString::at( int pos ) { detach(); Q_D(GeoDataLineString); d->m_dirtyRange = true; d->m_dirtyBox = true; return d->m_vector[pos]; } const GeoDataCoordinates& GeoDataLineString::at( int pos ) const { Q_D(const GeoDataLineString); return d->m_vector.at(pos); } GeoDataCoordinates& GeoDataLineString::operator[]( int pos ) { detach(); Q_D(GeoDataLineString); d->m_dirtyRange = true; d->m_dirtyBox = true; return d->m_vector[pos]; } GeoDataLineString GeoDataLineString::mid(int pos, int length) const { GeoDataLineString substring; auto d = substring.d_func(); d->m_vector = d_func()->m_vector.mid(pos, length); d->m_dirtyBox = true; d->m_dirtyRange = true; d->m_tessellationFlags = d_func()->m_tessellationFlags; d->m_extrude = d_func()->m_extrude; return substring; } const GeoDataCoordinates& GeoDataLineString::operator[]( int pos ) const { Q_D(const GeoDataLineString); return d->m_vector[pos]; } GeoDataCoordinates& GeoDataLineString::last() { detach(); Q_D(GeoDataLineString); d->m_dirtyRange = true; d->m_dirtyBox = true; return d->m_vector.last(); } GeoDataCoordinates& GeoDataLineString::first() { detach(); Q_D(GeoDataLineString); return d->m_vector.first(); } const GeoDataCoordinates& GeoDataLineString::last() const { Q_D(const GeoDataLineString); return d->m_vector.last(); } const GeoDataCoordinates& GeoDataLineString::first() const { Q_D(const GeoDataLineString); return d->m_vector.first(); } QVector::Iterator GeoDataLineString::begin() { detach(); Q_D(GeoDataLineString); return d->m_vector.begin(); } QVector::ConstIterator GeoDataLineString::begin() const { Q_D(const GeoDataLineString); return d->m_vector.constBegin(); } QVector::Iterator GeoDataLineString::end() { detach(); Q_D(GeoDataLineString); return d->m_vector.end(); } QVector::ConstIterator GeoDataLineString::end() const { Q_D(const GeoDataLineString); return d->m_vector.constEnd(); } QVector::ConstIterator GeoDataLineString::constBegin() const { Q_D(const GeoDataLineString); return d->m_vector.constBegin(); } QVector::ConstIterator GeoDataLineString::constEnd() const { Q_D(const GeoDataLineString); return d->m_vector.constEnd(); } void GeoDataLineString::insert( int index, const GeoDataCoordinates& value ) { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; d->m_vector.insert( index, value ); } void GeoDataLineString::append ( const GeoDataCoordinates& value ) { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; d->m_vector.append( value ); } void GeoDataLineString::reserve(int size) { Q_D(GeoDataLineString); d->m_vector.reserve(size); } void GeoDataLineString::append(const QVector& values) { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; #if QT_VERSION >= 0x050500 d->m_vector.append(values); #else d->m_vector.reserve(d->m_vector.size() + values.size()); for (const GeoDataCoordinates &coordinates: values) { d->m_vector.append(coordinates); } #endif } GeoDataLineString& GeoDataLineString::operator << ( const GeoDataCoordinates& value ) { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; d->m_vector.append( value ); return *this; } GeoDataLineString& GeoDataLineString::operator << ( const GeoDataLineString& value ) { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; QVector::const_iterator itCoords = value.constBegin(); QVector::const_iterator itEnd = value.constEnd(); d->m_vector.reserve(d->m_vector.size() + value.size()); for( ; itCoords != itEnd; ++itCoords ) { d->m_vector.append( *itCoords ); } return *this; } bool GeoDataLineString::operator==( const GeoDataLineString &other ) const { if ( !GeoDataGeometry::equals(other) || size() != other.size() || tessellate() != other.tessellate() ) { return false; } Q_D(const GeoDataLineString); const GeoDataLineStringPrivate* other_d = other.d_func(); QVector::const_iterator itCoords = d->m_vector.constBegin(); QVector::const_iterator otherItCoords = other_d->m_vector.constBegin(); QVector::const_iterator itEnd = d->m_vector.constEnd(); QVector::const_iterator otherItEnd = other_d->m_vector.constEnd(); for ( ; itCoords != itEnd && otherItCoords != otherItEnd; ++itCoords, ++otherItCoords ) { if ( *itCoords != *otherItCoords ) { return false; } } Q_ASSERT ( itCoords == itEnd && otherItCoords == otherItEnd ); return true; } bool GeoDataLineString::operator!=( const GeoDataLineString &other ) const { return !this->operator==(other); } void GeoDataLineString::clear() { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; d->m_vector.clear(); } bool GeoDataLineString::isClosed() const { return false; } bool GeoDataLineString::tessellate() const { Q_D(const GeoDataLineString); return d->m_tessellationFlags.testFlag(Tessellate); } void GeoDataLineString::setTessellate( bool tessellate ) { detach(); Q_D(GeoDataLineString); // According to the KML reference the tesselation of line strings in Google Earth // is generally done along great circles. However for subsequent points that share // the same latitude the latitude circles are followed. Our Tesselate and RespectLatitude // Flags provide this behaviour. For true polygons the latitude circles don't get considered. if (tessellate) { d->m_tessellationFlags |= (Tessellate | RespectLatitudeCircle); } else { d->m_tessellationFlags &= ~(Tessellate | RespectLatitudeCircle); } } TessellationFlags GeoDataLineString::tessellationFlags() const { Q_D(const GeoDataLineString); return d->m_tessellationFlags; } void GeoDataLineString::setTessellationFlags( TessellationFlags f ) { detach(); Q_D(GeoDataLineString); d->m_tessellationFlags = f; } void GeoDataLineString::reverse() { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; std::reverse(begin(), end()); } GeoDataLineString GeoDataLineString::toNormalized() const { Q_D(const GeoDataLineString); GeoDataLineString normalizedLineString; normalizedLineString.setTessellationFlags( tessellationFlags() ); qreal lon; qreal lat; // FIXME: Think about how we can avoid unnecessary copies // if the linestring stays the same. QVector::const_iterator end = d->m_vector.constEnd(); for( QVector::const_iterator itCoords = d->m_vector.constBegin(); itCoords != end; ++itCoords ) { itCoords->geoCoordinates( lon, lat ); qreal alt = itCoords->altitude(); GeoDataCoordinates::normalizeLonLat( lon, lat ); GeoDataCoordinates normalizedCoords( *itCoords ); normalizedCoords.set( lon, lat, alt ); normalizedLineString << normalizedCoords; } return normalizedLineString; } GeoDataLineString GeoDataLineString::toRangeCorrected() const { Q_D(const GeoDataLineString); if (d->m_dirtyRange) { delete d->m_rangeCorrected; if( isClosed() ) { d->m_rangeCorrected = new GeoDataLinearRing(toPoleCorrected()); } else { d->m_rangeCorrected = new GeoDataLineString(toPoleCorrected()); } d->m_dirtyRange = false; } return *d->m_rangeCorrected; } QVector GeoDataLineString::toDateLineCorrected() const { Q_D(const GeoDataLineString); QVector lineStrings; d->toDateLineCorrected(*this, lineStrings); return lineStrings; } GeoDataLineString GeoDataLineString::toPoleCorrected() const { Q_D(const GeoDataLineString); if( isClosed() ) { GeoDataLinearRing poleCorrected; d->toPoleCorrected(*this, poleCorrected); return poleCorrected; } else { GeoDataLineString poleCorrected; d->toPoleCorrected(*this, poleCorrected); return poleCorrected; } } void GeoDataLineStringPrivate::toPoleCorrected( const GeoDataLineString& q, GeoDataLineString& poleCorrected ) const { poleCorrected.setTessellationFlags( q.tessellationFlags() ); GeoDataCoordinates previousCoords; GeoDataCoordinates currentCoords; if ( q.isClosed() ) { if ( !( m_vector.first().isPole() ) && ( m_vector.last().isPole() ) ) { qreal firstLongitude = ( m_vector.first() ).longitude(); GeoDataCoordinates modifiedCoords( m_vector.last() ); modifiedCoords.setLongitude( firstLongitude ); poleCorrected << modifiedCoords; } } QVector::const_iterator itCoords = m_vector.constBegin(); QVector::const_iterator itEnd = m_vector.constEnd(); for( ; itCoords != itEnd; ++itCoords ) { currentCoords = *itCoords; if ( itCoords == m_vector.constBegin() ) { previousCoords = currentCoords; } if ( currentCoords.isPole() ) { if ( previousCoords.isPole() ) { continue; } else { qreal previousLongitude = previousCoords.longitude(); GeoDataCoordinates currentModifiedCoords( currentCoords ); currentModifiedCoords.setLongitude( previousLongitude ); poleCorrected << currentModifiedCoords; } } else { if ( previousCoords.isPole() ) { qreal currentLongitude = currentCoords.longitude(); GeoDataCoordinates previousModifiedCoords( previousCoords ); previousModifiedCoords.setLongitude( currentLongitude ); poleCorrected << previousModifiedCoords; poleCorrected << currentCoords; } else { // No poles at all. Nothing special to handle poleCorrected << currentCoords; } } previousCoords = currentCoords; } if ( q.isClosed() ) { if ( ( m_vector.first().isPole() ) && !( m_vector.last().isPole() ) ) { qreal lastLongitude = ( m_vector.last() ).longitude(); GeoDataCoordinates modifiedCoords( m_vector.first() ); modifiedCoords.setLongitude( lastLongitude ); poleCorrected << modifiedCoords; } } } void GeoDataLineStringPrivate::toDateLineCorrected( const GeoDataLineString & q, QVector & lineStrings ) const { const bool isClosed = q.isClosed(); const QVector::const_iterator itStartPoint = q.constBegin(); const QVector::const_iterator itEndPoint = q.constEnd(); QVector::const_iterator itPoint = itStartPoint; QVector::const_iterator itPreviousPoint = itPoint; TessellationFlags f = q.tessellationFlags(); GeoDataLineString * unfinishedLineString = nullptr; GeoDataLineString * dateLineCorrected = isClosed ? new GeoDataLinearRing( f ) : new GeoDataLineString( f ); qreal currentLon = 0.0; qreal previousLon = 0.0; int previousSign = 1; bool unfinished = false; for (; itPoint != itEndPoint; ++itPoint ) { currentLon = itPoint->longitude(); int currentSign = ( currentLon < 0.0 ) ? -1 : +1 ; if( itPoint == q.constBegin() ) { previousSign = currentSign; previousLon = currentLon; } // If we are crossing the date line ... if ( previousSign != currentSign && fabs(previousLon) + fabs(currentLon) > M_PI ) { unfinished = !unfinished; GeoDataCoordinates previousTemp; GeoDataCoordinates currentTemp; interpolateDateLine( *itPreviousPoint, *itPoint, previousTemp, currentTemp, q.tessellationFlags() ); *dateLineCorrected << previousTemp; if ( isClosed && unfinished ) { // If it's a linear ring and if it crossed the IDL only once then // store the current string inside the unfinishedLineString for later use ... unfinishedLineString = dateLineCorrected; // ... and start a new linear ring for now. dateLineCorrected = new GeoDataLinearRing( f ); } else { // Now it can only be a (finished) line string or a finished linear ring. // Store it in the vector if the size is not zero. if ( dateLineCorrected->size() > 0 ) { lineStrings << dateLineCorrected; } else { // Or delete it. delete dateLineCorrected; } // If it's a finished linear ring restore the "remembered" unfinished String if ( isClosed && !unfinished && unfinishedLineString ) { dateLineCorrected = unfinishedLineString; } else { // if it's a line string just create a new line string. dateLineCorrected = new GeoDataLineString( f ); } } *dateLineCorrected << currentTemp; *dateLineCorrected << *itPoint; } else { *dateLineCorrected << *itPoint; } previousSign = currentSign; previousLon = currentLon; itPreviousPoint = itPoint; } // If the line string doesn't cross the dateline an even number of times // then need to take care of the data stored in the unfinishedLineString if ( unfinished && unfinishedLineString && !unfinishedLineString->isEmpty() ) { *dateLineCorrected << *unfinishedLineString; delete unfinishedLineString; } lineStrings << dateLineCorrected; } const GeoDataLatLonAltBox& GeoDataLineString::latLonAltBox() const { Q_D(const GeoDataLineString); // GeoDataLatLonAltBox::fromLineString is very expensive // that's why we recreate it only if the m_dirtyBox // is TRUE. // DO NOT REMOVE THIS CONSTRUCT OR MARBLE WILL BE SLOW. if (d->m_dirtyBox) { d->m_latLonAltBox = GeoDataLatLonAltBox::fromLineString(*this); d->m_dirtyBox = false; } return d->m_latLonAltBox; } qreal GeoDataLineString::length( qreal planetRadius, int offset ) const { if( offset < 0 || offset >= size() ) { return 0; } Q_D(const GeoDataLineString); qreal length = 0.0; QVector const & vector = d->m_vector; int const start = qMax(offset+1, 1); int const end = d->m_vector.size(); for( int i=start; i::Iterator GeoDataLineString::erase ( const QVector::Iterator& pos ) { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; return d->m_vector.erase( pos ); } QVector::Iterator GeoDataLineString::erase ( const QVector::Iterator& begin, const QVector::Iterator& end ) { detach(); Q_D(GeoDataLineString); delete d->m_rangeCorrected; d->m_rangeCorrected = nullptr; d->m_dirtyRange = true; d->m_dirtyBox = true; return d->m_vector.erase( begin, end ); } void GeoDataLineString::remove ( int i ) { detach(); Q_D(GeoDataLineString); d->m_dirtyRange = true; d->m_dirtyBox = true; d->m_vector.remove( i ); } GeoDataLineString GeoDataLineString::optimized () const { Q_D(const GeoDataLineString); if( isClosed() ) { GeoDataLinearRing linearRing(*this); d->optimize(linearRing); return linearRing; } else { GeoDataLineString lineString(*this); d->optimize(lineString); return lineString; } } void GeoDataLineString::pack( QDataStream& stream ) const { Q_D(const GeoDataLineString); GeoDataGeometry::pack( stream ); stream << size(); stream << (qint32)(d->m_tessellationFlags); for( QVector::const_iterator iterator = d->m_vector.constBegin(); iterator != d->m_vector.constEnd(); ++iterator ) { mDebug() << "innerRing: size" << d->m_vector.size(); GeoDataCoordinates coord = ( *iterator ); coord.pack( stream ); } } void GeoDataLineString::unpack( QDataStream& stream ) { detach(); Q_D(GeoDataLineString); GeoDataGeometry::unpack( stream ); qint32 size; qint32 tessellationFlags; stream >> size; stream >> tessellationFlags; d->m_tessellationFlags = (TessellationFlags)(tessellationFlags); d->m_vector.reserve(d->m_vector.size() + size); for(qint32 i = 0; i < size; i++ ) { GeoDataCoordinates coord; coord.unpack( stream ); d->m_vector.append( coord ); } } } diff --git a/src/lib/marble/projections/AzimuthalProjection.cpp b/src/lib/marble/projections/AzimuthalProjection.cpp index ee94c9b5e..c97379b3f 100644 --- a/src/lib/marble/projections/AzimuthalProjection.cpp +++ b/src/lib/marble/projections/AzimuthalProjection.cpp @@ -1,681 +1,678 @@ // // This file is part of the Marble Virtual Globe. // // This program is free software licensed under the GNU LGPL. You can // find a copy of this license in LICENSE.txt in the top directory of // the source code. // // Copyright 2014 Gábor Péterffy // // Local #include "AzimuthalProjection.h" #include "AzimuthalProjection_p.h" #include "AbstractProjection_p.h" // Marble #include "GeoDataLinearRing.h" #include "GeoDataLineString.h" #include "GeoDataCoordinates.h" #include "GeoDataLatLonAltBox.h" #include "ViewportParams.h" #include namespace Marble { qreal AzimuthalProjection::maxValidLat() const { return +90.0 * DEG2RAD; } qreal AzimuthalProjection::minValidLat() const { return -90.0 * DEG2RAD; } bool AzimuthalProjection::isClippedToSphere() const { return true; } qreal AzimuthalProjection::clippingRadius() const { return 1; } bool AzimuthalProjection::screenCoordinates( const GeoDataLineString &lineString, const ViewportParams *viewport, QVector &polygons ) const { Q_D( const AzimuthalProjection ); // Compare bounding box size of the line string with the angularResolution // Immediately return if the latLonAltBox is smaller. if ( !viewport->resolves( lineString.latLonAltBox() ) ) { // mDebug() << "Object too small to be resolved"; return false; } d->lineStringToPolygon( lineString, viewport, polygons ); return true; } bool AzimuthalProjection::mapCoversViewport( const ViewportParams *viewport ) const { qint64 radius = viewport->radius() * viewport->currentProjection()->clippingRadius(); qint64 width = viewport->width(); qint64 height = viewport->height(); // This first test is a quick one that will catch all really big // radii and prevent overflow in the real test. if ( radius > width + height ) return true; // This is the real test. The 4 is because we are really // comparing to width/2 and height/2. if ( 4 * radius * radius >= width * width + height * height ) return true; return false; } GeoDataLatLonAltBox AzimuthalProjection::latLonAltBox( const QRect& screenRect, const ViewportParams *viewport ) const { // For the case where the whole viewport gets covered there is a // pretty dirty and generic detection algorithm: GeoDataLatLonAltBox latLonAltBox = AbstractProjection::latLonAltBox( screenRect, viewport ); // If the whole globe is visible we can easily calculate // analytically the lon-/lat- range. qreal pitch = GeoDataCoordinates::normalizeLat( viewport->planetAxis().pitch() ); if ( 2.0 * viewport->radius() <= viewport->height() && 2.0 * viewport->radius() <= viewport->width() ) { // Unless the planetaxis is in the screen plane the allowed longitude range // covers full -180 deg to +180 deg: if ( pitch > 0.0 && pitch < +M_PI ) { latLonAltBox.setWest( -M_PI ); latLonAltBox.setEast( +M_PI ); latLonAltBox.setNorth( +fabs( M_PI / 2.0 - fabs( pitch ) ) ); latLonAltBox.setSouth( -M_PI / 2.0 ); } if ( pitch < 0.0 && pitch > -M_PI ) { latLonAltBox.setWest( -M_PI ); latLonAltBox.setEast( +M_PI ); latLonAltBox.setNorth( +M_PI / 2.0 ); latLonAltBox.setSouth( -fabs( M_PI / 2.0 - fabs( pitch ) ) ); } // Last but not least we deal with the rare case where the // globe is fully visible and pitch = 0.0 or pitch = -M_PI or // pitch = +M_PI if ( pitch == 0.0 || pitch == -M_PI || pitch == +M_PI ) { qreal yaw = viewport->planetAxis().yaw(); latLonAltBox.setWest( GeoDataCoordinates::normalizeLon( yaw - M_PI / 2.0 ) ); latLonAltBox.setEast( GeoDataCoordinates::normalizeLon( yaw + M_PI / 2.0 ) ); latLonAltBox.setNorth( +M_PI / 2.0 ); latLonAltBox.setSouth( -M_PI / 2.0 ); } return latLonAltBox; } // Now we check whether maxLat (e.g. the north pole) gets displayed // inside the viewport to get more accurate values for east and west. // We need a point on the screen at maxLat that definitely gets displayed: qreal averageLongitude = ( latLonAltBox.west() + latLonAltBox.east() ) / 2.0; GeoDataCoordinates maxLatPoint( averageLongitude, maxLat(), 0.0, GeoDataCoordinates::Radian ); GeoDataCoordinates minLatPoint( averageLongitude, minLat(), 0.0, GeoDataCoordinates::Radian ); qreal dummyX, dummyY; // not needed bool dummyVal; if ( screenCoordinates( maxLatPoint, viewport, dummyX, dummyY, dummyVal ) || screenCoordinates( minLatPoint, viewport, dummyX, dummyY, dummyVal ) ) { latLonAltBox.setWest( -M_PI ); latLonAltBox.setEast( +M_PI ); } return latLonAltBox; } QPainterPath AzimuthalProjection::mapShape( const ViewportParams *viewport ) const { int radius = viewport->radius() * viewport->currentProjection()->clippingRadius(); int imgWidth = viewport->width(); int imgHeight = viewport->height(); QPainterPath fullRect; fullRect.addRect( 0 , 0 , imgWidth, imgHeight ); // If the globe covers the whole image, then the projected region represents // all of the image. // Otherwise the active region has got the shape of the visible globe. if ( !viewport->mapCoversViewport() ) { QPainterPath mapShape; mapShape.addEllipse( imgWidth / 2 - radius, imgHeight / 2 - radius, 2 * radius, 2 * radius ); return mapShape.intersected( fullRect ); } return fullRect; } AzimuthalProjection::AzimuthalProjection(AzimuthalProjectionPrivate * dd) : AbstractProjection( dd ) { } AzimuthalProjection::~AzimuthalProjection() { } void AzimuthalProjectionPrivate::tessellateLineSegment( const GeoDataCoordinates &aCoords, qreal ax, qreal ay, const GeoDataCoordinates &bCoords, qreal bx, qreal by, QVector &polygons, const ViewportParams *viewport, TessellationFlags f, bool allowLatePolygonCut ) const { // We take the manhattan length as a distance approximation // that can be too big by a factor of sqrt(2) qreal distance = fabs((bx - ax)) + fabs((by - ay)); #ifdef SAFE_DISTANCE // Interpolate additional nodes if the line segment that connects the // current or previous nodes might cross the viewport. // The latter can pretty safely be excluded for most projections if both points // are located on the same side relative to the viewport boundaries and if they are // located more than half the line segment distance away from the viewport. const qreal safeDistance = - 0.5 * distance; if ( !( bx < safeDistance && ax < safeDistance ) || !( by < safeDistance && ay < safeDistance ) || !( bx + safeDistance > viewport->width() && ax + safeDistance > viewport->width() ) || !( by + safeDistance > viewport->height() && ay + safeDistance > viewport->height() ) ) { #endif int maxTessellationFactor = viewport->radius() < 20000 ? 10 : 20; int const finalTessellationPrecision = qBound(2, viewport->radius()/200, maxTessellationFactor) * tessellationPrecision; // Let the line segment follow the spherical surface // if the distance between the previous point and the current point // on screen is too big if ( distance > finalTessellationPrecision ) { const int tessellatedNodes = qMin( distance / finalTessellationPrecision, maxTessellationNodes ); processTessellation( aCoords, bCoords, tessellatedNodes, polygons, viewport, f, allowLatePolygonCut); } else { crossHorizon( bCoords, polygons, viewport, allowLatePolygonCut ); } #ifdef SAFE_DISTANCE } #endif } void AzimuthalProjectionPrivate::processTessellation( const GeoDataCoordinates &previousCoords, const GeoDataCoordinates ¤tCoords, int tessellatedNodes, QVector &polygons, const ViewportParams *viewport, TessellationFlags f, bool allowLatePolygonCut ) const { const bool clampToGround = f.testFlag( FollowGround ); const bool followLatitudeCircle = f.testFlag( RespectLatitudeCircle ) && previousCoords.latitude() == currentCoords.latitude(); // Calculate steps for tessellation: lonDiff and altDiff qreal lonDiff = 0.0; if ( followLatitudeCircle ) { const int previousSign = previousCoords.longitude() > 0 ? 1 : -1; const int currentSign = currentCoords.longitude() > 0 ? 1 : -1; lonDiff = currentCoords.longitude() - previousCoords.longitude(); if ( previousSign != currentSign && fabs(previousCoords.longitude()) + fabs(currentCoords.longitude()) > M_PI ) { if ( previousSign > currentSign ) { // going eastwards -> lonDiff += 2 * M_PI ; } else { // going westwards -> lonDiff -= 2 * M_PI; } } } - const qreal altDiff = currentCoords.altitude() - previousCoords.altitude(); - // Create the tessellation nodes. GeoDataCoordinates previousTessellatedCoords = previousCoords; for ( int i = 1; i <= tessellatedNodes; ++i ) { const qreal t = (qreal)(i) / (qreal)( tessellatedNodes + 1 ); - // interpolate the altitude, too - const qreal altitude = clampToGround ? 0 : altDiff * t + previousCoords.altitude(); + GeoDataCoordinates currentTessellatedCoords; - qreal lon = 0.0; - qreal lat = 0.0; if ( followLatitudeCircle ) { // To tessellate along latitude circles use the // linear interpolation of the longitude. - lon = lonDiff * t + previousCoords.longitude(); - lat = previousTessellatedCoords.latitude(); + const qreal altDiff = currentCoords.altitude() - previousCoords.altitude(); + const qreal altitude = altDiff * t + previousCoords.altitude(); + const qreal lon = lonDiff * t + previousCoords.longitude(); + const qreal lat = previousTessellatedCoords.latitude(); + + currentTessellatedCoords = GeoDataCoordinates(lon, lat, altitude); } else { // To tessellate along great circles use the // normalized linear interpolation ("NLERP") for latitude and longitude. - const Quaternion itpos = Quaternion::nlerp( previousCoords.quaternion(), currentCoords.quaternion(), t ); - itpos. getSpherical( lon, lat ); + currentTessellatedCoords = previousCoords.nlerp(currentCoords, t); + } + + if (clampToGround) { + currentTessellatedCoords.setAltitude(0); } - const GeoDataCoordinates currentTessellatedCoords( lon, lat, altitude ); crossHorizon( currentTessellatedCoords, polygons, viewport, allowLatePolygonCut ); previousTessellatedCoords = currentTessellatedCoords; } // For the clampToGround case add the "current" coordinate after adding all other nodes. GeoDataCoordinates currentModifiedCoords( currentCoords ); if ( clampToGround ) { currentModifiedCoords.setAltitude( 0.0 ); } crossHorizon( currentModifiedCoords, polygons, viewport, allowLatePolygonCut ); } void AzimuthalProjectionPrivate::crossHorizon( const GeoDataCoordinates & bCoord, QVector &polygons, const ViewportParams *viewport, bool allowLatePolygonCut ) const { qreal x, y; bool globeHidesPoint; Q_Q( const AbstractProjection ); q->screenCoordinates( bCoord, viewport, x, y, globeHidesPoint ); if( !globeHidesPoint ) { *polygons.last() << QPointF( x, y ); } else { if ( allowLatePolygonCut && !polygons.last()->isEmpty() ) { QPolygonF *path = new QPolygonF; polygons.append( path ); } } } bool AzimuthalProjectionPrivate::lineStringToPolygon( const GeoDataLineString &lineString, const ViewportParams *viewport, QVector &polygons ) const { Q_Q( const AzimuthalProjection ); const TessellationFlags f = lineString.tessellationFlags(); bool const tessellate = lineString.tessellate(); const bool noFilter = f.testFlag(PreventNodeFiltering); qreal x = 0; qreal y = 0; bool globeHidesPoint = false; qreal previousX = -1.0; qreal previousY = -1.0; bool previousGlobeHidesPoint = false; qreal horizonX = -1.0; qreal horizonY = -1.0; QPolygonF * polygon = new QPolygonF; if (!tessellate) { polygon->reserve(lineString.size()); } polygons.append( polygon ); GeoDataLineString::ConstIterator itCoords = lineString.constBegin(); GeoDataLineString::ConstIterator itPreviousCoords = lineString.constBegin(); // Some projections display the earth in a way so that there is a // foreside and a backside. // The horizon is the line (usually a circle) which separates both // sides from each other and resembles the map shape. GeoDataCoordinates horizonCoords; // A horizon pair is a pair of two subsequent horizon crossings: // The first one describes the point where a line string disappears behind the horizon. // and where horizonPair is set to true. // The second one describes the point where the line string reappears. // In this case the two points are connected and horizonPair is set to false again. bool horizonPair = false; GeoDataCoordinates horizonDisappearCoords; // If the first horizon crossing in a line string describes the appearance of // a line string then we call it a "horizon orphan" and horizonOrphan is set to true. // In this case once the last horizon crossing in the line string is reached // it needs to be connected to the orphan. bool horizonOrphan = false; GeoDataCoordinates horizonOrphanCoords; GeoDataLineString::ConstIterator itBegin = lineString.constBegin(); GeoDataLineString::ConstIterator itEnd = lineString.constEnd(); bool processingLastNode = false; // We use a while loop to be able to cover linestrings as well as linear rings: // Linear rings require to tessellate the path from the last node to the first node // which isn't really convenient to achieve with a for loop ... const bool isLong = lineString.size() > 10; const int maximumDetail = levelForResolution(viewport->angularResolution()); // The first node of optimized linestrings has a non-zero detail value. const bool hasDetail = itBegin->detail() != 0; while ( itCoords != itEnd ) { // Optimization for line strings with a big amount of nodes bool skipNode = (hasDetail ? itCoords->detail() > maximumDetail : itCoords != itBegin && isLong && !processingLastNode && !viewport->resolves( *itPreviousCoords, *itCoords ) ); if ( !skipNode || noFilter) { q->screenCoordinates( *itCoords, viewport, x, y, globeHidesPoint ); // Initializing variables that store the values of the previous iteration if ( !processingLastNode && itCoords == itBegin ) { previousGlobeHidesPoint = globeHidesPoint; itPreviousCoords = itCoords; previousX = x; previousY = y; } // Check for the "horizon case" (which is present e.g. for the spherical projection const bool isAtHorizon = ( globeHidesPoint || previousGlobeHidesPoint ) && ( globeHidesPoint != previousGlobeHidesPoint ); if ( isAtHorizon ) { // Handle the "horizon case" horizonCoords = findHorizon( *itPreviousCoords, *itCoords, viewport, f ); if ( lineString.isClosed() ) { if ( horizonPair ) { horizonToPolygon( viewport, horizonDisappearCoords, horizonCoords, polygons.last() ); horizonPair = false; } else { if ( globeHidesPoint ) { horizonDisappearCoords = horizonCoords; horizonPair = true; } else { horizonOrphanCoords = horizonCoords; horizonOrphan = true; } } } q->screenCoordinates( horizonCoords, viewport, horizonX, horizonY ); // If the line appears on the visible half we need // to add an interpolated point at the horizon as the previous point. if ( previousGlobeHidesPoint ) { *polygons.last() << QPointF( horizonX, horizonY ); } } // This if-clause contains the section that tessellates the line // segments of a linestring. If you are about to learn how the code of // this class works you can safely ignore this section for a start. if ( lineString.tessellate() /* && ( isVisible || previousIsVisible ) */ ) { if ( !isAtHorizon ) { tessellateLineSegment( *itPreviousCoords, previousX, previousY, *itCoords, x, y, polygons, viewport, f, !lineString.isClosed() ); } else { // Connect the interpolated point at the horizon with the // current or previous point in the line. if ( previousGlobeHidesPoint ) { tessellateLineSegment( horizonCoords, horizonX, horizonY, *itCoords, x, y, polygons, viewport, f, !lineString.isClosed() ); } else { tessellateLineSegment( *itPreviousCoords, previousX, previousY, horizonCoords, horizonX, horizonY, polygons, viewport, f, !lineString.isClosed() ); } } } else { if ( !globeHidesPoint ) { *polygons.last() << QPointF( x, y ); } else { if ( !previousGlobeHidesPoint && isAtHorizon ) { *polygons.last() << QPointF( horizonX, horizonY ); } } } if ( globeHidesPoint ) { if ( !previousGlobeHidesPoint && !lineString.isClosed() ) { polygons.append( new QPolygonF ); } } previousGlobeHidesPoint = globeHidesPoint; itPreviousCoords = itCoords; previousX = x; previousY = y; } // Here we modify the condition to be able to process the // first node after the last node in a LinearRing. if ( processingLastNode ) { break; } ++itCoords; if ( itCoords == itEnd && lineString.isClosed() ) { itCoords = itBegin; processingLastNode = true; } } // In case of horizon crossings, make sure that we always get a // polygon closed correctly. if ( horizonOrphan && lineString.isClosed() ) { horizonToPolygon( viewport, horizonCoords, horizonOrphanCoords, polygons.last() ); } if ( polygons.last()->size() <= 1 ){ delete polygons.last(); polygons.pop_back(); // Clean up "unused" empty polygon instances } return polygons.isEmpty(); } void AzimuthalProjectionPrivate::horizonToPolygon( const ViewportParams *viewport, const GeoDataCoordinates & disappearCoords, const GeoDataCoordinates & reappearCoords, QPolygonF * polygon ) const { qreal x, y; const qreal imageHalfWidth = viewport->width() / 2; const qreal imageHalfHeight = viewport->height() / 2; bool dummyGlobeHidesPoint = false; Q_Q( const AzimuthalProjection ); // Calculate the angle of the position vectors of both coordinates q->screenCoordinates( disappearCoords, viewport, x, y, dummyGlobeHidesPoint ); qreal alpha = atan2( y - imageHalfHeight, x - imageHalfWidth ); q->screenCoordinates( reappearCoords, viewport, x, y, dummyGlobeHidesPoint ); qreal beta = atan2( y - imageHalfHeight, x - imageHalfWidth ); // Calculate the difference between both qreal diff = GeoDataCoordinates::normalizeLon( beta - alpha ); qreal sgndiff = diff < 0 ? -1 : 1; const qreal arcradius = q->clippingRadius() * viewport->radius(); const int itEnd = fabs(diff * RAD2DEG); // Create a polygon that resembles an arc between the two position vectors polygon->reserve(polygon->size() + itEnd); for ( int it = 1; it <= itEnd; ++it ) { const qreal angle = alpha + DEG2RAD * sgndiff * it; const qreal itx = imageHalfWidth + arcradius * cos( angle ); const qreal ity = imageHalfHeight + arcradius * sin( angle ); *polygon << QPointF( itx, ity ); } } GeoDataCoordinates AzimuthalProjectionPrivate::findHorizon( const GeoDataCoordinates & previousCoords, const GeoDataCoordinates & currentCoords, const ViewportParams *viewport, TessellationFlags f) const { bool currentHide = globeHidesPoint( currentCoords, viewport ) ; return doFindHorizon(previousCoords, currentCoords, viewport, f, currentHide, 0); } GeoDataCoordinates AzimuthalProjectionPrivate::doFindHorizon( const GeoDataCoordinates & previousCoords, const GeoDataCoordinates & currentCoords, const ViewportParams *viewport, TessellationFlags f, bool currentHide, int recursionCounter ) const { if ( recursionCounter > 20 ) { return currentHide ? previousCoords : currentCoords; } ++recursionCounter; bool followLatitudeCircle = false; // Calculate steps for tessellation: lonDiff and altDiff qreal lonDiff = 0.0; qreal previousLongitude = 0.0; qreal previousLatitude = 0.0; if ( f.testFlag( RespectLatitudeCircle ) ) { previousCoords.geoCoordinates( previousLongitude, previousLatitude ); qreal previousSign = previousLongitude > 0 ? 1 : -1; qreal currentLongitude = 0.0; qreal currentLatitude = 0.0; currentCoords.geoCoordinates( currentLongitude, currentLatitude ); qreal currentSign = currentLongitude > 0 ? 1 : -1; if ( previousLatitude == currentLatitude ) { followLatitudeCircle = true; lonDiff = currentLongitude - previousLongitude; if ( previousSign != currentSign && fabs(previousLongitude) + fabs(currentLongitude) > M_PI ) { if ( previousSign > currentSign ) { // going eastwards -> lonDiff += 2 * M_PI ; } else { // going westwards -> lonDiff -= 2 * M_PI; } } } else { // mDebug() << "Don't FollowLatitudeCircle"; } } - qreal lon = 0.0; - qreal lat = 0.0; - - qreal altDiff = currentCoords.altitude() - previousCoords.altitude(); + GeoDataCoordinates horizonCoords; if ( followLatitudeCircle ) { // To tessellate along latitude circles use the // linear interpolation of the longitude. - lon = lonDiff * 0.5 + previousLongitude; - lat = previousLatitude; + const qreal altDiff = currentCoords.altitude() - previousCoords.altitude(); + const qreal altitude = previousCoords.altitude() + 0.5 * altDiff; + const qreal lon = lonDiff * 0.5 + previousLongitude; + const qreal lat = previousLatitude; + + horizonCoords = GeoDataCoordinates(lon, lat, altitude); } else { // To tessellate along great circles use the // normalized linear interpolation ("NLERP") for latitude and longitude. - const Quaternion itpos = Quaternion::nlerp( previousCoords.quaternion(), currentCoords.quaternion(), 0.5 ); - itpos. getSpherical( lon, lat ); + horizonCoords = previousCoords.nlerp(currentCoords, 0.5); } - qreal altitude = previousCoords.altitude() + 0.5 * altDiff; - - GeoDataCoordinates horizonCoords( lon, lat, altitude ); - bool horizonHide = globeHidesPoint( horizonCoords, viewport ); if ( horizonHide != currentHide ) { return doFindHorizon(horizonCoords, currentCoords, viewport, f, currentHide, recursionCounter); } return doFindHorizon(previousCoords, horizonCoords, viewport, f, horizonHide, recursionCounter); } bool AzimuthalProjectionPrivate::globeHidesPoint( const GeoDataCoordinates &coordinates, const ViewportParams *viewport ) const { bool globeHidesPoint; qreal dummyX, dummyY; Q_Q( const AzimuthalProjection ); q->screenCoordinates(coordinates, viewport, dummyX, dummyY, globeHidesPoint); return globeHidesPoint; } } diff --git a/src/lib/marble/projections/CylindricalProjection.cpp b/src/lib/marble/projections/CylindricalProjection.cpp index 4d0079476..3eb62d8b4 100644 --- a/src/lib/marble/projections/CylindricalProjection.cpp +++ b/src/lib/marble/projections/CylindricalProjection.cpp @@ -1,481 +1,483 @@ // // This file is part of the Marble Virtual Globe. // // This program is free software licensed under the GNU LGPL. You can // find a copy of this license in LICENSE.txt in the top directory of // the source code. // // Copyright 2007 Inge Wallin // Copyright 2007-2012 Torsten Rahn // Copyright 2012 Cezar Mocan // // Local #include "CylindricalProjection.h" #include "CylindricalProjection_p.h" // Marble #include "GeoDataLinearRing.h" #include "GeoDataLineString.h" #include "GeoDataCoordinates.h" #include "GeoDataLatLonAltBox.h" #include "ViewportParams.h" #include // Maximum amount of nodes that are created automatically between actual nodes. static const int maxTessellationNodes = 200; namespace Marble { CylindricalProjection::CylindricalProjection() : AbstractProjection( new CylindricalProjectionPrivate( this ) ) { } CylindricalProjection::CylindricalProjection( CylindricalProjectionPrivate* dd ) : AbstractProjection( dd ) { } CylindricalProjection::~CylindricalProjection() { } CylindricalProjectionPrivate::CylindricalProjectionPrivate( CylindricalProjection * parent ) : AbstractProjectionPrivate( parent ), q_ptr( parent ) { } QPainterPath CylindricalProjection::mapShape( const ViewportParams *viewport ) const { // Convenience variables int width = viewport->width(); int height = viewport->height(); qreal yTop; qreal yBottom; qreal xDummy; // Get the top and bottom coordinates of the projected map. screenCoordinates( 0.0, maxLat(), viewport, xDummy, yTop ); screenCoordinates( 0.0, minLat(), viewport, xDummy, yBottom ); // Don't let the map area be outside the image if ( yTop < 0 ) yTop = 0; if ( yBottom > height ) yBottom = height; QPainterPath mapShape; mapShape.addRect( 0, yTop, width, yBottom - yTop ); return mapShape; } bool CylindricalProjection::screenCoordinates( const GeoDataLineString &lineString, const ViewportParams *viewport, QVector &polygons ) const { Q_D( const CylindricalProjection ); // Compare bounding box size of the line string with the angularResolution // Immediately return if the latLonAltBox is smaller. if ( !viewport->resolves( lineString.latLonAltBox() ) ) { // mDebug() << "Object too small to be resolved"; return false; } QVector subPolygons; d->lineStringToPolygon( lineString, viewport, subPolygons ); polygons << subPolygons; return polygons.isEmpty(); } int CylindricalProjectionPrivate::tessellateLineSegment( const GeoDataCoordinates &aCoords, qreal ax, qreal ay, const GeoDataCoordinates &bCoords, qreal bx, qreal by, QVector &polygons, const ViewportParams *viewport, TessellationFlags f, int mirrorCount, qreal repeatDistance) const { // We take the manhattan length as a distance approximation // that can be too big by a factor of sqrt(2) qreal distance = fabs((bx - ax)) + fabs((by - ay)); #ifdef SAFE_DISTANCE // Interpolate additional nodes if the line segment that connects the // current or previous nodes might cross the viewport. // The latter can pretty safely be excluded for most projections if both points // are located on the same side relative to the viewport boundaries and if they are // located more than half the line segment distance away from the viewport. const qreal safeDistance = - 0.5 * distance; if ( !( bx < safeDistance && ax < safeDistance ) || !( by < safeDistance && ay < safeDistance ) || !( bx + safeDistance > viewport->width() && ax + safeDistance > viewport->width() ) || !( by + safeDistance > viewport->height() && ay + safeDistance > viewport->height() ) ) { #endif int maxTessellationFactor = viewport->radius() < 20000 ? 10 : 20; int const finalTessellationPrecision = qBound(2, viewport->radius()/200, maxTessellationFactor) * tessellationPrecision; // Let the line segment follow the spherical surface // if the distance between the previous point and the current point // on screen is too big if ( distance > finalTessellationPrecision ) { const int tessellatedNodes = qMin( distance / finalTessellationPrecision, maxTessellationNodes ); mirrorCount = processTessellation( aCoords, bCoords, tessellatedNodes, polygons, viewport, f, mirrorCount, repeatDistance ); } else { mirrorCount = crossDateLine( aCoords, bCoords, bx, by, polygons, mirrorCount, repeatDistance ); } #ifdef SAFE_DISTANCE } #endif return mirrorCount; } int CylindricalProjectionPrivate::processTessellation( const GeoDataCoordinates &previousCoords, const GeoDataCoordinates ¤tCoords, int tessellatedNodes, QVector &polygons, const ViewportParams *viewport, TessellationFlags f, int mirrorCount, qreal repeatDistance) const { const bool clampToGround = f.testFlag( FollowGround ); const bool followLatitudeCircle = f.testFlag( RespectLatitudeCircle ) && previousCoords.latitude() == currentCoords.latitude(); // Calculate steps for tessellation: lonDiff and altDiff qreal lonDiff = 0.0; if ( followLatitudeCircle ) { const int previousSign = previousCoords.longitude() > 0 ? 1 : -1; const int currentSign = currentCoords.longitude() > 0 ? 1 : -1; lonDiff = currentCoords.longitude() - previousCoords.longitude(); if ( previousSign != currentSign && fabs(previousCoords.longitude()) + fabs(currentCoords.longitude()) > M_PI ) { if ( previousSign > currentSign ) { // going eastwards -> lonDiff += 2 * M_PI ; } else { // going westwards -> lonDiff -= 2 * M_PI; } } if ( fabs( lonDiff ) == 2 * M_PI ) { return mirrorCount; } } - const qreal altDiff = currentCoords.altitude() - previousCoords.altitude(); - // Create the tessellation nodes. GeoDataCoordinates previousTessellatedCoords = previousCoords; for ( int i = 1; i <= tessellatedNodes; ++i ) { const qreal t = (qreal)(i) / (qreal)( tessellatedNodes + 1 ); - // interpolate the altitude, too - const qreal altitude = clampToGround ? 0 : altDiff * t + previousCoords.altitude(); + GeoDataCoordinates currentTessellatedCoords; - qreal lon = 0.0; - qreal lat = 0.0; if ( followLatitudeCircle ) { // To tessellate along latitude circles use the // linear interpolation of the longitude. - lon = lonDiff * t + previousCoords.longitude(); - lat = previousTessellatedCoords.latitude(); + // interpolate the altitude, too + const qreal altDiff = currentCoords.altitude() - previousCoords.altitude(); + const qreal altitude = altDiff * t + previousCoords.altitude(); + const qreal lon = lonDiff * t + previousCoords.longitude(); + const qreal lat = previousTessellatedCoords.latitude(); + + currentTessellatedCoords = GeoDataCoordinates(lon, lat, altitude); } else { // To tessellate along great circles use the // normalized linear interpolation ("NLERP") for latitude and longitude. - const Quaternion itpos = Quaternion::nlerp( previousCoords.quaternion(), currentCoords.quaternion(), t ); - itpos. getSpherical( lon, lat ); + currentTessellatedCoords = previousCoords.nlerp(currentCoords, t); + } + + if (clampToGround) { + currentTessellatedCoords.setAltitude(0); } - const GeoDataCoordinates currentTessellatedCoords( lon, lat, altitude ); Q_Q(const CylindricalProjection); qreal bx, by; q->screenCoordinates( currentTessellatedCoords, viewport, bx, by ); mirrorCount = crossDateLine( previousTessellatedCoords, currentTessellatedCoords, bx, by, polygons, mirrorCount, repeatDistance ); previousTessellatedCoords = currentTessellatedCoords; } // For the clampToGround case add the "current" coordinate after adding all other nodes. GeoDataCoordinates currentModifiedCoords( currentCoords ); if ( clampToGround ) { currentModifiedCoords.setAltitude( 0.0 ); } Q_Q(const CylindricalProjection); qreal bx, by; q->screenCoordinates( currentModifiedCoords, viewport, bx, by ); mirrorCount = crossDateLine( previousTessellatedCoords, currentModifiedCoords, bx, by, polygons, mirrorCount, repeatDistance ); return mirrorCount; } int CylindricalProjectionPrivate::crossDateLine( const GeoDataCoordinates & aCoord, const GeoDataCoordinates & bCoord, qreal bx, qreal by, QVector &polygons, int mirrorCount, qreal repeatDistance ) { qreal aLon = aCoord.longitude(); qreal aSign = aLon > 0 ? 1 : -1; qreal bLon = bCoord.longitude(); qreal bSign = bLon > 0 ? 1 : -1; qreal delta = 0; if( aSign != bSign && fabs(aLon) + fabs(bLon) > M_PI ) { int sign = aSign > bSign ? 1 : -1; mirrorCount += sign; } delta = repeatDistance * mirrorCount; *polygons.last() << QPointF( bx + delta, by ); return mirrorCount; } bool CylindricalProjectionPrivate::lineStringToPolygon( const GeoDataLineString &lineString, const ViewportParams *viewport, QVector &polygons ) const { const TessellationFlags f = lineString.tessellationFlags(); bool const tessellate = lineString.tessellate(); const bool noFilter = f.testFlag(PreventNodeFiltering); qreal x = 0; qreal y = 0; qreal previousX = -1.0; qreal previousY = -1.0; int mirrorCount = 0; qreal distance = repeatDistance( viewport ); QPolygonF * polygon = new QPolygonF; if (!tessellate) { polygon->reserve(lineString.size()); } polygons.append( polygon ); GeoDataLineString::ConstIterator itCoords = lineString.constBegin(); GeoDataLineString::ConstIterator itPreviousCoords = lineString.constBegin(); GeoDataLineString::ConstIterator itBegin = lineString.constBegin(); GeoDataLineString::ConstIterator itEnd = lineString.constEnd(); bool processingLastNode = false; // We use a while loop to be able to cover linestrings as well as linear rings: // Linear rings require to tessellate the path from the last node to the first node // which isn't really convenient to achieve with a for loop ... const bool isLong = lineString.size() > 10; const int maximumDetail = levelForResolution(viewport->angularResolution()); // The first node of optimized linestrings has a non-zero detail value. const bool hasDetail = itBegin->detail() != 0; bool isStraight = lineString.latLonAltBox().height() == 0 || lineString.latLonAltBox().width() == 0; Q_Q( const CylindricalProjection ); bool const isClosed = lineString.isClosed(); while ( itCoords != itEnd ) { // Optimization for line strings with a big amount of nodes bool skipNode = (hasDetail ? itCoords->detail() > maximumDetail : isLong && !processingLastNode && itCoords != itBegin && !viewport->resolves( *itPreviousCoords, *itCoords ) ); if ( !skipNode || noFilter) { q->screenCoordinates( *itCoords, viewport, x, y ); // Initializing variables that store the values of the previous iteration if ( !processingLastNode && itCoords == itBegin ) { itPreviousCoords = itCoords; previousX = x; previousY = y; } // This if-clause contains the section that tessellates the line // segments of a linestring. If you are about to learn how the code of // this class works you can safely ignore this section for a start. if ( tessellate && !isStraight) { mirrorCount = tessellateLineSegment( *itPreviousCoords, previousX, previousY, *itCoords, x, y, polygons, viewport, f, mirrorCount, distance ); } else { // special case for polys which cross dateline but have no Tesselation Flag // the expected rendering is a screen coordinates straight line between // points, but in projections with repeatX things are not smooth mirrorCount = crossDateLine( *itPreviousCoords, *itCoords, x, y, polygons, mirrorCount, distance ); } itPreviousCoords = itCoords; previousX = x; previousY = y; } // Here we modify the condition to be able to process the // first node after the last node in a LinearRing. if ( processingLastNode ) { break; } ++itCoords; if (isClosed && itCoords == itEnd) { itCoords = itBegin; processingLastNode = true; } } // Closing e.g. in the Antarctica case. // This code makes the assumption that // - the first node is located at 180 E // - and the last node is located at 180 W // TODO: add a similar pattern in the crossDateLine() code. /* GeoDataLatLonAltBox box = lineString.latLonAltBox(); if( lineString.isClosed() && box.width() == 2*M_PI ) { QPolygonF *poly = polygons.last(); if( box.containsPole( NorthPole ) ) { qreal topMargin = 0.0; qreal dummy = 0.0; q_ptr->screenCoordinates(0.0, q_ptr->maxLat(), viewport, topMargin, dummy ); poly->push_back( QPointF( poly->last().x(), topMargin ) ); poly->push_back( QPointF( poly->first().x(), topMargin ) ); } else { qreal bottomMargin = 0.0; qreal dummy = 0.0; q_ptr->screenCoordinates(0.0, q_ptr->minLat(), viewport, bottomMargin, dummy ); poly->push_back( QPointF( poly->last().x(), bottomMargin ) ); poly->push_back( QPointF( poly->first().x(), bottomMargin ) ); } } */ repeatPolygons( viewport, polygons ); return polygons.isEmpty(); } void CylindricalProjectionPrivate::translatePolygons( const QVector &polygons, QVector &translatedPolygons, qreal xOffset ) { // mDebug() << "Translation: " << xOffset; translatedPolygons.reserve(polygons.size()); QVector::const_iterator itPolygon = polygons.constBegin(); QVector::const_iterator itEnd = polygons.constEnd(); for( ; itPolygon != itEnd; ++itPolygon ) { QPolygonF * polygon = new QPolygonF; *polygon = **itPolygon; polygon->translate( xOffset, 0 ); translatedPolygons.append( polygon ); } } void CylindricalProjectionPrivate::repeatPolygons( const ViewportParams *viewport, QVector &polygons ) const { Q_Q( const CylindricalProjection ); qreal xEast = 0; qreal xWest = 0; qreal y = 0; // Choose a latitude that is inside the viewport. const qreal centerLatitude = viewport->viewLatLonAltBox().center().latitude(); const GeoDataCoordinates westCoords(-M_PI, centerLatitude); const GeoDataCoordinates eastCoords(+M_PI, centerLatitude); q->screenCoordinates( westCoords, viewport, xWest, y ); q->screenCoordinates( eastCoords, viewport, xEast, y ); if ( xWest <= 0 && xEast >= viewport->width() - 1 ) { // mDebug() << "No repeats"; return; } const qreal repeatXInterval = xEast - xWest; const int repeatsLeft = (xWest > 0 ) ? (int)(xWest / repeatXInterval) + 1 : 0; const int repeatsRight = (xEast < viewport->width()) ? (int)((viewport->width() - xEast) / repeatXInterval) + 1 : 0; QVector repeatedPolygons; for (int it = repeatsLeft; it > 0; --it) { const qreal xOffset = -it * repeatXInterval; QVector translatedPolygons; translatePolygons( polygons, translatedPolygons, xOffset ); repeatedPolygons << translatedPolygons; } repeatedPolygons << polygons; for (int it = 1; it <= repeatsRight; ++it) { const qreal xOffset = +it * repeatXInterval; QVector translatedPolygons; translatePolygons( polygons, translatedPolygons, xOffset ); repeatedPolygons << translatedPolygons; } polygons = repeatedPolygons; // mDebug() << Q_FUNC_INFO << "Coordinates: " << xWest << xEast // << "Repeats: " << repeatsLeft << repeatsRight; } qreal CylindricalProjectionPrivate::repeatDistance( const ViewportParams *viewport ) const { // Choose a latitude that is inside the viewport. qreal centerLatitude = viewport->viewLatLonAltBox().center().latitude(); GeoDataCoordinates westCoords( -M_PI, centerLatitude ); GeoDataCoordinates eastCoords( +M_PI, centerLatitude ); qreal xWest, xEast, dummyY; Q_Q( const AbstractProjection ); q->screenCoordinates( westCoords, viewport, xWest, dummyY ); q->screenCoordinates( eastCoords, viewport, xEast, dummyY ); return xEast - xWest; } }