구의 지점에서 라인 세그먼트까지의 거리를 계산하는 방법은 무엇입니까?
https://stackoverflow.com/questions/1299567
-
18-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian문제
지구상에 라인 세그먼트 (큰 원 부분)가 있습니다. 라인 세그먼트는 끝의 좌표로 정의됩니다. 분명히, 두 점은 두 개의 줄 세그먼트를 정의하므로 짧은 세그먼트에 관심이 있다고 가정합니다.
나는 세 번째 포인트가 주어지고 선과 지점 사이의 (가장 짧은) 거리를 찾고 있습니다.
모든 좌표는 경도 위도 (WGS 84)로 제공됩니다.
거리를 어떻게 계산합니까?
합리적인 프로그래밍 언어의 솔루션이 수행됩니다.
해결책
다음은 아이디어를 기반으로 내 자신의 솔루션입니다. Dr. Math에게 물어보십시오. 나는 당신의 의견을 보게되어 기쁩니다.
먼저 면책 조항. 이 솔루션은 구체에 맞습니다. 지구는 구체가 아니며 좌표 시스템 (WGS 84)은 그것이 구라고 가정하지 않습니다. 이것은 단지 근사치 일 뿐이며 실제로 오류라고 추정 할 수는 없습니다. 또한 매우 작은 거리의 경우 모든 것이 단지 코플라나라고 가정하여 좋은 근사치를 얻는 것도 가능합니다. 다시 한 번 나는 거리가 얼마나 "작은"지 모르겠습니다.
이제 사업에. 나는 줄 A, B 및 세 번째 지점 C를 호출하겠습니다. 기본적으로 알고리즘은 다음과 같습니다.
- 먼저 좌표를 직교 좌표로 변환하십시오 (지구 중심에 원점이있는) - 예 : 여기.
다음 3 개의 벡터 제품을 사용하여 C에 가장 가까운 라인 AB의 지점을 계산하십시오.
g = a x b
f = c x g
t = g x f
T를 정규화하고 지구의 반경을 곱하십시오.
- t를 경도 위도로 변환하십시오.
- t와 c 사이의 거리를 계산하십시오. 예 : 여기.
이 단계는 A와 B로 정의 된 C와 Great Circle 사이의 거리를 찾고 있다면 충분합니다. 나와 같이 C와 짧은 선 세그먼트 사이의 거리에 관심이 있으시면 추가 단계를 수행해야합니다. T는 실제로이 세그먼트에 있습니다. 그렇지 않은 경우, 반드시 가장 가까운 지점은 끝 a 또는 b 중 하나입니다. 가장 쉬운 방법은 어느 쪽을 확인하는 것입니다.
일반적으로 세 벡터 제품의 아이디어는 다음과 같습니다. 첫 번째 (g)는 우리에게 A와 B의 큰 원의 평면을 제공합니다 (따라서 A, B 및 원점을 포함하는 평면). 두 번째 (f)는 우리에게 큰 원이 C를 통과하고 G에 수직입니다. T는 F와 G로 정의 된 큰 원의 교차점이며, R에 의한 정규화 및 곱셈에 의해 올바른 길이로 가져옵니다.
다음은이를위한 부분 자바 코드입니다.
큰 원에서 가장 가까운 지점을 찾습니다. 입력 및 출력은 길이 -2 배열입니다. 중간 배열의 길이는 3입니다.
double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
normalize(t);
multiplyByScalar(t, R_EARTH);
return fromCartsian(t);
}
세그먼트에서 가장 가까운 지점 찾기 :
double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (distance(a,c) < distance(b,c)) ? a : c;
}
이것은 우리가 알고있는 지점 T가 A 및 B와 동일한 큰 원에있는 경우이 위대한 원의 짧은 부분에있는 경우 간단한 테스트 방법입니다. 그러나 더 효율적인 방법이 있습니다.
boolean onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.abs(distance(a,b)-distance(a,t)-distance(b,t)) < PRECISION;
}
다른 팁
노력하다 지점에서 큰 원으로의 거리, Dr. Math에서 물어보세요. 여전히 경도/위도를 지구 반경의 구형 좌표로 바꿔야하지만 이것은 좋은 방향처럼 보입니다.
이것은 IDEONE 바이올린으로 허용 된 답변에 대한 완전한 코드입니다 (찾기 여기):
import java.util.*;
import java.lang.*;
import java.io.*;
/* Name of the class has to be "Main" only if the class is public. */
class Ideone
{
private static final double _eQuatorialEarthRadius = 6378.1370D;
private static final double _d2r = (Math.PI / 180D);
private static double PRECISION = 0.1;
// Haversine Algorithm
// source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates
private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
return (1000D * HaversineInKM(lat1, long1, lat2, long2));
}
private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
* Math.pow(Math.sin(dlong / 2D), 2D);
double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;
return d;
}
// Distance between a point and a line
public static void pointLineDistanceTest() {
//line
//double [] a = {50.174315,19.054743};
//double [] b = {50.176019,19.065042};
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
//double [] c = {50.184373,19.054657};
double [] c = {52.008308, 17.542927};
double[] nearestNode = nearestPointGreatCircle(a, b, c);
System.out.println("nearest node: " + Double.toString(nearestNode[0]) + "," + Double.toString(nearestNode[1]));
double result = HaversineInM(c[0], c[1], nearestNode[0], nearestNode[1]);
System.out.println("result: " + Double.toString(result));
}
// source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
private static double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
}
@SuppressWarnings("unused")
private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
}
private static boolean onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
}
// source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
private static double[] toCartsian(double[] coord) {
double[] result = new double[3];
result[0] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.cos(Math.toRadians(coord[1]));
result[1] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.sin(Math.toRadians(coord[1]));
result[2] = _eQuatorialEarthRadius * Math.sin(Math.toRadians(coord[0]));
return result;
}
private static double[] fromCartsian(double[] coord){
double[] result = new double[2];
result[0] = Math.toDegrees(Math.asin(coord[2] / _eQuatorialEarthRadius));
result[1] = Math.toDegrees(Math.atan2(coord[1], coord[0]));
return result;
}
// Basic functions
private static double[] vectorProduct (double[] a, double[] b){
double[] result = new double[3];
result[0] = a[1] * b[2] - a[2] * b[1];
result[1] = a[2] * b[0] - a[0] * b[2];
result[2] = a[0] * b[1] - a[1] * b[0];
return result;
}
private static double[] normalize(double[] t) {
double length = Math.sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
double[] result = new double[3];
result[0] = t[0]/length;
result[1] = t[1]/length;
result[2] = t[2]/length;
return result;
}
private static double[] multiplyByScalar(double[] normalize, double k) {
double[] result = new double[3];
result[0] = normalize[0]*k;
result[1] = normalize[1]*k;
result[2] = normalize[2]*k;
return result;
}
public static void main(String []args){
System.out.println("Hello World");
Ideone.pointLineDistanceTest();
}
}
주석 데이터에 적합합니다.
//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};
가장 가까운 노드는 다음과 같습니다. 50.17493121381319,19.05846668493702
그러나이 데이터에 문제가 있습니다.
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
가장 가까운 노드는 52.00834987257176,17.542691313436357입니다.
두 점으로 지정된 라인은 닫힌 세그먼트가 아니라고 생각합니다.
누군가가 필요하다면 이것은 c#에 포팅 된 loleksy 답변입니다.
private static double _eQuatorialEarthRadius = 6378.1370D;
private static double _d2r = (Math.PI / 180D);
private static double PRECISION = 0.1;
// Haversine Algorithm
// source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates
private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
return (1000D * HaversineInKM(lat1, long1, lat2, long2));
}
private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r)
* Math.Pow(Math.Sin(dlong / 2D), 2D);
double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;
return d;
}
// Distance between a point and a line
static double pointLineDistanceGEO(double[] a, double[] b, double[] c)
{
double[] nearestNode = nearestPointGreatCircle(a, b, c);
double result = HaversineInKM(c[0], c[1], nearestNode[0], nearestNode[1]);
return result;
}
// source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
private static double[] nearestPointGreatCircle(double[] a, double[] b, double [] c)
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
}
private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
}
private static bool onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.Abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
}
// source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
private static double[] toCartsian(double[] coord) {
double[] result = new double[3];
result[0] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Cos(deg2rad(coord[1]));
result[1] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Sin(deg2rad(coord[1]));
result[2] = _eQuatorialEarthRadius * Math.Sin(deg2rad(coord[0]));
return result;
}
private static double[] fromCartsian(double[] coord){
double[] result = new double[2];
result[0] = rad2deg(Math.Asin(coord[2] / _eQuatorialEarthRadius));
result[1] = rad2deg(Math.Atan2(coord[1], coord[0]));
return result;
}
// Basic functions
private static double[] vectorProduct (double[] a, double[] b){
double[] result = new double[3];
result[0] = a[1] * b[2] - a[2] * b[1];
result[1] = a[2] * b[0] - a[0] * b[2];
result[2] = a[0] * b[1] - a[1] * b[0];
return result;
}
private static double[] normalize(double[] t) {
double length = Math.Sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
double[] result = new double[3];
result[0] = t[0]/length;
result[1] = t[1]/length;
result[2] = t[2]/length;
return result;
}
private static double[] multiplyByScalar(double[] normalize, double k) {
double[] result = new double[3];
result[0] = normalize[0]*k;
result[1] = normalize[1]*k;
result[2] = normalize[2]*k;
return result;
}
최대 수천 미터의 거리에서는 구에서 비행기까지 문제를 단순화 할 것입니다. 그런 다음 쉽게 삼각형 계산을 사용할 수 있으므로 문제는 매우 간단합니다.
우리는 포인트 A와 B를 가지고 있으며 라인 AB까지의 거리 X를 찾습니다. 그 다음에:
Location a;
Location b;
Location x;
double ax = a.distanceTo(x);
double alfa = (Math.abs(a.bearingTo(b) - a.bearingTo(x))) / 180
* Math.PI;
double distance = Math.sin(alfa) * ax;
구에서 두 지점 사이의 가장 짧은 거리는 두 지점을 통과하는 큰 원의 작은면입니다. 나는 당신이 이미 이것을 알고 있다고 확신합니다. 여기에는 비슷한 질문이 있습니다 http://www.physicsforums.com/archive/index.php/t-178252.html 그것은 당신이 그것을 수학적으로 모델링하는 데 도움이 될 수 있습니다.
솔직히 말해서 당신이 이것의 코딩 된 예를 얻을 가능성이 확실하지 않습니다.
나는 기본적으로 지금 똑같은 것을 찾고 있습니다. 제가 큰 원의 세그먼트를 갖는 것에 신경 쓰지 않고 오히려 전체 원의 어느 지점까지의 거리를 원한다는 점을 제외하고.
현재 조사중인 두 가지 링크 :
이 페이지 기본적으로 당신이 찾고있는 것 같습니다.
또한 PostGIS 메일 링리스트의 다음 스레드에서 시도는 (1) 2D 평면 (postgis의 line_locate_point)에서 선거에 사용 된 것과 동일한 공식으로 큰 원의 가장 가까운 지점을 결정한 것 같습니다. (2) 그 구형의 거리와 세 번째 지점 사이의 거리를 계산합니다. 수학적으로 단계 (1)이 정확한지 모르겠지만 놀랄 것입니다.
http://postgis.refractions.net/pipermail/postgis-users/2009-july/023903.html
마지막으로, 나는 다음이 "관련"에 링크 된 것을 보았습니다.
