我想得到距一点的垂直距离(t)到一条线段(p, q)。垂线不能与直线相交[p, q]。在这种情况下我想延长线路(p, q)假设，然后绘制垂线以获得距离。 p、q、t 都是 GPS 坐标。我正在使用增强几何。
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 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;
这两个距离是相同的。这意味着它不计算垂直距离。而是计算了shortest点到线内的距离bounds.
如何计算垂直距离，使其无论边界如何都必须垂直？
工作示例在cpp.sh
这个答案完成所有计算，无升压.
考虑半径 R = 1 的球体。
A、B点在a上大圆。这个大圈子gcAB也穿过中心点O球体的直径（大圆所需）。积分A, B, O定义一个平面PL1.
Point P也位于一个大圆圈内。
的最小距离（沿着大圆弧测量，而不是沿着 3D 直线测量）P到大圆gcAB是弧的长度PC.
Plane PL2大圆的gcPC垂直于平面PL1.
我们想要点C，位于该行OC，它是两个提到的平面的交集。
.
Plane PL1由其垂直向量定义pp1。该向量通过以下方式获得叉积向量数OA and OB.
因为飞机PL2垂直于平面PL1，它必须包含向量pp1。所以垂直向量pp2顶级车道PL2可以通过叉积得到OP and pp1.
一个向量ppi在行中OC两个平面的交线由以下叉积获得pp1 and pp2.
If we 正常化 vector ppi并将其分量乘以半径R地球的坐标，我们得到点的坐标C.
叉积不可交换。这意味着如果我们交换点 A、B，我们会得到相反的点C'在球体中。我们可以测试距离PC and PC'并得到他们的最小值。
计算大圆距离 维基百科链接为两点A, B，它依赖于角度a字里行间OA and OB.
为了在所有角度上获得最佳精度，我们使用a = atan2(y, x)其中，使用半径 1，y= sin(a) and x= cos(a). sin(a) and cos(a)可以通过叉积 (OA, OB) 计算点积分别为（OA、OB）。
将所有内容放在一起，我们得到了以下 C++ 代码：
#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。
运行代码here