Revision control
++ // Copyright (C) 2020 Paolo Angelelli <paolo.angelelli@gmail.com>
++ // Copyright (C) 2022 The Qt Company Ltd.
++ // SPDX-License-Identifier: LicenseRef-Qt-Commercial OR LGPL-3.0-only OR GPL-2.0-only OR GPL-3.0-only
++ // Qt-Security score:significant reason:default
++
++ #include "qdeclarativegeomapitemutils_p.h"
++
++ #include <QPointF>
++ #include <QMatrix4x4>
++ #include <QPainterPath>
++ #include <QPainterPathStroker>
++ #include <QtPositioning/QGeoCoordinate>
++
++ #include <QtPositioning/private/qclipperutils_p.h>
++
++ QT_BEGIN_NAMESPACE
++
++ namespace QDeclarativeGeoMapItemUtils {
++
++ double distanceSqrPointLine(double p0_x
++ , double p0_y
++ , double p1_x
++ , double p1_y
++ , double p2_x
++ , double p2_y)
++ {
++ const double t_x = p2_x - p1_x;
++ const double t_y = p2_y - p1_y;
++ const double p_x = p0_x - p1_x;
++ const double p_y = p0_y - p1_y;
++ const double tsqr = t_x * t_x + t_y * t_y;
++
++ if (tsqr == 0)
++ return qInf();
++
++ double alpha = (p_x * t_x + p_y * t_y) / tsqr;
++ alpha = qBound<double>(0, alpha, 1);
++
++ const double dx = p_x - t_x * alpha;
++ const double dy = p_y - t_y * alpha;
++
++ return dx * dx + dy * dy;
++ }
++
++ void wrapPath(const QList<QGeoCoordinate> &perimeter,
++ const QGeoCoordinate &geoLeftBound,
++ const QGeoProjectionWebMercator &p,
++ QList<QDoubleVector2D> &wrappedPath,
++ QList<QDoubleVector2D> &wrappedPathMinus1,
++ QList<QDoubleVector2D> &wrappedPathPlus1,
++ QDoubleVector2D *leftBoundWrapped)
++ {
++ QList<QDoubleVector2D> path;
++ for (const QGeoCoordinate &c : perimeter)
++ path << p.geoToMapProjection(c);
++ const QDoubleVector2D leftBound = p.geoToMapProjection(geoLeftBound);
++ wrappedPath.clear();
++ wrappedPathPlus1.clear();
++ wrappedPathMinus1.clear();
++ // compute 3 sets of "wrapped" coordinates: one w regular mercator, one w regular mercator +- 1.0
++ for (QDoubleVector2D coord : path) {
++ // We can get NaN if the map isn't set up correctly, or the projection
++ // is faulty -- probably best thing to do is abort
++ if (!qIsFinite(coord.x()) || !qIsFinite(coord.y()))
++ return;
++
++ const bool isPointLessThanUnwrapBelowX = (coord.x() < leftBound.x());
++ // unwrap x to preserve geometry if moved to border of map
++ if (isPointLessThanUnwrapBelowX)
++ coord.setX(coord.x() + 1.0);
++
++ QDoubleVector2D coordP1(coord.x() + 1.0, coord.y());
++ QDoubleVector2D coordM1(coord.x() - 1.0, coord.y());
++
++ wrappedPath.append(coord);
++ wrappedPathPlus1.append(coordP1);
++ wrappedPathMinus1.append(coordM1);
++ }
++ if (leftBoundWrapped)
++ *leftBoundWrapped = leftBound;
++ }
++
++ void wrapPath(const QList<QGeoCoordinate> &perimeter,
++ const QGeoCoordinate &geoLeftBound,
++ const QGeoProjectionWebMercator &p,
++ QList<QDoubleVector2D> &wrappedPath,
++ QDoubleVector2D *leftBoundWrapped)
++ {
++ QList<QDoubleVector2D> path;
++ for (const QGeoCoordinate &c : perimeter)
++ path << p.geoToMapProjection(c);
++ const QDoubleVector2D leftBound = p.geoToMapProjection(geoLeftBound);
++ wrapPath(path, leftBound,wrappedPath);
++ if (leftBoundWrapped)
++ *leftBoundWrapped = leftBound;
++ }
++
++ void wrapPath(const QList<QDoubleVector2D> &path,
++ const QDoubleVector2D &geoLeftBound,
++ QList<QDoubleVector2D> &wrappedPath)
++ {
++ wrappedPath.clear();
++ // compute 3 sets of "wrapped" coordinates: one w regular mercator, one w regular mercator +- 1.0
++ for (QDoubleVector2D coord : path) {
++ // We can get NaN if the map isn't set up correctly, or the projection
++ // is faulty -- probably best thing to do is abort
++ if (!qIsFinite(coord.x()) || !qIsFinite(coord.y()))
++ return;
++
++ const bool isPointLessThanUnwrapBelowX = (coord.x() < geoLeftBound.x());
++ // unwrap x to preserve geometry if moved to border of map
++ if (isPointLessThanUnwrapBelowX)
++ coord.setX(coord.x() + 1.0);
++
++ wrappedPath.append(coord);
++ }
++ }
++
++ void clipPolygon(const QList<QDoubleVector2D> &wrappedPath,
++ const QGeoProjectionWebMercator &p,
++ QList<QList<QDoubleVector2D>> &clippedPaths,
++ QDoubleVector2D *leftBoundWrapped,
++ bool closed)
++ {
++ // 2) Clip bounding box
++ clippedPaths.clear();
++ const QList<QDoubleVector2D> &visibleRegion = p.projectableGeometry();
++ if (visibleRegion.size()) {
++ QClipperUtils clipper;
++ clipper.addSubjectPath(wrappedPath, closed);
++ clipper.addClipPolygon(visibleRegion);
++ clippedPaths = clipper.execute(QClipperUtils::Intersection, QClipperUtils::pftEvenOdd,
++ QClipperUtils::pftEvenOdd);
++
++ if (leftBoundWrapped) {
++ // 2.1) update srcOrigin_ and leftBoundWrapped with the point with minimum X
++ QDoubleVector2D lb(qInf(), qInf());
++ for (const QList<QDoubleVector2D> &path : clippedPaths) {
++ for (const QDoubleVector2D &p : path) {
++ if (p.x() < lb.x() || (p.x() == lb.x() && p.y() < lb.y()))
++ // y-minimization needed to find the same point on polygon and border
++ lb = p;
++ }
++ }
++
++ if (qIsInf(lb.x())) // e.g., when the polygon is clipped entirely
++ return;
++
++ // 2.2) Prevent the conversion to and from clipper from introducing tiny negative offsets which,
++ // in turn will make the geometry wrap around.
++ lb.setX(qMax(leftBoundWrapped->x(), lb.x()));
++
++ *leftBoundWrapped = lb;
++ // srcOrigin_ = p.mapProjectionToGeo(p.unwrapMapProjection(lb));
++ }
++ } else {
++ clippedPaths.append(wrappedPath);
++ }
++ }
++
++ void projectBbox(const QList<QDoubleVector2D> &clippedBbox,
++ const QGeoProjectionWebMercator &p,
++ QPainterPath &projectedBbox)
++ {
++ projectedBbox.clear();
++ bool first = true;
++ for (const auto &coord : clippedBbox) {
++ QDoubleVector2D point = p.wrappedMapProjectionToItemPosition(coord);
++ if (first) {
++ first = false;
++ projectedBbox.moveTo(point.toPointF());
++ } else {
++ projectedBbox.lineTo(point.toPointF());
++ }
++ }
++ projectedBbox.closeSubpath();
++ }
++
++
++ QRectF boundingRectangleFromList(const QList<QDoubleVector2D> &list)
++ {
++ double xMin = qInf();
++ double xMax = -qInf();
++ double yMin = qInf();
++ double yMax = -qInf();
++ for (const auto &coord : list) {
++ xMin = qMin(xMin, coord.x());
++ xMax = qMax(xMax, coord.x());
++ yMin = qMin(yMin, coord.y());
++ yMax = qMax(yMax, coord.y());
++ }
++ return QRectF(xMin, yMin, xMax - xMin, yMax - yMin);
++ }
++
++ QList<QGeoCoordinate> greaterCirclePath(const QList<QGeoCoordinate> &cornerPoints,
++ greaterCirclePathForm form, int N)
++ {
++ QList<QGeoCoordinate> path;
++ if (cornerPoints.empty())
++ return path;
++ path.reserve(N);
++ //If the path has to be closed we include the first coordinate again at the end of the procedure
++ qsizetype lineCount = cornerPoints.size() - ((form == OpenPath) ? 1 : 0);
++ for (qsizetype i = 0; i < lineCount; i++) {
++ path.append(cornerPoints.at(i));
++ const double p_lat = cornerPoints.at(i).latitude() / 180.0 * M_PI; //latitude point p in rad
++ const double p_lon = cornerPoints.at(i).longitude() / 180.0 * M_PI; //longitude point p in rad
++ const double q_lat = cornerPoints.at((i + 1) % cornerPoints.size()).latitude() / 180.0 * M_PI; //latitude point q in rad
++ const double q_lon = cornerPoints.at((i + 1) % cornerPoints.size()).longitude() / 180.0 * M_PI;//longitude point q in rad
++
++ const double c_p_lat = cos(p_lat);
++ const double c_p_lon = cos(p_lon);
++ const double s_p_lat = sin(p_lat);
++ const double s_p_lon = sin(p_lon);
++ const double c_q_lat = cos(q_lat);
++ const double c_q_lon = cos(q_lon);
++ const double s_q_lat = sin(q_lat);
++ const double s_q_lon = sin(q_lon);
++
++ const QDoubleVector3D p(c_p_lat * c_p_lon,
++ c_p_lat * s_p_lon,
++ s_p_lat);
++ const QDoubleVector3D q(c_q_lat * c_q_lon,
++ c_q_lat * s_q_lon,
++ s_q_lat);
++ const double qp = QDoubleVector3D::dotProduct(q, p); //scalar product of q p
++
++ if (qp <= -1.0 || qp >= 1.0)
++ continue;
++
++ // the path from q to p will be parametrized as a combination of points p and q:
++ // x = s*p + t*q.
++ // We will then directly calculate the latitude and longitude of the interpolated point x
++ // from the latitude and longitude of cornerpoints q and p.
++ // The parameters s and t depend on each other to ensure that x is located on the
++ // surface of the earth: s*s + 2*s*t*qp + t*t = 1
++ // Their minimum value is 0, as negative values indicate a path that goes around earth
++ // long way. Their maximum value follows as
++ const double paramMax = sqrt(1 / (1 - qp * qp));
++
++ double t = 0.0;
++ for (int sign = 1; sign > -2.0; sign -= 2.0) {
++ for (; t <= paramMax && t >= 0.0; t += sign * paramMax / N / 2) {
++ // dependence between s and t requires a plus/minus
++ // therefore we solve this equation two times with sign = -1/+1
++ const double s = - t * qp + sign * sqrt(t * t * (qp * qp - 1.0) + 1.0);
++ if (s < 0.0) //s > paramMax will never happen. If s < 0 we are done.
++ break;
++ const double lat = asin(s * s_p_lat + t * s_q_lat);
++ const double lon = atan2(s * c_p_lat * s_p_lon + t * c_q_lat * s_q_lon,
++ s * c_p_lat * c_p_lon + t * c_q_lat * c_q_lon);
++ path.append(QGeoCoordinate(lat * 180.0 / M_PI, lon * 180.0 / M_PI));
++ }
++ t -= paramMax / N / 2;
++ }
++ }
++ if (form == OpenPath)
++ path.append(cornerPoints.last());
++ path.squeeze();
++ return path;
++ }
++
++ } // namespace QDeclarativeGeoMapItemUtils
++
++ QT_END_NAMESPACE