File indexing completed on 2024-05-05 03:49:51

0001 // SPDX-License-Identifier: LGPL-2.1-or-later
0002 //
0003 // SPDX-FileCopyrightText: 2014 Torsten Rahn <rahn@kde.org>
0004 //
0005 
0006 // Local
0007 #include "AzimuthalEquidistantProjection.h"
0008 #include "AbstractProjection_p.h"
0009 
0010 #include "MarbleDebug.h"
0011 
0012 // Marble
0013 #include "ViewportParams.h"
0014 #include "GeoDataPoint.h"
0015 #include "GeoDataLineString.h"
0016 #include "GeoDataCoordinates.h"
0017 #include "MarbleGlobal.h"
0018 #include "AzimuthalProjection_p.h"
0019 
0020 #include <QIcon>
0021 #include <qmath.h>
0022 
0023 #define SAFE_DISTANCE
0024 
0025 namespace Marble
0026 {
0027 
0028 class AzimuthalEquidistantProjectionPrivate : public AzimuthalProjectionPrivate
0029 {
0030   public:
0031     explicit AzimuthalEquidistantProjectionPrivate( AzimuthalEquidistantProjection * parent );
0032 
0033     Q_DECLARE_PUBLIC( AzimuthalEquidistantProjection )
0034 };
0035 
0036 AzimuthalEquidistantProjection::AzimuthalEquidistantProjection()
0037     : AzimuthalProjection( new AzimuthalEquidistantProjectionPrivate( this ) )
0038 {
0039     setMinLat( minValidLat() );
0040     setMaxLat( maxValidLat() );
0041 }
0042 
0043 AzimuthalEquidistantProjection::AzimuthalEquidistantProjection( AzimuthalEquidistantProjectionPrivate *dd )
0044         : AzimuthalProjection( dd )
0045 {
0046     setMinLat( minValidLat() );
0047     setMaxLat( maxValidLat() );
0048 }
0049 
0050 AzimuthalEquidistantProjection::~AzimuthalEquidistantProjection()
0051 {
0052 }
0053 
0054 QString AzimuthalEquidistantProjection::name() const
0055 {
0056     return QObject::tr( "Azimuthal Equidistant" );
0057 }
0058 
0059 QString AzimuthalEquidistantProjection::description() const
0060 {
0061     return QObject::tr( "<p><b>Azimuthal Equidistant Projection</b> (\"fish eye\")</p><p>Applications: Display of seismic and radio data and for use in digital planetariums.</p>" );
0062 }
0063 
0064 QIcon AzimuthalEquidistantProjection::icon() const
0065 {
0066     return QIcon(QStringLiteral(":/icons/map-globe.png"));
0067 }
0068 
0069 AzimuthalEquidistantProjectionPrivate::AzimuthalEquidistantProjectionPrivate( AzimuthalEquidistantProjection * parent )
0070         : AzimuthalProjectionPrivate( parent )
0071 {
0072 }
0073 
0074 qreal AzimuthalEquidistantProjection::clippingRadius() const
0075 {
0076     return 1;
0077 }
0078 
0079 bool AzimuthalEquidistantProjection::screenCoordinates( const GeoDataCoordinates &coordinates,
0080                                              const ViewportParams *viewport,
0081                                              qreal &x, qreal &y, bool &globeHidesPoint ) const
0082 {
0083     const qreal lambda = coordinates.longitude();
0084     const qreal phi = coordinates.latitude();
0085     const qreal lambdaPrime = viewport->centerLongitude();
0086     const qreal phi1 = viewport->centerLatitude();
0087 
0088     qreal cosC = qSin( phi1 ) * qSin( phi ) + qCos( phi1 ) * qCos( phi ) * qCos( lambda - lambdaPrime );
0089     // Prevent division by zero
0090     if (cosC <= 0) {
0091         globeHidesPoint = true;
0092         return false;
0093     }
0094 
0095     qreal c = qAcos(cosC);
0096 
0097     qreal k = cosC == 1 ? 1 : c / qSin( c );
0098 
0099     // Let (x, y) be the position on the screen of the placemark..
0100     x = ( qCos( phi ) * qSin( lambda - lambdaPrime ) ) * k;
0101     y = ( qCos( phi1 ) * qSin( phi ) - qSin( phi1 ) * qCos( phi ) * qCos( lambda - lambdaPrime ) ) * k;
0102 
0103     x *= 2 * viewport->radius() / M_PI;
0104     y *= 2 * viewport->radius() / M_PI;
0105 
0106     const qint64  radius  = clippingRadius() * viewport->radius();
0107 
0108     if (x*x + y*y > radius * radius) {
0109         globeHidesPoint = true;
0110         return false;
0111     }
0112 
0113     globeHidesPoint = false;
0114 
0115     x += viewport->width() / 2;
0116     y = viewport->height() / 2 - y;
0117 
0118     // Skip placemarks that are outside the screen area
0119     return !(x < 0 || x >= viewport->width() || y < 0 || y >= viewport->height());
0120 }
0121 
0122 bool AzimuthalEquidistantProjection::screenCoordinates( const GeoDataCoordinates &coordinates,
0123                                              const ViewportParams *viewport,
0124                                              qreal *x, qreal &y,
0125                                              int &pointRepeatNum,
0126                                              const QSizeF& size,
0127                                              bool &globeHidesPoint ) const
0128 {
0129     pointRepeatNum = 0;
0130     globeHidesPoint = false;
0131 
0132     bool visible = screenCoordinates( coordinates, viewport, *x, y, globeHidesPoint );
0133 
0134     // Skip placemarks that are outside the screen area
0135     if ( *x + size.width() / 2.0 < 0.0 || *x >= viewport->width() + size.width() / 2.0
0136          || y + size.height() / 2.0 < 0.0 || y >= viewport->height() + size.height() / 2.0 )
0137     {
0138         return false;
0139     }
0140 
0141     // This projection doesn't have any repetitions,
0142     // so the number of screen points referring to the geopoint is one.
0143     pointRepeatNum = 1;
0144     return visible;
0145 }
0146 
0147 
0148 bool AzimuthalEquidistantProjection::geoCoordinates( const int x, const int y,
0149                                           const ViewportParams *viewport,
0150                                           qreal& lon, qreal& lat,
0151                                           GeoDataCoordinates::Unit unit ) const
0152 {
0153     const qint64  radius  = viewport->radius();
0154     // Calculate how many degrees are being represented per pixel.
0155     const qreal rad2Pixel = ( 2 * radius ) / M_PI;
0156     const qreal centerLon = viewport->centerLongitude();
0157     const qreal centerLat = viewport->centerLatitude();
0158     const qreal rx = ( - viewport->width()  / 2 + x ) / rad2Pixel;
0159     const qreal ry = (   viewport->height() / 2 - y ) / rad2Pixel;
0160     const qreal c = qMax( qSqrt( rx*rx + ry*ry ), qreal(0.0001) ); // ensure we don't divide by zero
0161     const qreal sinc = qSin(c);
0162 
0163     lon = centerLon + qAtan2( rx*sinc , ( c*qCos( centerLat )*qCos( c ) - ry*qSin( centerLat )*sinc  ) );
0164 
0165     while ( lon < -M_PI ) lon += 2 * M_PI;
0166     while ( lon >  M_PI ) lon -= 2 * M_PI;
0167 
0168     lat = qAsin( qCos(c)*qSin(centerLat) + (ry*sinc*qCos(centerLat))/c );
0169 
0170     if ( unit == GeoDataCoordinates::Degree ) {
0171         lon *= RAD2DEG;
0172         lat *= RAD2DEG;
0173     }
0174 
0175     return true;
0176 }
0177 
0178 }