performance - 実際にその角度を計算せずにベクトルを角度でソートする最速の方法

Question

多くのアルゴリズム ( Graham scanなど) では、ポイントまたはベクトルを角度でソートする必要があります (おそらく、他のポイントから見た角度、つまり差分ベクトルを使用)。この順序は本質的に循環的であり、線形値を計算するためにこの循環がどこで中断されるかは、多くの場合、それほど重要ではありません。しかし、周期的な順序が維持されている限り、実際の角度の値もあまり重要ではありません。そのため、すべてのポイントに対して呼び出しを行うのはatan2無駄かもしれません。角度が厳密に単調な値を計算するためのより高速な方法はありatan2ますか? このような関数は、明らかに「疑似角度」と呼ばれていることがあります。

Answer 1

私はこれをいじり始め、仕様が不完全であることに気付きました。 atan2dx と dy が変化すると、atan2-pi と +pi の間でジャンプするポイントがあるため、不連続性があります。下のグラフは @MvG が提案した 2 つの式を示しており、実際にはどちらも とは異なる場所で不連続性を持っていatan2ます。(注: 線がグラフ上で重ならないように、最初の式に 3 を追加し、別の式に 4 を追加しました)。そのグラフに追加atan2すると、直線 y=x になります。ですから、不連続をどこに置きたいかによって、さまざまな答えがあり得るように私には思えます。本当に を複製したい場合atan2、(このジャンルでの) 答えは

# Input:  dx, dy: coordinates of a (difference) vector.
# Output: a number from the range [-2 .. 2] which is monotonic
#         in the angle this vector makes against the x axis.
#         and with the same discontinuity as atan2
def pseudoangle(dx, dy):
    p = dx/(abs(dx)+abs(dy)) # -1 .. 1 increasing with x
    if dy < 0: return p - 1  # -2 .. 0 increasing with x
    else:      return 1 - p  #  0 .. 2 decreasing with x

これは、使用している言語に符号関数がある場合、sign(dy)(1-p) を返すことで分岐を回避できることを意味します。これにより、-2 を返すことと+2。そして、@MvG の元の方法論でも同じトリックが機能し、sign(dx)(p-1) を返すことができます。

更新以下のコメントで、@MvG はこれを 1 行で C に実装することを提案しています。

pseudoangle = copysign(1. - dx/(fabs(dx)+fabs(dy)),dy)

@MvGはそれがうまく機能すると言い、私には良さそうです:-)。

Answer 2

そのような機能の可能性を 1 つ知っているので、ここで説明します。

# Input:  dx, dy: coordinates of a (difference) vector.
# Output: a number from the range [-1 .. 3] (or [0 .. 4] with the comment enabled)
#         which is monotonic in the angle this vector makes against the x axis.
def pseudoangle(dx, dy):
    ax = abs(dx)
    ay = abs(dy)
    p = dy/(ax+ay)
    if dx < 0: p = 2 - p
    # elif dy < 0: p = 4 + p
    return p

では、なぜこれが機能するのでしょうか。注意すべきことの 1 つは、すべての入力の長さをスケーリングしても出力には影響しないことです。したがって、ベクトルの長さ(dx, dy)は関係なく、方向のみが重要です。最初の象限に集中すると、当面はdx == 1. 次に、 は 0 から無限(つまり) の 1dy/(1+dy)まで単調に増加します。ここで、他の象限も同様に処理する必要があります。が負の場合、初期の も負です。したがって、正の場合、角度で単調な範囲が既にあります。符号を変更して 2 つ追加します。これにより、 の範囲が得られ、全体として の範囲が得られます。何らかの理由で負の数が望ましくない場合、dy == 0dydx == 0dypdx-1 <= p <= 1dx < 01 <= p <= 3dx < 0-1 <= p <= 3elifコメント行を含めることができます。これにより、第 4 象限が から-1…0にシフトし3…4ます。

上記の関数に確立された名前があるかどうか、また誰が最初に公開したかはわかりません。私はかなり前にそれを入手し、あるプロジェクトから次のプロジェクトにコピーしました。しかし、私はこれがウェブ上で発生していることを発見したので、この切り取られた公開は再利用するのに十分であると考えています.

さらなる大文字と小文字の区別を導入せずに範囲 [0 … 4] (実角度 [0 … 2π]) を取得する方法があります。

# Input:  dx, dy: coordinates of a (difference) vector.
# Output: a number from the range [0 .. 4] which is monotonic
#         in the angle this vector makes against the x axis.
def pseudoangle(dx, dy):
    p = dx/(abs(dx)+abs(dy)) # -1 .. 1 increasing with x
    if dy < 0: return 3 + p  #  2 .. 4 increasing with x
    else:      return 1 - p  #  0 .. 2 decreasing with x

Answer 3

角度がそれ自体では必要なく、並べ替えのためだけに必要な場合は、@jjrv アプローチが最適です。ここにジュリアの比較があります

using StableRNGs
using BenchmarkTools

# Definitions
struct V{T}
  x::T
  y::T
end

function pseudoangle(v)
    copysign(1. - v.x/(abs(v.x)+abs(v.y)), v.y)
end

function isangleless(v1, v2)
    a1 = abs(v1.x) + abs(v1.y)
    a2 = abs(v2.x) + abs(v2.y)

    a2*copysign(a1 - v1.x, v1.y) < a1*copysign(a2 - v2.x, v2.y)
end

# Data
rng = StableRNG(2021)
vectors = map(x -> V(x...), zip(rand(rng, 1000), rand(rng, 1000)))

# Comparison
res1 = sort(vectors, by = x -> pseudoangle(x));
res2 = sort(vectors, lt = (x, y) -> isangleless(x, y));

@assert res1 == res2

@btime sort($vectors, by = x -> pseudoangle(x));
  # 110.437 μs (3 allocations: 23.70 KiB)

@btime sort($vectors, lt = (x, y) -> isangleless(x, y));
  # 65.703 μs (3 allocations: 23.70 KiB)

したがって、除算を回避することで、結果の品質を失うことなく、時間がほぼ半分になります。もちろん、より正確な計算のために、を時々isangleless装備する必要がありますが、 についても同じことが言えます。bigfloatpseudoangle

Answer 4

クロス積関数を使用するだけです。一方のセグメントを他方に対して相対的に回転させる方向によって、正または負の数値が得られます。三角関数も除算もありません。速くて簡単。グーグルで検索してください。