MapProjection.h

#pragma once

#include <filesystem>
#include "Visual.h"
#include "WebCache.h"
#include "GraphicsModel.h"
#include "ImageStore.h"
#include "Texture.h"
#include "TimerTick.h"
#include "SatelliteModel.h"
#include "Surface.h"
#include "AntiAliasedDrawing.h"

// https://earthobservatory.nasa.gov/features/NightLights/page3.php
// https://visibleearth.nasa.gov/images/57752/blue-marble-land-surface-shallow-water-and-shaded-topography
// https://visibleearth.nasa.gov/images/57752/blue-marble-land-surface-shallow-water-and-shaded-topography/57754l
// https://commons.wikimedia.org/wiki/File:Large_World_Topo_Map_2.png
// https://commons.wikimedia.org/wiki/File:The_earth_at_night.jpg
// https://commons.wikimedia.org/wiki/File:The_earth_at_night_(2).jpg
// https://edwilliams.org/avform147.htm

namespace rose {

#if SDL_VERSION_ATLEAST(2, 0, 10)
    using MapPositionType = float;
#else
    using MapPositionType = int;
#endif

    namespace set {
        static constexpr std::string_view ChronoMapProjection{"MapProjection"};
        static constexpr std::string_view ChronoMapDepiction{"MapDepiction"};
    }

    enum class MapProjectionType {
        Mercator,           
        StationMercator,    
        StationAzimuthal,   
    };

    enum class MapDepiction {
        Terrain,            
        Countries,          
    };

    enum class MapIllumination {
        Day,                
        Night,              
    };

    enum class MapSize {
        Small,              
        Medium,             
        Large,              
        ExtraLarge,         
        Last
    };

    static constexpr Size MapImageSize(MapSize mapSize) {
        switch (mapSize) {
            case MapSize::Small:
                return Size{660,330};
            case MapSize::Medium:
                return Size{1320, 660};
            case MapSize::Large:
                return Size{1980, 990};
            case MapSize::ExtraLarge:
                return Size{2640,1320};
            default:
                break;
        }
        return Size{};
    }

    static constexpr uint32_t MapImageId(MapDepiction mapDepiction, MapSize mapSize, MapIllumination illumination) {
        return (static_cast<uint32_t>(mapSize) << 2u) | (static_cast<uint32_t>(mapDepiction) << 1u) | static_cast<uint32_t>(illumination);
    }

    static std::string MapFileName(MapDepiction mapDepiction, MapSize mapSize, MapIllumination illumination) {
        std::stringstream sstrm{};

        switch (mapDepiction) {
            case MapDepiction::Terrain:
                sstrm << "World_Topo_";
                break;
            case MapDepiction::Countries:
                sstrm << "Countries";
                break;
            default:
                break;
        }

        switch (illumination) {
            case MapIllumination::Day:
                sstrm << "D_";
                break;
            case MapIllumination::Night:
                sstrm << "N_";
                break;
            default:
                break;
        }

        auto size = MapImageSize(mapSize);
        sstrm << size.w << 'x' << size.h;

        sstrm << ".png";
        return sstrm.str();
    }

    template<typename T>
    static constexpr T rad2deg(T r) noexcept {
        return 180. * (r/M_PI);
    }

    template<typename T>
    static constexpr T deg2rad(T d) noexcept {
        return M_PI * (d/180.);
    }

    struct GeoPosition {
        double lat{0.},         
        lon{0.};        
        bool radians{false};    
        bool end{false};        

        constexpr GeoPosition() = default;

        constexpr explicit GeoPosition(bool endPoint) : GeoPosition() {
            end = endPoint;
        }

        constexpr GeoPosition(double latitude, double longitude, bool rad = false) {
            radians = rad;
            if (radians) {
                lat = std::clamp(latitude, -M_PI_2, M_PI_2);
                lon = std::clamp(longitude, -M_PI, M_PI);
            } else {
                lat = std::clamp(latitude, -90., 90.);
                lon = std::clamp(longitude, -180., 180.);
            }
        }

        [[nodiscard]] GeoPosition constexpr toRadians() const {
            if (radians)
                return GeoPosition{*this};
            else {
                GeoPosition g{};
                g.radians = true;
                g.lat = deg2rad(lat);
                g.lon = deg2rad(lon);
                return g;
            }
        }

        [[nodiscard]] GeoPosition constexpr toDegrees() const {
            if (radians) {
                GeoPosition g{};
                g.radians = false;
                g.lat = rad2deg(lat);
                g.lon = rad2deg(lon);
                return g;
            } else
                return GeoPosition{*this};
        }

        [[nodiscard]] double distance(const GeoPosition &other) const {
            auto r = toRadians();
            auto o = other.toRadians();
            return acos(sin(r.lat) * sin(o.lat) + cos(r.lat) * cos(o.lat) * cos(r.lon - o.lon));
        }

        [[nodiscard]] GeoPosition midpoint(const GeoPosition &other, double distance, double fraction = 0.5) const {
            auto r = toRadians();
            auto o = other.toRadians();
            auto A = sin((1. - fraction) * distance) / sin(distance);
            auto B = sin(fraction * distance) / sin(distance);
            auto x = A * cos(r.lat) * cos(r.lon) + B * cos(o.lat) * cos(o.lon);
            auto y = A * cos(r.lat) * sin(r.lon) + B * cos(o.lat) * sin(o.lon);
            auto z = A * sin(r.lat) + B * sin(o.lat);

            return GeoPosition{atan2(z, sqrt(x * x + y * y)), atan2(y, x), true};
        }

        [[nodiscard]] GeoPosition midpoint(const GeoPosition &other, double fraction = 0.5) const {
            return midpoint(other, distance(other), fraction);
        }

        std::ostream &printOn(std::ostream &strm) const {
            auto g = toDegrees();
            return strm << '(' << g.lat << ',' << g.lon << ')';
        }
    };

    static constexpr double EquatorLatitude = 0.;
    static constexpr double TropicLatitude = 23.4365;
    static constexpr double ArcticCircle = 66.5635;
    static constexpr double PrimeMeridian = 0.;
    static constexpr std::array<GeoPosition, 21> InternationalDateLine = {
            {
                    {90., 180.},
                    {75., 180.},
                    {67.7356, -169.25},
                    {65.0189, -169.25},
                    {52.6863, 170.05},
                    {47.8353, 180.},
                    {-0.9, 180.},
                    {-0.9, -159.65},
                    {2.9, -159.65},
                    {2.9, -161.85},
                    {5., -161.85},
                    {5., -155.95},
                    {-7.8, -150.65},
                    {-10., -150.65},
                    {-10., -156.05},
                    {-7.8, -156.05},
                    {-7.8, -178.05},
                    {-15., -172.75},
                    {-45., -172.75},
                    {-51.1815, 180.},
                    {-90., 180.},
            }
    };

    enum class MapOverLayImage : size_t {
        Sun,
        Moon,
        Count,
    };

    struct OverlayImageSpec {
        MapOverLayImage mapOverLayImage;
        std::string_view fileName;
    };

    class MapProjection : public Manager {
    public:
        enum ShortCutCode : uint32_t {
            MercatorProjection,
            StationMercatorProjection,
            AzimuthalProjection,
            TerrainMap,
            CountryMap,
        };

    protected:
        std::shared_ptr<TimerTick> mTimerTick{};

        TickProtocol::slot_type mMapIlluminationTimer{};

//        std::unique_ptr<WebCache> mMapCache{};

        WebCacheProtocol::slot_type mMapSlot{};

        MapProjectionType mProjection{};

        MapDepiction mMapDepiction{MapDepiction::Terrain};

        MapSize mMapSize{MapSize::Small};

        bool mNewSurfaces{};

        Size mMapImgSize{};

        std::array<gm::Surface,2> mMapSurface{};
        std::array<gm::Surface,2> mAzSurface{};

        std::array<gm::Surface,2> mMercatorTemp{};
        std::array<gm::Surface,2> mAzimuthalTemp{};

        std::array<gm::Texture,2> mMercator{};     
        std::array<gm::Texture,2> mAzimuthal{};    

        std::atomic_bool mAbortFuture{};           

        GeoPosition mQth{45.,-75.};

        GeoPosition mQthRad{mQth.toRadians()};

        static constexpr std::array<double, 3> GrayLineCos = {-0.105,
                                                              -0.208,
                                                              -0.309};   
        static constexpr double GrayLinePow = .80;      

        static void
        azimuthalProjection(gm::Surface &projectedSurface, const gm::Surface &mapSurface, Position<int> projected,
                            Position<int> map);

        static std::tuple<bool, double, double>
        xyToAzLatLong(int x, int y, const Size &mapSize, const GeoPosition &location, double sinY, double cosY);

        template<size_t N>
        static bool
        azimuthalProjectionSet(std::atomic_bool &abort, const GeoPosition &qthRad, const Size &mapImageSize, std::array<gm::Surface,N> &projected, std::array<gm::Surface,N> &map) {
            // Compute Azimuthal maps from the Mercator maps
            auto sinY = sin(qthRad.lat);
            auto cosY = cos(qthRad.lat);
            for (int y = 0; y < mapImageSize.h; y += 1) {
                for (int x = 0; x < mapImageSize.w; x += 1) {
                    if (abort) {
                        abort = false;
                        return false;
                    }

                    auto[valid, lat, lon] = xyToAzLatLong(x, y, mapImageSize, qthRad, sinY, cosY);
                    if (valid) {
                        auto xx = std::min(mapImageSize.w - 1,
                                           (int) round((double) mapImageSize.w * ((lon + M_PI) / (2 * M_PI))));
                        auto yy = std::min(mapImageSize.h - 1,
                                           (int) round((double) mapImageSize.h * ((M_PI_2 - lat) / M_PI)));

                        for (auto idx = 0; idx < N; ++idx)
                            azimuthalProjection(projected[idx], map[idx], Position<int>(x, y), Position<int>(xx, yy));
                    }
                }
            }
            return true;
        }

        bool computeAzimuthalMaps() {
            return azimuthalProjectionSet(mAbortFuture, mQthRad, mMapImgSize, mAzSurface, mMapSurface);
        }

        std::future<bool> mComputeAzimuthalMapsFuture;

        bool mMapProjectionsInvalid{true};

        std::future<bool> mForegroundBackgroundFuture;

        bool setForegroundBackground();

        Position<MapPositionType> geoToMap(GeoPosition geo, MapProjectionType projection, int splitPixel, Rectangle &mapRect) const;

    public:
        MapProjection() = delete;

        explicit MapProjection(std::shared_ptr<TimerTick> timerTick, std::filesystem::path& xdgDataPath);

        ~MapProjection() override = default;

        MapProjection(const MapProjection &) = delete;

        MapProjection(MapProjection &&) = delete;

        MapProjection &operator=(const MapProjection &) = delete;

        MapProjection &operator=(MapProjection &&) = delete;

        static constexpr std::string_view id = "MapProjection";
        std::string_view nodeId() const noexcept override {
            return id;
        }

        void draw(gm::Context &context, const Position<int>& containerPosition) override;

        Rectangle layout(gm::Context &context, const Rectangle &screenRect) override;

        void addedToContainer() override;

        void cacheCurrentMaps();

        static std::tuple<double, double> subSolar();

        /*
         * @brief Accessor for Qth
         */
        auto getQth() const {
            return mQth;
        }

        bool mapProjectionsValid() const {
            return !mMapProjectionsInvalid; // && !mForegroundBackgroundFuture.valid();
        }

        auto getProjection() const {
            return mProjection;
        }

        void drawMapItem(const ImageId &mapItem, gm::Context& context, Rectangle mapRectangle, GeoPosition& geoPosition,
                         MapProjectionType projection, int splitPixel);

        int projectionSplitPixel(Size drawSize) const {
            int splitPixel = util::roundToInt((double) drawSize.w * ((mQth.lon) / 360.));
            if (splitPixel < 0)
                splitPixel += drawSize.w;
            return splitPixel;
        }

        void drawMapLine(gm::Context &context, AntiAliasedDrawing &drawing, GeoPosition begin, Rectangle mapRectangle,
                         const std::function<GeoPosition(GeoPosition &, bool fine)> &increment) {
            int splitPixel = 0;
            std::function<bool(const Position<MapPositionType>& p0, const Position<MapPositionType>& p1)> gapTest;

            switch (mProjection) {
                case MapProjectionType::StationAzimuthal:
                    gapTest = [&](const Position<MapPositionType>& p0, const Position<MapPositionType>& p1) -> bool {
                        auto split = mapRectangle.w / 2 + mapRectangle.x;
                        return (p0.x < split && p1.x < split) || (p0.x > split && p1.x > split);
                    };
                    break;
                case MapProjectionType::StationMercator:
                    splitPixel = projectionSplitPixel(mapRectangle.size());
                case MapProjectionType::Mercator:
                    gapTest = [&](const Position<MapPositionType>& p0, const Position<MapPositionType>& p1) -> bool {
                        return abs(p0.x - p1.x) < static_cast<MapPositionType>(mapRectangle.w) / 4 &&
                               abs(p0.y - p1.y) < static_cast<MapPositionType>(mapRectangle.h) / 4;
                    };
                    break;
            }

            auto p0 = geoToMap(begin.toRadians(), mProjection, splitPixel, mapRectangle) + mapRectangle.position();
            auto g0 = begin;
            do {
                auto g1 = increment(g0, false);
                auto p1 = geoToMap(g1.toRadians(), mProjection, splitPixel, mapRectangle) + mapRectangle.position();
                if (gapTest(p0, p1)) {
                    // Draw up to a plotting gap.
                    drawing.renderLine(context, p0, p1);
                }
                else {
                    // Switch to fine increment until the gap is encountered again.
                    for (g1 = increment(g0, true); !g1.end; g1 = increment(g1, true)) {
                        p1 = geoToMap(g1.toRadians(), mProjection, splitPixel, mapRectangle) + mapRectangle.position();
                        if (gapTest(p0, p1)) {
                            drawing.renderLine(context, p0, p1);
                        } else {
                            break;
                        }
                        p0 = p1;
                        g0 = g1;
                    }
                }
                p0 = p1;
                g0 = g1;
            } while (!g0.end);
        }

        void drawInterpolate(gm::Context &context, AntiAliasedDrawing &drawing, Rectangle mapRect, GeoPosition &geo0,
                             GeoPosition &geo1) {
            auto r0 = geo0.toRadians();
            auto r1 = geo1.toRadians();
            int splitPixel = 0;
            std::function<bool(const Position<MapPositionType>& p0, const Position<MapPositionType>& p1)> gapTest;

            auto plotPoints = [this,&splitPixel,&mapRect,&gapTest,&drawing,&context](GeoPosition &g0, GeoPosition &g1) -> bool {
                auto p0 = geoToMap(g0, mProjection, splitPixel, mapRect) + mapRect.position();
                auto p1 = geoToMap(g1, mProjection, splitPixel, mapRect) + mapRect.position();
                if (gapTest(p0, p1)) {
                    drawing.renderLine(context, p0, p1);
                    return false;
                }
                return true;
            };

            switch (mProjection) {
                case MapProjectionType::StationAzimuthal:
                    gapTest = [&](const Position<MapPositionType>& p0, const Position<MapPositionType>& p1) -> bool {
                        auto split = mapRect.w / 2 + mapRect.x;
                        return (p0.x < split && p1.x < split) || (p0.x > split && p1.x > split);
                    };
                    break;
                case MapProjectionType::StationMercator:
                    splitPixel = projectionSplitPixel(mapRect.size());
                case MapProjectionType::Mercator:
                    gapTest = [&](const Position<MapPositionType>& p0, const Position<MapPositionType>& p1) -> bool {
                        return abs(p0.x - p1.x) < mapRect.w / 4 &&
                               abs(p0.y - p1.y) < mapRect.h / 4;
                    };
                    break;
            }

            for (auto r = r0.distance(r1); r > deg2rad(0.25); r = r0.distance(r1)) {
                if (plotPoints(r0, r1)) {
                    auto midPoint = r0.midpoint(r1, r, 0.5);
                    if (plotPoints(r0, midPoint)) {
                        plotPoints(midPoint, r1);
                        r1 = midPoint;
                    } else {
                        r0 = midPoint;
                    }
                } else {
                    break;
                }
            }
        }

        template<typename Iterator>
        void drawMapLine(gm::Context &context, AntiAliasedDrawing &drawing, Rectangle mapRect, Iterator first, Iterator last) {

            double StepSize = deg2rad(3.);
            for (auto idx = first; idx != last - 1; ++idx) {
                auto g1 = (idx + 1)->toRadians();
                auto g0 = idx->toRadians();

                auto dist = g0.distance(g1);
                auto steps = std::max(util::roundToInt(dist / StepSize), 1);
                double f = 0.;
                double fInc = 1. / static_cast<double>(steps);

                auto r0 = g0;
                for (int fIdx = 0; fIdx <= steps; fIdx++) {
                    f = fInc * fIdx;

                    auto r1 = g0.midpoint(g1, dist, f);
                    drawInterpolate(context, drawing, mapRect, r0, r1);
                    r0 = r1;
                }
            }
        }

        void drawLongitude(gm::Context &context, AntiAliasedDrawing &drawing, double longitude, double latitudeBound,
                           Rectangle mapRect) {
            static constexpr double fineInc = 1.;
            static constexpr double coarseInc = 3.;
            double begin = -abs(latitudeBound);
            double end = abs(latitudeBound);
            drawMapLine(context, drawing, GeoPosition{begin, longitude}, mapRect,
                        [&end](GeoPosition &g0, bool fine) -> GeoPosition {
                            auto r = g0;
                            if (fine) {
                                r.lat += fineInc;
                            } else {
                                r.lat += coarseInc;
                            }
                            if (r.lat > end) {
                                r.end = true;
                                r.lat = end;
                            }
                            return r;
                        });
        }

        void drawLatitude(gm::Context &context, AntiAliasedDrawing &drawing, double latitude, Rectangle mapRect) {
            static constexpr double begin = -180.;
            static constexpr double end = 180.;
            static constexpr double fineInc = 1.;
            static constexpr double coarseInc = 3.;
            drawMapLine(context, drawing, GeoPosition{latitude, begin}, mapRect,
                        [](GeoPosition &g0, bool fine) -> GeoPosition {
                            auto r = g0;
                            if (fine) {
                                r.lon += fineInc;
                            } else {
                                r.lon += coarseInc;
                            }
                            r.end = r.lon > end + coarseInc;
                            return r;
                        });
        }
    };
}

inline std::ostream &operator<<(std::ostream &strm, rose::GeoPosition &geoPosition) {
    return geoPosition.printOn(strm);
}