本文介绍了使用Boost几何图形从点到线的垂直地理距离的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我想获得从点(t)到线段(p, q)的垂直距离.垂直线可能不与线[p, q]相交.在那种情况下,我想假设延伸线(p, q),然后绘制垂直线以获取距离. p,q,t均为gps坐标.我正在使用升压几何.

I want to get perpendicular distance from a point (t) to a line segment (p, q). The perpendicular may not intersect the line [p, q]. In that case I want to extend the line (p, q) hypothetically and then draw the perpendicular to get the distance. p, q, t are all gps coordinates. I am using boost geometry.

typedef boost::geometry::model::point<
    double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree>
> geo_point;
typedef boost::geometry::model::segment<geo_point> geo_segment;

geo_point p(88.41253929999999, 22.560206299999997);
geo_point q(88.36928063300775, 22.620867969497795);
geo_point t(88.29580956367181, 22.71558662052875);

我已在地图上绘制了这三个位置

I have plotted these three locations on map

我测量两个距离qt和从tpq

I measure two distances qt and distance from t to pq

double dist_qt = boost::geometry::distance(q, t);
std::cout << dist_qt*earth_radius << std::endl;

geo_segment line(p, q);
double perp_dist = boost::geometry::distance(t, line);
std::cout << perp_dist*earth_radius << std::endl;

这两个距离相同.这意味着它不计算垂直距离.相反,它计算了从点到bounds内的直线的shortest距离.

Both of these distances are same. This means it doesn't calculate the perpendicular distance. Rather it calculated shortest distance from a point to a line within bounds.

我如何以这样的方式计算垂直距离,使其无论边界如何都必须垂直?

How can I calculate the perpendicular distance in such a way so that it is necessarily perpendicular irrespective of bounds ?

cpp.sh

推荐答案

此答案可以进行所有计算,无需提升.

This answer do all calculations, without boost.

考虑半径为R = 1的球体.

Consider a sphere of radius R = 1.

A,B点位于大圆圈上.这个大圆gcAB也穿过球体的中心点 O (大圆是必需的).点 A B O 定义了平面PL1.

Points A, B are on a great circle. This great circle gcAB goes also through the center point O of the sphere (required for great circles). Points A, B, O define a plane PL1.

P 也位于一个很大的圆圈中.

Point P also lies in a great circle.

P 到大圆gcAB的最小距离(沿着大圆弧而不是3D直线测量)是圆弧的长度 PC .
大圆gcPC的平面 PL2 垂直于平面 PL1 .

The minimum distance (measured along the arc of a great circle, not along a 3D straight line) from P to the great circle gcAB is the length of the arc PC.
Plane PL2 of great circle gcPC is perpendicular to the plane PL1.

我们想要点 C ,它位于线 OC 中,该线是上述两个平面的交点.

We want point C, which lies in the line OC, which is the intersection of the two mentioned planes.

.
平面 PL1 由其垂直向量pp1定义.此向量是通过向量OAOB叉积获得的.

.
Plane PL1 is defined by its perpendicular vector pp1. This vector is obtained by the cross product of vectors OA and OB.

由于平面 PL2 与平面 PL1 垂直,因此它必须包含矢量pp1.因此,可以通过OPpp1的叉积获得与平面 PL2 垂直的矢量pp2.

Because the plane PL2 is perpendicular to plane PL1, it must contain the vector pp1. So the perpendicular vector pp2 to plane PL2 can be obtained by the cross product of OP and pp1.

通过pp1pp2的叉积获得两个平面的线OC交点中的向量ppi.

A vector ppi in the line OC intersection of both planes is obtained by the cross product of pp1 and pp2.

如果我们归一化矢量ppi,并将其分量乘以地球的半径R,我们将获得点 C 的坐标.
跨产品不是可交换的.这意味着,如果我们交换点A,B,我们将在球体中得到相反的点 C'.我们可以测试距离PCPC'并获得它们的最小值.

If we normalize vector ppi and multiply its components by the radius R of Earth, we get the coordinates of point C.
Cross product is not commutative. This means that if we interchange points A,B we get the opposite point C' in the sphere. We can test distances PC and PC' and get their minimum.

计算大圆距离 Wikipedia链接对于两个点 A B ,它依赖于线OAOB之间的角度a.
为了获得所有角度的最佳精度,我们使用a = atan2(y, x),其中使用半径1,y= sin(a)x= cos(a). sin(a)cos(a)可以分别通过叉积(OA,OB)和点积(OA,OB)计算.

To calculate the great-circle distance Wikipedia link for two points A, B, it relies on the angle a between the lines OA and OB.
For best accuracy on all angles we use a = atan2(y, x) where, using radius 1, y= sin(a) and x= cos(a). sin(a) and cos(a) can be calculated by cross product (OA, OB) and dot product (OA, OB) respectively.

总的来说,我们有以下C ++代码:

Putting all together we have this C++ code:

#include <iostream>
#include <cmath>

const double degToRad = std::acos(-1) / 180;

struct vec3
{
    double x, y, z;
    vec3(double xd, double yd, double zd) : x(xd), y(yd), z(zd) {}
    double length()
    {
        return std::sqrt(x*x + y*y + z*z);
    }
    void normalize()
    {
        double len = length();
        x = x / len;
        y = y / len;
        z = z / len;
    }
};

vec3 cross(const vec3& v1, const vec3& v2)
{
    return vec3( v1.y * v2.z - v2.y * v1.z,
                 v1.z * v2.x - v2.z * v1.x,
                 v1.x * v2.y - v2.x * v1.y );
}

double dot(const vec3& v1, const vec3& v2)
{
    return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}

double GCDistance(const vec3& v1, const vec3& v2, double R)
{
    //normalize, so we can pass any vectors
    vec3 v1n = v1;
    v1n.normalize();
    vec3 v2n = v2;
    v2n.normalize();
    vec3 tmp = cross(v1n, v2n);
    //minimum distance may be in one direction or the other
    double d1 = std::abs(R * std::atan2(tmp.length() , dot(v1n, v2n)));
    double d2 = std::abs(R * std::atan2(tmp.length() , -dot(v1n, v2n)));

    return std::min(std::abs(d1), std::abs(d2));
}

int main()
{
    //Points A, B, and P
    double lon1 = 88.41253929999999  * degToRad;
    double lat1 = 22.560206299999997 * degToRad;
    double lon2 = 88.36928063300775  * degToRad;
    double lat2 = 22.620867969497795 * degToRad;
    double lon3 = 88.29580956367181  * degToRad;
    double lat3 = 22.71558662052875  * degToRad;

    //Let's work with a sphere of R = 1
    vec3 OA(std::cos(lat1) * std::cos(lon1), std::cos(lat1) * std::sin(lon1), std::sin(lat1));
    vec3 OB(std::cos(lat2) * std::cos(lon2), std::cos(lat2) * std::sin(lon2), std::sin(lat2));
    vec3 OP(std::cos(lat3) * std::cos(lon3), std::cos(lat3) * std::sin(lon3), std::sin(lat3));
    //plane OAB, defined by its perpendicular vector pp1
    vec3 pp1 = cross(OA, OB);
    //plane OPC
    vec3 pp2 = cross(pp1, OP);
    //planes intersection, defined by a line whose vector is ppi
    vec3 ppi = cross(pp1, pp2);
    ppi.normalize(); //unitary vector

    //Radious or Earth
    double R = 6371000; //mean value. For more precision, data from a reference ellipsoid is required

    std::cout << "Distance AP = " << GCDistance(OA, OP, R) << std::endl;
    std::cout << "Distance BP = " << GCDistance(OB, OP, R) << std::endl;
    std::cout << "Perpendicular distance (on arc) = " << GCDistance(OP, ppi, R) << std::endl;
}

给出距离对于给定的三个点,AP = 21024.4 BP = 12952.1,PC = 499.493.

Wich gives the distancesAP = 21024.4 BP = 12952.1 and PC= 499.493 for the given three points.

运行代码此处

这篇关于使用Boost几何图形从点到线的垂直地理距离的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!

07-30 04:57