toys/nds/cpp/source/nds.h

#pragma once

#include <cmath>
#include <cstdint>
#include <limits>
#include <tuple>
#include <vector>

constexpr auto DOUBLE_MAX = std::numeric_limits<double>::max();
constexpr auto DOUBLE_MIN = std::numeric_limits<double>::min();

constexpr auto EARTH_EQUATOR_RADIUS = 6378137.0; //   赤道半径 m
constexpr auto EARTH_POLAR_RADIUS = 6356725.0;   //   极半径   m
//
template<typename Integer,
         typename std::enable_if<std::is_integral<Integer>::value, bool>::type = true>
bool equal(const Integer& val1, const Integer& val2, const Integer& eps = 0)
{
    return val1 == val2;
}

template<typename Floating,
         typename std::enable_if<std::is_floating_point<Floating>::value, bool>::type = true>
bool equal(const Floating& val1, const Floating& val2,
           const Floating& eps = std::numeric_limits<Floating>::min())
{
    return std::abs(val1 - val2) < eps;
}

template<typename T>
struct PointT
{
    T x;
    T y;
    T z;

    PointT()
    {
        x = std::numeric_limits<T>::infinity();
        y = std::numeric_limits<T>::infinity();
        z = std::numeric_limits<T>::infinity();
    }

    PointT(const T& _x, const T& _y, const T& _z = std::numeric_limits<T>::infinity())
    {
        x = _x;
        y = _y;
        z = _z;
    }

    template<typename S>
    PointT<S> cast() const
    {
        return PointT<S>(static_cast<S>(x), static_cast<S>(y), static_cast<S>(z));
    }

    bool operator==(const PointT<T>& other) const
    {
        if (is3D())
            return equal(this->x, other.x) && equal(this->y, other.y) && equal(this->z, other.z);
        else
            return equal(this->x, other.x) && equal(this->y, other.y);
    }

    T distance(const PointT<T>& other) const
    {
        return std::hypot(this->x - other.x, this->y - other.y);
    }

    bool is3D() const
    {
        return !(std::isinf(z) || std::isnan(z));
    }

    bool isValid() const
    {
        return !std::isinf(x) && !std::isinf(y);
    }
};

using PointD = PointT<double>;

template<typename T>
struct RectT
{
    T minX;
    T maxX;
    T minY;
    T maxY;
    T minZ;
    T maxZ;

    RectT()
    {
        minX = minY = minZ = std::numeric_limits<T>::infinity();
        maxX = maxY = maxZ = -std::numeric_limits<T>::infinity();
    }

    RectT(const T& minx, const T& maxx, const T& miny, const T& maxy,
          const T& minz = std::numeric_limits<T>::infinity(),
          const T& maxz = -std::numeric_limits<T>::infinity())
    {
        minX = minx;
        maxX = maxx;
        minY = miny;
        maxY = maxy;
        minZ = minz;
        maxZ = maxz;
    }

    template<typename S>
    RectT<S> cast() const
    {
        return RectT<S>(static_cast<S>(minX), static_cast<S>(maxX), static_cast<S>(minY),
                        static_cast<S>(maxY), static_cast<S>(minZ), static_cast<S>(maxZ));
    }

    bool isInvalid() const
    {
        return this->minX >= this->maxX || this->minY >= this->maxY;
    }

    bool is3D() const
    {
        return !std::isinf(minZ) && !std::isinf(maxZ);
    }

    bool operator==(const RectT<T>& other) const
    {
        if (is3D() && other.is3D())
            return equal(this->minX, other.minX) && equal(this->minY, other.minY)
                   && equal(this->minX, other.minX) && equal(this->maxY, other.maxY)
                   && equal(this->minZ, other.minZ) && equal(this->maxZ, other.maxZ);
        else
            return equal(this->minX, other.minX) && equal(this->minY, other.minY)
                   && equal(this->minX, other.minX) && equal(this->maxY, other.maxY);
    }

    bool contain(const PointT<T>& pt) const
    {
        return pt.x <= this->maxX && pt.x >= this->minX && pt.y <= this->maxY && pt.y >= this->minY;
    }

    bool contain(const RectT<T>& other) const
    {
        return this->minX < other.minX && this->minY <= other.minY && this->maxX >= other.maxX
               && this->maxY >= other.maxY;
    }

    bool intersection(const RectT<T>& other) const
    {
        return std::max(this->minX, other.minX) <= std::min(this->maxX, other.maxX)
               && std::max(this->minY, other.minY) <= std::min(this->maxY, other.maxY);
    }

    bool contain3d(const PointT<T>& pt) const
    {
        return pt.x <= this->maxX && pt.x >= this->minX && pt.y <= this->maxY && pt.y >= this->minY
               && pt.z <= this->maxZ && pt.z >= this->minZ;
    }

    bool contain3d(const RectT<T>& other) const
    {
        return this->minX < other.minX && this->minY <= other.minY && this->maxX >= other.maxX
               && this->maxY >= other.maxY && this->minZ < other.minZ && this->maxZ >= other.maxZ;
    }

    void clip(const RectT<T>& other)
    {
        this->minX = std::max(this->minX, other.minX);
        this->minY = std::max(this->minY, other.minY);
        this->minZ = std::max(this->minZ, other.minZ);
        this->maxX = std::min(this->maxX, other.maxX);
        this->maxY = std::min(this->maxY, other.maxY);
        this->maxZ = std::min(this->maxZ, other.maxZ);
    }

    void combine(const RectT<T>& other)
    {
        this->minX = std::min(this->minX, other.minX);
        this->minY = std::min(this->minY, other.minY);
        this->minZ = std::min(this->minZ, other.minZ);
        this->maxX = std::max(this->maxX, other.maxX);
        this->maxY = std::max(this->maxY, other.maxY);
        this->maxZ = std::max(this->maxZ, other.maxZ);
    }

    void combine(const PointT<T>& other)
    {
        this->minX = std::min(this->minX, other.x);
        this->minY = std::min(this->minY, other.y);
        this->minZ = std::min(this->minZ, other.z);
        this->maxX = std::max(this->maxX, other.x);
        this->maxY = std::max(this->maxY, other.y);
        this->maxZ = std::max(this->maxZ, other.z);
    }

    void expand(const T& dis)
    {
        this->minX -= dis;
        this->maxX += dis;
        this->minY -= dis;
        this->maxY += dis;
    }
};

using RectD = RectT<double>;

namespace nds {

// intN转换为int32,主要注意负数情况下前面的符号位
#define INT_N_2_INT32(v, n) (((v) << (32 - (n))) >> (32 - (n)))

// Navinfo mesh id
// see "tile id" in NDS specification or "Mesh" in Navinfo HAD mifg specification
using MeshId = std::uint32_t;

// NDS level id, see NDS specification
using LevelId = std::uint32_t;

// NDS grid id, see NDS specification
using NdsGridId = std::uint32_t;

// In NDS, a coordinate unit corresponds to 360°/(2^32) degrees of longitude or latitude.
// Longitudes range: (-180°, 180°] <==> (–2^31, 2^31]
// Latitudes range: [-90°, 90°] <==> [-2^30, 2^30]
using NdsPoint = std::tuple<std::int32_t, std::int32_t>;

///< Invalid NDS related id
constexpr auto NDS_INVALID_ID = std::numeric_limits<std::uint32_t>::max();

constexpr auto NDS_INVALID_VALUE = std::numeric_limits<double>::infinity();

// The latitude and longitude value corresponding to the NDS coordinate unit
// 360 / math.pow(2, 32)
constexpr auto NDS_PRECISION = 0.00000008381903171539306640625;

constexpr auto MIN_LON = -180.0;
constexpr auto MAX_LON = 180 - NDS_PRECISION;
constexpr auto MIN_LAT = -90.0;
constexpr auto MAX_LAT = 90.0 - NDS_PRECISION;

// NDS level number, highest: 0, lowest: 15, total: 16
constexpr auto NDS_LEVEL_COUNT = 16;

/**
 * @brief get the width of the grid at the specified level in latitude and
 * longitude coordinate, unit: degree
 *
 * @param level
 * @return grid width in degree
 */
double getGridSize(uint32_t level);

/**
 * @brief convert (lon.lat) into NDS point
 *
 * @param lon
 * @param lat
 * @return NdsPoint
 */
NdsPoint lonLat2NdsPoint(double lon, double lat);

/**
 * @brief convert NDS point into (lon,lat)
 *
 * @param ndsPoint
 * @return std::tuple<double, double>: (x,y)
 */
std::tuple<double, double> ndsPoint2lonLat(const NdsPoint& ndsPoint);

/**
 * @brief Get the Mesh Id object
 *
 * @param lon
 * @param lat
 * @param level
 * @return MeshId
 */
MeshId getMeshId(double lon, double lat, uint32_t level);

MeshId getMeshId(const NdsPoint& ndsPt, uint32_t level);

NdsGridId getGridId(MeshId meshId, uint32_t level);

/**
 * @brief Get the Grid Id object
 *
 * @param lon
 * @param lat
 * @param level
 * @return NdsGridId
 */
NdsGridId getGridId(double lon, double lat, uint32_t level);

NdsGridId getGridId(const NdsPoint& ndsPt, uint32_t level);

// 获取和线段p1-->p2相交的meshids
std::vector<MeshId> getSegmentMeshIds(const PointD& p1, const PointD& p2, uint32_t level);

std::tuple<MeshId, std::int32_t> gridId2MeshIdAndLevel(NdsGridId gridId);

/**
 * @brief Get left/top/right/bottom longitude or latitude of specified grid
 *
 * @param gridId
 * @return basic::RectD
 */
RectD getGridRect(NdsGridId gridId);
RectD getMeshRect(MeshId meshId, uint32_t level);
PointD getMeshCenter(MeshId meshId, uint32_t level);

/**
 * @brief 给定一系列NDS grid id,获取所有瓦片的外接矩形
 *
 * @param gridIds
 * @return basic::RectD
 */
RectD getGridsRect(std::vector<NdsGridId> gridIds);

/**
 * @brief 给定一系列mesh ids,获取这些meshs的外接矩形
 *
 * @param meshIds
 * @return basic::RectD
 */
RectD getMeshsRect(std::vector<MeshId> meshIds, uint32_t level);

/**
 * @brief  given a rectangle, get all NDS grid ids of the specified level that
 * intersect with it
 *
 * @param rct 左/下边界为闭区间, 右/上边界为开区间(右上边界值不用来计算meshid)
 * @param level
 * @return std::vector<NdsGridId>: NDS grid ids of specified level
 */
std::vector<NdsGridId> getGridIdInRect(const RectD& rct, uint32_t level);

/**
 * @brief 获取和指定经纬度矩形框相交的mesh ids
 *
 * @param rct 左/下边界为闭区间, 右/上边界为开区间(右上边界值不用来计算meshid)
 * @param level
 * @return std::vector<MeshId>
 */
std::vector<MeshId> getMeshIdInRect(const RectD& rct, uint32_t level);

/**
 * @brief 获取以指定mesh id为中心的 (2n+1)*(2n+1) 的mesh id
 *      注意:返回的顺序从左到右从上到下,即经度从小到大,但是纬度是从大到小
 * @param id
 * @param level
 * @param n
 * @return std::vector<MeshId>
 */
std::vector<MeshId> getAdjacentMeshIds(MeshId id, uint32_t level, uint32_t n);

/**
 * @brief Given a Mesh id and it's level, return 8 mesh ids of the same level
 * around id
 *
 * @param meshId
 * @param level
 * @return std::vector<MeshId>
 */
std::vector<MeshId> get8AroundMeshId(MeshId meshId, uint32_t level);

/**
 * @brief Given a NDS grid id, get 8 grid ids around it
 *      (not valid for level 0 and 1)
 *   Note:: around NDS grid is sequence:
 *      0       1       2
 *      7     gridId    3
 *      6       5       4
 *
 * @param gridId
 * @return std::vector<NdsGridId>: if level ==0 or 1, return all NDS_INVALID_ID
 */
std::vector<NdsGridId> get8AroundGridId(NdsGridId gridId);

} // namespace nds