不同平面中不同幅度的两条线段之间最近的两个 3D 点(已解决) [英] Closest Two 3D Point between two Line Segment of varied Magnitude in Different Plane(SOLVED)
问题描述
比方说AB1, AB2, CD1, CD2.AB1&AB2 和 CD1&CD2 3D 点构成线段.并且所述线段不在同一平面.
Let's say AB1, AB2, CD1, CD2. AB1&AB2 and CD1&CD2 3D Points makes a Line Segment. And the Said Line segments are Not in the same Plane.
AP是一个点线段AB1&AB2，BP是一个点线段CD1&CD2.
AP is a point Line segment AB1&AB2,
BP is a point Line segment CD1&CD2.
Point1 和 Point2 彼此最近(两条线段之间的最短距离)
现在，我怎样才能找到上述两个点Point1 和Point2?我应该使用什么方法?
Now, how can I Find the said two points Point1 and Point2? What method should I use?
以下只是部分解决对于完整的解决方案，请在此处查看此答案...因为此功能在以下情况下不起作用两条线在同一平面上...感谢@MBo，我遇到了代码和解释的几何金矿！他们有很多源代码贡献者！我从那里挑选了一个它很干净而且很棒！
Below is only partially solved For full solution please See this answer here... because This function does not work when Two Line is on the same plane... Thanks to @MBo I have come across Geometry GoldMine of Code and Explanations! They have Many Source Code Contributors! i picked one from there here it is clean and great!
bool CalculateLineLineIntersection(Vector3D p1, Vector3D p2, Vector3D p3, Vector3D p4, Vector3D& resultSegmentPoint1, Vector3D& resultSegmentPoint2)
{
// Algorithm is ported from the C algorithm of
// Paul Bourke at http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline3d/
resultSegmentPoint1 = { 0,0,0 };
resultSegmentPoint2 = { 0,0,0 };
Vector3D p13 = VectorMinus(p1, p3);
Vector3D p43 = VectorMinus(p4, p3);
/*if (p43.LengthSq() < Math.Epsilon) {
return false;
}*/
Vector3D p21 = VectorMinus(p2, p1);
/*if (p21.LengthSq() < Math.Epsilon) {
return false;
}*/
double d1343 = p13.x * (double)p43.x + (double)p13.y * p43.y + (double)p13.z * p43.z;
double d4321 = p43.x * (double)p21.x + (double)p43.y * p21.y + (double)p43.z * p21.z;
double d1321 = p13.x * (double)p21.x + (double)p13.y * p21.y + (double)p13.z * p21.z;
double d4343 = p43.x * (double)p43.x + (double)p43.y * p43.y + (double)p43.z * p43.z;
double d2121 = p21.x * (double)p21.x + (double)p21.y * p21.y + (double)p21.z * p21.z;
double denom = d2121 * d4343 - d4321 * d4321;
/*if (Math.Abs(denom) < Math.Epsilon) {
return false;
}*/
double numer = d1343 * d4321 - d1321 * d4343;
double mua = numer / denom;
double mub = (d1343 + d4321 * (mua)) / d4343;
resultSegmentPoint1.x = (float)(p1.x + mua * p21.x);
resultSegmentPoint1.y = (float)(p1.y + mua * p21.y);
resultSegmentPoint1.z = (float)(p1.z + mua * p21.z);
resultSegmentPoint2.x = (float)(p3.x + mub * p43.x);
resultSegmentPoint2.y = (float)(p3.y + mub * p43.y);
resultSegmentPoint2.z = (float)(p3.z + mub * p43.z);
return true;
}
到目前为止，我已经尝试了以下所有这些，仅当两条线段具有相同的幅度时才有效......
So Far I have Tried All these Below which works only when both Line segments have the same Magnitude...
我尝试计算两条线段的质心并从中点计算线段上最近的点.(我知道如何计算另一个点的最近点线段)
I tried Calculating the centroid of both line segments and calculating the nearest Point on Segment From the midpoint. (I know how to calculate the Closest Point line segment from another Point)
但这仅适用于两条线段长度相等且线段的中点彼此垂直且质心...
But This only works when Both Line segments are of equal length AND each of Both the Linesegment's MidPoint is perpendicular to Each other and the centroid...
注意:视觉几何 Geogbra3D 用于这些点的视觉表示注意:AB1CD 表示从点 AB1 到线 CD(非线段)
NOTE:Visual Geometry Geogbra3D for a visual representation of these Points
NOTE:AB1CD means From Point AB1 to Line CD(not segment)
AB1 = (6.550000, -7.540000, 0.000000 )
AB2 = (4.540000, -3.870000, 6.000000 )
CD1 = (0.000000, 8.000000, 3.530000 )
CD2 = (0.030000, -7.240000, -1.340000 )
PointCD1AB = (3.117523, -1.272742, 10.246199 )
PointCD2AB = (6.318374, -7.117081, 0.691420 )
PointAB1CD = (0.029794, -7.135321, -1.306549 )
PointAB2CD = (0.019807, -2.062110, 0.314614 )
Magntidue of PointCD1AB - P1LineSegmentCD = 11.866340
Magntidue of PointCD2AB - P2LineSegmentCD = 6.609495
Magntidue of PointAB1CD - P1LineSegmentAB = 6.662127
Magntidue of PointAB2CD - P2LineSegmentAB = 9.186399
Magntidue of PointCD1AB - PointAB1CD = 13.318028
Magntidue of PointCD2AB - PointAB2CD = 8.084965
Magntidue of PointCD1AB - PointAB2CD = 10.433375
Magntidue of PointCD2AB - PointAB1CD = 6.598368
Actual Shortest Point are
Point1 = (0.01, 1.59, 1.48 )
Point2 = (-1.23, 1.11, 3.13 )
Magnitude of Point1 And Point2 = 2.1190799890518526
对于上面的数据，我使用了下面的函数
For the Above Data, I used this Below Function
void NearestPointBetweenTwoLineSegmentOfVariedLength(Vector3D P1LineSegmentAB, Vector3D P2LineSegmentAB, Vector3D P1LineSegmentCD, Vector3D P2LineSegmentCD, Vector3D Testing)
{
/* float Line1Mag = Magnitude(VectorMinus(P1LineSegmentAB, P2LineSegmentAB));
float Line2Mag = Magnitude(VectorMinus(P1LineSegmentCD, P2LineSegmentCD));
P2LineSegmentAB = VectorMinus(P2LineSegmentAB, P1LineSegmentAB);
P1LineSegmentCD = VectorMinus(P1LineSegmentCD, P1LineSegmentAB);
P2LineSegmentCD = VectorMinus(P2LineSegmentCD, P1LineSegmentAB);
P1LineSegmentAB = VectorMinus(P1LineSegmentAB, P1LineSegmentAB);
Vector3D P1P2UnitDirection = GetUnitVector(P2LineSegmentAB, { 0,0,0 });
AngleBetweenTwoVectorsWithCommonUnitVectorAngleOfSecondArgument(P1LineSegmentAB, P2LineSegmentAB, P1P2UnitDirection);*/
Vector3D ReturnVal;
Vector3D PointCD1AB;
Vector3D PointCD2AB;
Vector3D PointAB1CD;
Vector3D PointAB2CD;
NearestPointOnLineFromPoint(P1LineSegmentCD, P1LineSegmentAB, P2LineSegmentAB, PointCD1AB, false);
PrintVector3Dfor(VectorMinus(PointCD1AB, Testing), "PointCD1AB", true);
NearestPointOnLineFromPoint(P2LineSegmentCD, P1LineSegmentAB, P2LineSegmentAB, PointCD2AB, false);
PrintVector3Dfor(VectorMinus(PointCD2AB, Testing), "PointCD2AB", true);
NearestPointOnLineFromPoint(P1LineSegmentAB, P1LineSegmentCD, P2LineSegmentCD, PointAB1CD, false);
PrintVector3Dfor(VectorMinus(PointAB1CD, Testing), "PointAB1CD", true);
NearestPointOnLineFromPoint(P2LineSegmentAB, P1LineSegmentCD, P2LineSegmentCD, PointAB2CD, false);
PrintVector3Dfor(VectorMinus(PointAB2CD, Testing), "PointAB2CD", true);
float m1 = Magnitude(VectorMinus(PointCD1AB, P1LineSegmentCD));
float m2 = Magnitude(VectorMinus(PointCD2AB, P2LineSegmentCD));
float m3 = Magnitude(VectorMinus(PointAB1CD, P1LineSegmentAB));
float m4 = Magnitude(VectorMinus(PointAB1CD, P2LineSegmentAB));
float m5 = Magnitude(VectorMinus(PointCD1AB, PointAB1CD));
float m6 = Magnitude(VectorMinus(PointCD2AB, PointAB2CD));
float m7 = Magnitude(VectorMinus(PointCD1AB, PointAB2CD));
float m8 = Magnitude(VectorMinus(PointCD2AB, PointAB1CD));
Printfloatfor(m1, "Magntidue of PointCD1AB - P1LineSegmentCD");
Printfloatfor(m2, "Magntidue of PointCD2AB - P2LineSegmentCD");
Printfloatfor(m3, "Magntidue of PointAB1CD - P1LineSegmentAB");
Printfloatfor(m4, "Magntidue of PointAB2CD - P2LineSegmentAB");
Printfloatfor(m5, "Magntidue of PointCD1AB - PointAB1CD");
Printfloatfor(m6, "Magntidue of PointCD2AB - PointAB2CD");
Printfloatfor(m7, "Magntidue of PointCD1AB - PointAB2CD");
Printfloatfor(m8, "Magntidue of PointCD2AB - PointAB1CD");
//NearestPointBetweenTwoLineSegmentOfSameLength1(P1LineSegmentAB, P2LineSegmentAB, P1LineSegmentCD, P2LineSegmentCD);
//NearestPointBetweenTwoLineSegmentOfSameLength2(P1LineSegmentAB, P2LineSegmentAB, P1LineSegmentCD, P2LineSegmentCD);
//NearestPointBetweenTwoLineSegmentOfSameLength3(P1LineSegmentAB, P2LineSegmentAB, P1LineSegmentCD, P2LineSegmentCD);
}
void NearestPointOnLineFromPoint(Vector3D Point, Vector3D LineSegmentStart, Vector3D LineSegmentEnd, Vector3D& ReturnVector, bool ClampTheValue)
{
//Get Heading Direction of Capsule from Origin To End
Vector3D CapsuleHeading = VectorMinus(LineSegmentEnd, LineSegmentStart);
float MagnitudeOfLineSegment = Magnitude(CapsuleHeading);
CapsuleHeading = VectorDivide(CapsuleHeading, MagnitudeOfLineSegment);
// Project From Point to Origin
Vector3D Projection = VectorMinus(Point, LineSegmentStart);
float DotProd = DotProduct(Projection, CapsuleHeading);
if (ClampTheValue)
{
DotProd = Clamp(DotProd, 0.0f, MagnitudeOfLineSegment);
}
ReturnVector = VectorAdd(LineSegmentStart, VectorMultiply(CapsuleHeading, DotProd));
}
我已转换此代码 从 C# 到 C++，它没有按预期工作......我不知道是我的代码转换有问题还是代码本身有问题?
I have Converted This Code from C# to C++ and it is not working as intended... I don't know if there is a problem with my code conversion or a problem within the code itself?
Vector3D ClampPointToLine(Vector3D pointToClamp, Vector3D LineStart, Vector3D LineEnd)
{
Vector3D clampedPoint = {0,0,0};
double minX, minY, minZ, maxX, maxY, maxZ;
if (LineStart.x <= LineEnd.x)
{
minX = LineStart.x;
maxX = LineEnd.x;
}
else
{
minX = LineEnd.x;
maxX = LineStart.x;
}
if (LineStart.y <= LineEnd.y)
{
minY = LineStart.y;
maxY = LineEnd.y;
}
else
{
minY = LineEnd.y;
maxY = LineStart.y;
}
if (LineStart.z <= LineEnd.z)
{
minZ = LineStart.z;
maxZ = LineEnd.z;
}
else
{
minZ = LineEnd.z;
maxZ = LineStart.z;
}
clampedPoint.x = (pointToClamp.x < minX) ? minX : (pointToClamp.x > maxX) ? maxX : pointToClamp.x;
clampedPoint.y = (pointToClamp.y < minY) ? minY : (pointToClamp.y > maxY) ? maxY : pointToClamp.y;
clampedPoint.z = (pointToClamp.z < minZ) ? minZ : (pointToClamp.z > maxZ) ? maxZ : pointToClamp.z;
return clampedPoint;
}
void distBetweenLines(Vector3D p1, Vector3D p2, Vector3D p3, Vector3D p4, Vector3D& ClosestPointOnLineP1P2, Vector3D& ClosestPointOnLineP3P4)
{
Vector3D d1;
Vector3D d2;
d1 = VectorMinus(p2, p1);
d2 = VectorMinus(p4, p3);
double eq1nCoeff = (d1.x * d2.x) + (d1.y * d2.y) + (d1.z * d2.z);
double eq1mCoeff = (-(powf(d1.x, 2)) - (powf(d1.y, 2)) - (powf(d1.z, 2)));
double eq1Const = ((d1.x * p3.x) - (d1.x * p1.x) + (d1.y * p3.y) - (d1.y * p1.y) + (d1.z * p3.z) - (d1.z * p1.z));
double eq2nCoeff = ((powf(d2.x, 2)) + (powf(d2.y, 2)) + (powf(d2.z, 2)));
double eq2mCoeff = -(d1.x * d2.x) - (d1.y * d2.y) - (d1.z * d2.z);
double eq2Const = ((d2.x * p3.x) - (d2.x * p1.x) + (d2.y * p3.y) - (d2.y * p2.y) + (d2.z * p3.z) - (d2.z * p1.z));
double M[2][3] = { { eq1nCoeff, eq1mCoeff, -eq1Const }, { eq2nCoeff, eq2mCoeff, -eq2Const } };
int rowCount = 2;
// pivoting
for (int col = 0; col + 1 < rowCount; col++) if (M[col, col] == 0)
// check for zero coefficients
{
// find non-zero coefficient
int swapRow = col + 1;
for (; swapRow < rowCount; swapRow++) if (M[swapRow, col] != 0) break;
if (M[swapRow, col] != 0) // found a non-zero coefficient?
{
// yes, then swap it with the above
double tmp[2];
for (int i = 0; i < rowCount + 1; i++)
{
tmp[i] = M[swapRow][i];
M[swapRow][i] = M[col][i];
M[col][i] = tmp[i];
}
}
else
{
std::cout << "\n the matrix has no unique solution";
return; // no, then the matrix has no unique solution
}
}
// elimination
for (int sourceRow = 0; sourceRow + 1 < rowCount; sourceRow++)
{
for (int destRow = sourceRow + 1; destRow < rowCount; destRow++)
{
double df = M[sourceRow][sourceRow];
double sf = M[destRow][sourceRow];
for (int i = 0; i < rowCount + 1; i++)
M[destRow][i] = M[destRow][i] * df - M[sourceRow][i] * sf;
}
}
// back-insertion
for (int row = rowCount - 1; row >= 0; row--)
{
double f = M[row][row];
if (f == 0) return;
for (int i = 0; i < rowCount + 1; i++) M[row][i] /= f;
for (int destRow = 0; destRow < row; destRow++)
{
M[destRow][rowCount] -= M[destRow][row] * M[row][rowCount]; M[destRow][row] = 0;
}
}
double n = M[0][2];
double m = M[1][2];
Vector3D i1 = { p1.x + (m * d1.x), p1.y + (m * d1.y), p1.z + (m * d1.z) };
Vector3D i2 = { p3.x + (n * d2.x), p3.y + (n * d2.y), p3.z + (n * d2.z) };
Vector3D i1Clamped = ClampPointToLine(i1, p1, p2);
Vector3D i2Clamped = ClampPointToLine(i2, p3, p4);
ClosestPointOnLineP1P2 = i1Clamped;
ClosestPointOnLineP3P4 = i2Clamped;
return;
}
推荐答案
你的问题是求两条线段AB和CD之间的最短连接P1P2.让我们将 l1 定义为通过点 A 和 B 的线，将 l2 定义为通过点的线通过C和D.
Your problem is to find the shortest connection P1P2 between two line segments AB and CD. Let us define l1 as the line which goes through the points A and B and l2 as the line which goes through C and D.
您可以将此问题拆分为多个子任务:
You can split this problem up into several subtasks:
- 找到
l1和l2行之间的最短连接. - 找到从点
A、B到段CD的最短连接(同样对于C,D分割AB).
- finding the shortest connection between the lines
l1andl2. - finding the shortest connection from either of the points
A,Bto segmentCD(likewise forC,Dto segmentAB).
让我们从第一个子任务开始.行 l1，通过 A 和 B，可以由单个标量参数化，例如 sc,
Let's start with the first subtask. THe line l1, going through A and B, can be parametrised by a single scalar, say sc,
l1(sc) = u*sc + A
方向向量u=(B-A).因此，我们也有 l1(0) = A 和 l(1) = B.现在，我们想要找到这条线和另一条通过 C 和 D 的线之间的最小距离，即
with the direction vector u=(B-A).
As a consequence, we also have l1(0) = A and l(1) = B. Now, we want to find the minimal distance between this line and another line going through C and D, i.e.
l2(c) = v*tc + C
与 v = D-C.与另一行类似，我们有 l2(0) = C 和 l(1) = D.现在，我们定义
with v = D-C. In analogy to the other line, we have have l2(0) = C and l(1) = D. Now, we define
f(sc, tc) = 1/2*|l1(sc)-l2(tc)|^2
这只不过是两条线平方之间距离的一半.如果我们现在想最小化这个函数，我们需要满足
which is nothing but half the distance between the two lines squared. If we now want to minimise this function, we need to satisfy
df/dsc = 0 and df/dtc = 0
你会发现
df/dsc = [u*sc - v*tc + (A-C)]*u and df/dtc = [u*sc - v*tc + (A-C)]*(-v)
引入 w=A-C 并排列向量和矩阵产生:
Introducing w=A-C and arranging in vectors and matrices yields:
[ u*u -v*u] * [sc] = -[ w*u]
[-u*v v*v] [tc] [-w*v]
m * result = -rhs
线性系统的解是result = -m^(⁻1)* rhs，其中m^(⁻1)是的逆米.如果 a 和 c 小于 0 或大于 1，则线的最近点在 AB 和 CD 段之外.您也可以返回这些值.
The solution of the linear system is result = -m^(⁻1)* rhs, where m^(⁻1) is the inverse of m. If a and c are less than 0 or greater than 1, the closest point of the lines is outside the segments AB and CD. You might return these values as well.
第二个子任务与这个问题密切相关，但我们最小化
The second subtask is closely related to this problem, but we minimise
f(sc) = 1/2*|l1(sc)-P|^2 and g(tc) = 1/2*|l2(tc)-P|^2
直接产生
sc = -(A-P)*u/(u*u) and rc = -(C-P)*v/(v*v)
如果 sc <0 我们设置 sc = 0 或者如果 sc >1 我们设置 sc = 1(对于 tc 也是如此)，以便在片段上获得分数.
If sc < 0 we set sc = 0 or if sc > 1 we set sc = 1 (and likewise for tc) in order to get points on the segments.
这是我从这里 获取并修改的实现.首先，我们定义一些辅助函数，即向量和一些基本的数学关系.
Here is the implementation, which I took from here and modified it. First, we define some helpers, i.e. vectors and some basic mathematical relations.
template <int dim>
struct Vector
{
std::array<double, dim> components;
};
using Vector2D = Vector<2>;
using Vector3D = Vector<3>;
// subtract
template <int dim>
Vector<dim> operator-(const Vector<dim> &u, const Vector<dim> &v) {
Vector<dim> result(u);
for (int i = 0; i < dim; ++i)
result.components[i] -= v.components[i];
return result;
}
// add
template <int dim>
Vector<dim> operator+(const Vector<dim> &u, const Vector<dim> &v) {
Vector<dim> result(u);
for (int i = 0; i < dim; ++i)
result.components[i] += v.components[i];
return result;
}
// negate
template <int dim>
Vector<dim> operator-(const Vector<dim> &u) {
Vector<dim> result;
for (int i = 0; i < dim; ++i)
result.components[i] = -u.components[i];
return result;
}
// scalar product
template <int dim>
double operator*(const Vector<dim> &u, const Vector<dim> &v) {
double result = 0;
for (int i = 0; i < dim; ++i)
result += u.components[i] * v.components[i];
return result;
}
// scale
template <int dim>
Vector<dim> operator*(const Vector<dim> &u, const double s) {
Vector<dim> result(u);
for (int i = 0; i < dim; ++i)
result.components[i] *= s;
return result;
}
// scale
template <int dim>
Vector<dim> operator*(const double s, const Vector<dim> &u) {
return u*s;
}
// ostream
template <int dim>
std::ostream& operator<< (std::ostream& out, const Vector<dim> &u) {
out << "(";
for (auto c : u.components)
out << std::setw(15) << c ;
out << ")";
return out;
}
这个函数做实际的工作:
This function does the actual work:
std::pair<Vector3D, Vector3D>
shortest_connection_segment_to_segment(const Vector3D A, const Vector3D B,
const Vector3D C, const Vector3D D)
{
Vector3D u = B - A;
Vector3D v = D - C;
Vector3D w = A - C;
double a = u*u; // always >= 0
double b = u*v;
double c = v*v; // always >= 0
double d = u*w;
double e = v*w;
double sc, sN, sD = a*c - b*b; // sc = sN / sD, sD >= 0
double tc, tN, tD = a*c - b*b; // tc = tN / tD, tD >= 0
double tol = 1e-15;
// compute the line parameters of the two closest points
if (sD < tol) { // the lines are almost parallel
sN = 0.0; // force using point A on segment AB
sD = 1.0; // to prevent possible division by 0.0 later
tN = e;
tD = c;
}
else { // get the closest points on the infinite lines
sN = (b*e - c*d);
tN = (a*e - b*d);
if (sN < 0.0) { // sc < 0 => the s=0 edge is visible
sN = 0.0; // compute shortest connection of A to segment CD
tN = e;
tD = c;
}
else if (sN > sD) { // sc > 1 => the s=1 edge is visible
sN = sD; // compute shortest connection of B to segment CD
tN = e + b;
tD = c;
}
}
if (tN < 0.0) { // tc < 0 => the t=0 edge is visible
tN = 0.0;
// recompute sc for this edge
if (-d < 0.0) // compute shortest connection of C to segment AB
sN = 0.0;
else if (-d > a)
sN = sD;
else {
sN = -d;
sD = a;
}
}
else if (tN > tD) { // tc > 1 => the t=1 edge is visible
tN = tD;
// recompute sc for this edge
if ((-d + b) < 0.0) // compute shortest connection of D to segment AB
sN = 0;
else if ((-d + b) > a)
sN = sD;
else {
sN = (-d + b);
sD = a;
}
}
// finally do the division to get sc and tc
sc = (fabs(sN) < tol ? 0.0 : sN / sD);
tc = (fabs(tN) < tol ? 0.0 : tN / tD);
Vector3D P1 = A + (sc * u);
Vector3D P2 = C + (tc * v);
return {P1, P2}; // return the closest distance
}
用法:
int main() {
Vector3D A = {-7.54, 6.55, 0 };
Vector3D B = {4.54, -3.87, 6.0 };
Vector3D C = {0.0, 8.0, 3.53 };
Vector3D D = {0.03, -7.24, -1.34 };
auto [P1, P2] = shortest_connection_segment_to_segment (A, B, C, D);
std::cout << "P1 = " << P1 << std::endl;
std::cout << "P2 = " << P2 << std::endl;
return 0;
}
打印出来
P1 = ( -1.24635 1.1212 3.12599)
P2 = ( 0.0125125 1.64365 1.49881)
现场演示.
请注意，此代码仍需要更多测试.
这篇关于不同平面中不同幅度的两条线段之间最近的两个 3D 点(已解决)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
