File indexing completed on 2025-01-05 03:59:32

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