2つの2D線分AとBがある場合、AとBを接続する最短の2D線分Cの長さを計算するにはどうすればよいですか?
7 に答える
2 つの線分 A と B がそれぞれ 2 つの点で表されるとします。
A1(x,y)、A2(x,y)で表されるラインA
B1(x,y) B2(x,y) で表される直線 B
最初に、このアルゴリズムを使用して 2 つの線が交差するかどうかを確認します。
交差する場合、2 つの線の間の距離はゼロであり、2 つの線を結ぶ線分が交点になります。
交差しない場合は、次の方法を使用します: http://paulbourke.net/geometry/pointlineplane/間の最短距離を計算します:
- 点 A1 と線 B
- 点 A2 と線 B
- 点 B1 と線 A
- 点 B2 と線 A
これら 4 つの線分の中で最も短いものが答えです。
このページにはあなたが探しているかもしれない情報があります。
上記のAfterlifeとRobParkerのアルゴリズムの一般的な考え方を使用して、2つの任意の2Dセグメント間の最小距離を取得するための一連のメソッドのC++バージョンを次に示します。これは、重なり合うセグメント、平行なセグメント、交差するセグメントと交差しないセグメントを処理します。さらに、さまざまなイプシロン値を使用して、浮動小数点の不正確さから保護します。最後に、最小距離を返すことに加えて、このアルゴリズムは、セグメント2に最も近いセグメント1上のポイント(セグメントが交差する場合は交差ポイントでもあります)を提供します。必要に応じて、[p1、p2]に最も近い[p3、p4]のポイントを返すことも非常に簡単ですが、読者の練習問題として残しておきます:)
// minimum distance (squared) between vertices, i.e. minimum segment length (squared)
#define EPSILON_MIN_VERTEX_DISTANCE_SQUARED 0.00000001
// An arbitrary tiny epsilon. If you use float instead of double, you'll probably want to change this to something like 1E-7f
#define EPSILON_TINY 1.0E-14
// Arbitrary general epsilon. Useful for places where you need more "slop" than EPSILON_TINY (which is most places).
// If you use float instead of double, you'll likely want to change this to something like 1.192092896E-04
#define EPSILON_GENERAL 1.192092896E-07
bool AreValuesEqual(double val1, double val2, double tolerance)
{
if (val1 >= (val2 - tolerance) && val1 <= (val2 + tolerance))
{
return true;
}
return false;
}
double PointToPointDistanceSquared(double p1x, double p1y, double p2x, double p2y)
{
double dx = p2x - p1x;
double dy = p2y - p1y;
return (dx * dx) + (dy * dy);
}
double PointSegmentDistanceSquared( double px, double py,
double p1x, double p1y,
double p2x, double p2y,
double& t,
double& qx, double& qy)
{
double dx = p2x - p1x;
double dy = p2y - p1y;
double dp1x = px - p1x;
double dp1y = py - p1y;
const double segLenSquared = (dx * dx) + (dy * dy);
if (AreValuesEqual(segLenSquared, 0.0, EPSILON_MIN_VERTEX_DISTANCE_SQUARED))
{
// segment is a point.
qx = p1x;
qy = p1y;
t = 0.0;
return ((dp1x * dp1x) + (dp1y * dp1y));
}
else
{
t = ((dp1x * dx) + (dp1y * dy)) / segLenSquared;
if (t <= EPSILON_TINY)
{
// intersects at or to the "left" of first segment vertex (p1x, p1y). If t is approximately 0.0, then
// intersection is at p1. If t is less than that, then there is no intersection (i.e. p is not within
// the 'bounds' of the segment)
if (t >= -EPSILON_TINY)
{
// intersects at 1st segment vertex
t = 0.0;
}
// set our 'intersection' point to p1.
qx = p1x;
qy = p1y;
// Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
// we were doing PointLineDistanceSquared, then qx would be (p1x + (t * dx)) and qy would be (p1y + (t * dy)).
}
else if (t >= (1.0 - EPSILON_TINY))
{
// intersects at or to the "right" of second segment vertex (p2x, p2y). If t is approximately 1.0, then
// intersection is at p2. If t is greater than that, then there is no intersection (i.e. p is not within
// the 'bounds' of the segment)
if (t <= (1.0 + EPSILON_TINY))
{
// intersects at 2nd segment vertex
t = 1.0;
}
qx = p2x;
qy = p2y;
// Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
// we were doing PointLineDistanceSquared, then qx would be (p1x + (t * dx)) and qy would be (p1y + (t * dy)).
}
else
{
// The projection of the point to the point on the segment that is perpendicular succeeded and the point
// is 'within' the bounds of the segment. Set the intersection point as that projected point.
qx = ((1.0 - t) * p1x) + (t * p2x);
qy = ((1.0 - t) * p1y) + (t * p2y);
// for debugging
//ASSERT(AreValuesEqual(qx, p1x + (t * dx), EPSILON_TINY));
//ASSERT(AreValuesEqual(qy, p1y + (t * dy), EPSILON_TINY));
}
// return the squared distance from p to the intersection point.
double dpqx = px - qx;
double dpqy = py - qy;
return ((dpqx * dpqx) + (dpqy * dpqy));
}
}
double SegmentSegmentDistanceSquared( double p1x, double p1y,
double p2x, double p2y,
double p3x, double p3y,
double p4x, double p4y,
double& qx, double& qy)
{
// check to make sure both segments are long enough (i.e. verts are farther apart than minimum allowed vert distance).
// If 1 or both segments are shorter than this min length, treat them as a single point.
double segLen12Squared = PointToPointDistanceSquared(p1x, p1y, p2x, p2y);
double segLen34Squared = PointToPointDistanceSquared(p3x, p3y, p4x, p4y);
double t = 0.0;
double minDist2 = 1E+38;
if (segLen12Squared <= EPSILON_MIN_VERTEX_DISTANCE_SQUARED)
{
qx = p1x;
qy = p1y;
if (segLen34Squared <= EPSILON_MIN_VERTEX_DISTANCE_SQUARED)
{
// point to point
minDist2 = PointToPointDistanceSquared(p1x, p1y, p3x, p3y);
}
else
{
// point - seg
minDist2 = PointSegmentDistanceSquared(p1x, p1y, p3x, p3y, p4x, p4y);
}
return minDist2;
}
else if (segLen34Squared <= EPSILON_MIN_VERTEX_DISTANCE_SQUARED)
{
// seg - point
minDist2 = PointSegmentDistanceSquared(p3x, p3y, p1x, p1y, p2x, p2y, t, qx, qy);
return minDist2;
}
// if you have a point class and/or methods to do cross products, you can use those here.
// This is what we're actually doing:
// Point2D delta43(p4x - p3x, p4y - p3y); // dir of p3 -> p4
// Point2D delta12(p1x - p2x, p1y - p2y); // dir of p2 -> p1
// double d = delta12.Cross2D(delta43);
double d = ((p4y - p3y) * (p1x - p2x)) - ((p1y - p2y) * (p4x - p3x));
bool bParallel = AreValuesEqual(d, 0.0, EPSILON_GENERAL);
if (!bParallel)
{
// segments are not parallel. Check for intersection.
// Point2D delta42(p4x - p2x, p4y - p2y); // dir of p2 -> p4
// t = 1.0 - (delta42.Cross2D(delta43) / d);
t = 1.0 - ((((p4y - p3y) * (p4x - p2x)) - ((p4y - p2y) * (p4x - p3x))) / d);
double seg12TEps = sqrt(EPSILON_MIN_VERTEX_DISTANCE_SQUARED / segLen12Squared);
if (t >= -seg12TEps && t <= (1.0 + seg12TEps))
{
// inside [p1,p2]. Segments may intersect.
// double s = 1.0 - (delta12.Cross2D(delta42) / d);
double s = 1.0 - ((((p4y - p2y) * (p1x - p2x)) - ((p1y - p2y) * (p4x - p2x))) / d);
double seg34TEps = sqrt(EPSILON_MIN_VERTEX_DISTANCE_SQUARED / segLen34Squared);
if (s >= -seg34TEps && s <= (1.0 + seg34TEps))
{
// segments intersect!
minDist2 = 0.0;
qx = ((1.0 - t) * p1x) + (t * p2x);
qy = ((1.0 - t) * p1y) + (t * p2y);
// for debugging
//double qsx = ((1.0 - s) * p3x) + (s * p4x);
//double qsy = ((1.0 - s) * p3y) + (s * p4y);
//ASSERT(AreValuesEqual(qx, qsx, EPSILON_MIN_VERTEX_DISTANCE_SQUARED));
//ASSERT(AreValuesEqual(qy, qsy, EPSILON_MIN_VERTEX_DISTANCE_SQUARED));
return minDist2;
}
}
}
// Segments do not intersect. Find closest point and return dist. No other way at this
// point except to just brute-force check each segment end-point vs opposite segment. The
// minimum distance of those 4 tests is the closest point.
double tmpQx, tmpQy, tmpD2;
minDist2 = PointSegmentDistanceSquared(p3x, p3y, p1x, p1y, p2x, p2y, t, qx, qy);
tmpD2 = PointSegmentDistanceSquared(p4x, p4y, p1x, p1y, p2x, p2y, t, tmpQx, tmpQy);
if (tmpD2 < minDist2)
{
qx = tmpQx;
qy = tmpQy;
minDist2 = tmpD2;
}
tmpD2 = PointSegmentDistanceSquared(p1x, p1y, p3x, p3y, p4x, p4y, t, tmpQx, tmpQy);
if (tmpD2 < minDist2)
{
qx = p1x;
qy = p1y;
minDist2 = tmpD2;
}
tmpD2 = PointSegmentDistanceSquared(p2x, p2y, p3x, p3y, p4x, p4y, t, tmpQx, tmpQy);
if (tmpD2 < minDist2)
{
qx = p2x;
qy = p2y;
minDist2 = tmpD2;
}
return minDist2;
}
Gernot Hoffmannの論文(アルゴリズムとパスカルコード):
クイックヒント:ポイントに基づいて距離を比較する場合は、平方根を実行する必要はありません。
たとえば、P-to-QがQ-to-Rよりも小さいかどうかを確認するには、(擬似コード)を確認します。
square(P.x-Q.x) + square(P.y-Q.y) < square(Q.x-R.x) + square(Q.y-R.y)
Afterlife は、「まず、このアルゴリズムを使用して 2 つの線が交差するかどうかを確認してください」と述べましたが、どのアルゴリズムを意味するのかは明らかにしませんでした。明らかに、重要なのは延長線ではなく、線分の交点です。平行でない線分 (線を定義しない一致端点を除く) は交差しますが、線分間の距離は必ずしもゼロではありません。だから私は彼がそこで「線」ではなく「線分」を意味していたと思います。
Afterlife が提供したリンクは、線 (または線分、または光線) 上の別の任意の点に最も近い点を見つけるための非常に洗練されたアプローチです。これは、各端点から他の線分までの距離を見つけるために機能します (計算されたパラメーター u が、線分または光線の場合は 0 以上、線分では 1 以下になるように制約します) が、そうではありません。実際に交差するため、1 つの線分の内側の点がいずれかの端点よりも近い可能性を処理します。したがって、交差に関する追加のチェックが必要です。
線分との交点を決定するアルゴリズムについては、延長した線分の交点を見つけ (平行であれば完了です)、その点が両方の線分内にあるかどうかを決定する方法があります。交点 T から 2 つの端点までのベクトルのドット積:
((Tx - A1x) * (Tx - A2x)) + ((Ty - A1y) * (Ty - A2y))
これが負 (または「ゼロ」) の場合、T は A1 と A2 の間 (または 1 つの端点) にあります。もう一方の線分についても同様に確認します。いずれかが「ゼロ」より大きい場合、線分は交差しません。もちろん、これは最初に延長された線の交点を見つけることに依存します。これには、各線を方程式として表現し、ガウス縮約などによってシステムを解く必要がある場合があります。
しかし、ベクトル (B1-A1) と (B2-A1) の外積と、ベクトル (B1-A2) の外積と、交点を解く必要のない、より直接的な方法があるかもしれません。 (B2-A2)。これらの外積が同じ方向にある場合、A1 と A2 は直線 B の同じ側にあります。それらが反対方向にある場合、それらはライン B の反対側にあります (0 の場合、一方または両方がライン B上にあります)。同様に、ベクトル (A1-B1) と (A2-B1) および (A1-B2) と (A2-B2) の外積を確認します。これらの外積のいずれかが「ゼロ」である場合、または両方の線分の終点が他の線の反対側にある場合、線分自体が交差する必要があり、そうでない場合は交差しません。
もちろん、座標から 2 つのベクトルの外積を計算するための便利な式が必要です。または、角度 (正または負) を決定できる場合、実際の外積は必要ありません。これは、実際に気にするベクトル間の角度の方向 (または実際には角度のサイン) であるためです。 . しかし、クロス積 (2 次元) の式は単純だと思います。
Cross(V1,V2) = (V1x * V2y) - (V2x * V1y)
これは、3 次元の外積ベクトルの z 軸成分です (初期ベクトルは平面 z=0 にあるため、x 成分と y 成分はゼロでなければなりません)。そのため、符号 (または"ゼロ")。
したがって、これら 2 つの方法のいずれかを使用して、Afterlife が説明するアルゴリズム (リンクを参照) で線分の交差を確認できます。
このページには、2行間の最短距離を見つけるための簡単な説明がありますが、@ stragerのリンクにはいくつかのコードが含まれています(Fortranで!)