lng/lat 座標のセットがあります。セット内の任意の 2 点間の最大距離 (場合によっては「最大直径」) を計算する効率的な方法は何でしょうか?
単純な方法は、Haversine 式を使用して各 2 点間の距離を計算し、最大値を取得することですが、これは明らかにうまくスケーリングしません。
編集: ポイントは、モバイル デバイスを持っている人が 1 日のうちに活動したエリアを測定する、十分に小さいエリアにあります。
質問する
5754 次
4 に答える
11
定理 #1: 地球の表面に沿った任意の 2 つの大円距離の順序は、地球をトンネルで貫くポイント間の直線距離の順序と同じです。
したがって、任意の半径の球状地球または指定された形状パラメータの楕円体に基づいて、緯度経度を x、y、z に変換します。これは、ポイントごとに (ポイントのペアごとではなく) サイン/コサインのカップルです。
これで、Haversine 距離の計算に依存しない標準的な 3 次元問題ができました。ポイント間の距離はちょうどユークリッド (3 次元のピタゴラス) です。平方根といくつかの平方根が必要ですが、比較のみを気にする場合は平方根を省略できます。
これを支援するための高度な空間ツリー データ構造が存在する場合があります。またはhttp://www.tcs.fudan.edu.cn/rudolf/Courses/Algorithms/Alg_ss_07w/Webprojects/Qinbo_diameter/2d_alg.htmなどのアルゴリズム(3D メソッドの場合は [次へ] をクリックします)。または C++ コード: http://valis.cs.uiuc.edu/~sariel/papers/00/diameter/diam_prog.html
最大距離のペアが見つかったら、Haversine 式を使用して、そのペアのサーフェスに沿った距離を取得できます。
于 2013-06-01T07:40:59.367 に答える
10
以下は、ポイントの数に応じて二次関数ではなく線形にスケーリングされ、実装が非常に簡単な、有用な近似になると思います。
- 点の重心 M を計算する
- M までの距離が最大となる点 P 0を見つける
- P 0までの距離が最大となる点 P 1を見つける
- P 0と P 1の間の距離で最大直径を近似します
これは、ステップ 3 を N 回繰り返し、P N-1と P Nの間の距離を取ることで一般化できます。
ステップ 1 は、M を経度と緯度の平均として効率的に近似することで実行できます。これは、距離が「小さく」、極が十分に離れている場合に問題ありません。他のステップは正確な距離の式を使用して実行できますが、ポイントの座標が平面上にあると近似できる場合は、はるかに高速です。「遠いペア」(できれば最大距離のペア) が見つかったら、その距離を正確な式で再計算できます。
近似の例は次のようになります。 φ(M) と λ(M) が、Σφ(P)/n と Σλ(P)/n として計算された質量中心の緯度と経度である場合、
- x(P) = (λ(P) - λ(M) + C) cos(φ(P))
- y(P) = φ(P) - φ(M) [これはわかりやすくするためのものであり、単純に y(P) = φ(P) とすることもできます]
ここで、C は通常 0 ですが、点のセットが λ=±180° の線と交差する場合は ± 360° になる可能性があります。最大距離を見つけるには、単に見つける必要があります
- max((x(P N ) - x(P N-1 )) 2 + (y(P N ) - y(P N-1 )) 2 )
(単調なので平方根は必要ありません)
同じ座標変換を使用して、(新しい座標系で) ステップ 1 を繰り返し、開始点をより適切にすることができます。いくつかの条件が満たされている場合、上記の手順 (手順 3 を繰り返さない) は常に「真の遠隔ペア」(私の用語) につながると思います。条件さえわかれば…
編集：
私は他人のソリューションの上に構築するのは嫌いですが、誰かがそうしなければなりません。
上記の 4 つのステップを維持し、ステップ 3 のオプション (ただし、ポイントの典型的な分布によってはおそらく有益) の繰り返しを行い、 Spacedman の解に従って、3D で計算を行うと、極からの近さと距離の制限が克服されます。
- x(P) = sin(φ(P))
- y(P) = cos(φ(P)) sin(λ(P))
- z(P) = cos(φ(P)) cos(λ(P))
(唯一の近似は、これが完全な球に対してのみ当てはまるということです)
重心は x(M) = Σx(P)/n などで与えられ、探す必要のある最大値は
- max((x(P N ) - x(P N-1 )) 2 + (y(P N ) - y(P N-1 )) 2 + (z(P N ) - z(P N-1 ) ) 2 )
つまり、最初に球座標をデカルト座標に変換し、重心から開始して、少なくとも 2 つのステップ (ステップ 2 と 3) で、前の点から最も遠い点を見つけます。距離が長くなる限り、おそらく最大の繰り返し回数でステップ 3 を繰り返すことができますが、これによって極大値から離れることはありません。ポイントが地球全体に広がっている場合、重心から開始してもあまり役に立ちません。
編集2：
アルゴリズムのコアを書き留めるのに十分な R を学習しました (データ分析に適した言語です!)
平面近似の場合、λ=±180° ラインの周りの問題を無視します。
# input: lng, lat (vectors)
rad = pi / 180;
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i = which.max((x - mean(x))^2 + (y )^2)
j = which.max((x - x[i] )^2 + (y - y[i])^2)
# output: i, j (indices)
i私の PC では、インデックスを見つけてj1000000 ポイントを見つけるのに 1 秒もかかりません。
次の 3D バージョンは少し遅くなりますが、ポイントの任意の分布に対して機能します (λ=±180° の線が交差する場合は修正する必要はありません)。
# input: lng, lat
rad = pi / 180
x = sin(lat * rad)
f = cos(lat * rad)
y = sin(lng * rad) * f
z = cos(lng * rad) * f
i = which.max((x - mean(x))^2 + (y - mean(y))^2 + (z - mean(z))^2)
j = which.max((x - x[i] )^2 + (y - y[i] )^2 + (z - z[i] )^2)
k = which.max((x - x[j] )^2 + (y - y[j] )^2 + (z - z[j] )^2) # optional
# output: j, k (or i, j)
の計算は、データと要件に応じて省略kすることができます (つまり、結果は と で与えられiます)。j一方、私の実験では、それ以上の指数を計算しても意味がないことが示されています。
いずれにせよ、結果のポイント間の距離は、セットの「直径」の下限である推定値であることに注意してください。 )
編集3：
残念ながら、平面近似の相対誤差は、極端な場合には 1-1/√3 ≅ 42.3% にもなり、非常にまれではありますが、許容できない場合があります。このアルゴリズムは、私がコンパスと定規から導き出した約 20% の上限を持つように変更できます (分析ソリューションは面倒です)。変更されたアルゴリズムは、局所的に最大の距離を持つ点のペアを見つけ、同じ手順を繰り返しますが、今回は最初のペアの中点から開始し、別のペアを見つける可能性があります。
# input: lng, lat
rad = pi / 180
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i.n_1 = 1 # n_1: n-1
x.n_1 = mean(x)
y.n_1 = 0 # = mean(y)
s.n_1 = 0 # s: square of distance
repeat {
s = (x - x.n_1)^2 + (y - y.n_1)^2
i.n = which.max(s)
x.n = x[i.n]
y.n = y[i.n]
s.n = s[i.n]
if (s.n <= s.n_1) break
i.n_1 = i.n
x.n_1 = x.n
y.n_1 = y.n
s.n_1 = s.n
}
i.m_1 = 1
x.m_1 = (x.n + x.n_1) / 2
y.m_1 = (y.n + y.n_1) / 2
s.m_1 = 0
m_ok = TRUE
repeat {
s = (x - x.m_1)^2 + (y - y.m_1)^2
i.m = which.max(s)
if (i.m == i.n || i.m == i.n_1) { m_ok = FALSE; break }
x.m = x[i.m]
y.m = y[i.m]
s.m = s[i.m]
if (s.m <= s.m_1) break
i.m_1 = i.m
x.m_1 = x.m
y.m_1 = y.m
s.m_1 = s.m
}
if (m_ok && s.m > s.n) {
i = i.m
j = i.m_1
} else {
i = i.n
j = i.n_1
}
# output: i, j
3D アルゴリズムも同様の方法で変更できます。(2D と 3D の両方の場合で) 2 番目のポイント ペアの中点からやり直すことができます (見つかった場合)。この場合の上限は「読者の演習として残されています」:-)。
修正されたアルゴリズムを (あまりにも) 単純なアルゴリズムと比較すると、正規分布と正方形一様分布の場合、処理時間がほぼ 2 倍になり、平均誤差が 0.6% から 0.03% (桁違い) に減少しました。 . 中間点からさらに再開すると、平均誤差はわずかに改善されますが、最大誤差はほぼ同じになります。
編集4：
この記事はまだ勉強中ですが、コンパスと直定規で見つけた 20% は、実際には 1-1/√(5-2√3) ≅ 19.3% のようです。
于 2013-05-31T20:55:33.087 に答える
3
これは、あなたが言うように（あなたが言うように）うまくスケーリングしない単純な例ですが、Rでソリューションを構築するのに役立つかもしれません.
## lonlat points
n <- 100
d <- cbind(runif(n, -180, 180), runif(n, -90, 90))
library(sp)
## distances on WGS84 ellipsoid
x <- spDists(d, longlat = TRUE)
## row, then column index of furthest points
ind <- c(row(x)[which.max(x)], col(x)[which.max(x)])
## maps
library(maptools)
data(wrld_simpl)
plot(as(wrld_simpl, "SpatialLines"), col = "grey")
points(d, pch = 16, cex = 0.5)
## draw the points and a line between on the page
points(d[ind, ], pch = 16)
lines(d[ind, ], lwd = 2)
## for extra credit, draw the great circle on which the furthest points lie
library(geosphere)
lines(greatCircle(d[ind[1], ], d[ind[2], ]), col = "firebrick")
このgeosphereパッケージは、必要に応じて距離計算のオプションをさらに提供します。ここで使用されている詳細については、 を参照?spDistsしspてください。
于 2013-06-01T00:05:51.607 に答える
3
これらのポイントが地球の十分に小さい部分に配置されるかどうかはわかりません。真にグローバルなポイント セットの場合、私の最初の推測では、単純な O(n^2) アルゴリズムを実行し、空間インデックス (R* ツリー、8 進ツリーなど) を使用してパフォーマンスを向上させることができます。アイデアは、距離行列の三角形の n*(n-1) リストを事前に生成し、チャンクで高速距離ライブラリにフィードして、I/O を最小限に抑え、チャーンを処理することです。Haversine で問題ありません。Vincenty の方法を使用することもできます (実行時間の最大の要因は、Vincenty の式の (固定回数の) 反復ではなく、2 次の複雑さです)。ちなみに、実際には、このようなものに R は必要ありません。
編集 #2: Barequet-Har-Peledアルゴリズム(Spacedman の返信で指摘されているように) は、e>0 に対して O((n+1/(e^3))log(1/e)) の複雑さを持ち、探索する価値があります。
準平面問題の場合、これは「凸包の直径」として知られており、次の 3 つの部分があります。
- O(n*log(n)) であるグラハム スキャンで凸包を計算する- 実際、ポイントを横メルカトル図法に変換してみる必要があります (データ セット内のポイントの重心を使用)。
- Rotating Calipersアルゴリズムによる対蹠点の検出- 線形 O(n)。
- すべての対蹠ペアの中で最大の距離を見つける - 線形検索、O(n)。
疑似コードとディスカッションへのリンク: http://fredfsh.com/2013/05/03/convex-hull-and-its-diameter/
ここで関連する質問に関する議論も参照してください: https://gis.stackexchange.com/questions/17358/how-can-i-find-the-farthest-point-from-a-set-of-existing-points
編集: Spacedman のソリューションは、私にMalandain-Boissonnatアルゴリズムを指摘しました (こちらの pdf の論文を参照してください)。ただし、これはブルートフォース ナイーブ O(n^2) アルゴリズムより悪いか同じです。
于 2013-06-01T18:39:40.383 に答える