问题描述
我想获得从点(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
和从t
到pq
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
定义.此向量是通过向量OA
和OB
的叉积获得的.
.
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
.因此,可以通过OP
和pp1
的叉积获得与平面 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
.
通过pp1
和pp2
的叉积获得两个平面的线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'.我们可以测试距离PC
和PC'
并获得它们的最小值.
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 ,它依赖于线OA
和OB
之间的角度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几何图形从点到线的垂直地理距离的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!