CODEWITHSUNDEEP
c++mathvector
Punal Manalan
Punal Manalan

Reputation: 57

Closest Two 3D Point between two Line Segment of varied Magnitude in Different Plane

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 is a point Line segment AB1&AB2, BP is a point Line segment CD1&CD2.

Point1 and Point2 Closest To each other (Shortest distance between the two line segment)

Now, how can I Find the said two points Point1 and Point2? What method should I use?

THIS IS only partially solved... 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...

Link 1 Link 2

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...

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));
}

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;
}

Upvotes: 1

Views: 1311

Answers (2)

Punal Manalan
Punal Manalan

Reputation: 57

Below Is a "Compact" version of the code from @StefanKssmr which is Here, This "Compact" version can easily be ported to OpenCL

Many thanks to @StefanKssmr for posting the Correct Answer,

void NearestPointBetweenTwoLineSegment(Vector3D AB1, Vector3D AB2, Vector3D CD1, Vector3D CD2, Vector3D& resultSegmentPoint1, Vector3D& resultSegmentPoint2)
{
Vector3D u = VectorMinus(AB2, AB1);
Vector3D v = VectorMinus(CD2, CD1);
Vector3D w = VectorMinus(AB1, CD1);

double    a = DotProduct(u, u);         // always >= 0
double    b = DotProduct(u, v);
double    c = DotProduct(v, v);         // always >= 0
double    d = DotProduct(u, w);
double    e = DotProduct(v, w);
double    sN, sD = (a * c) - (b * b);  // sc = sN / sD, default sD = D >= 0
double    tN, tD = (a * c) - (b * b);  // tc = tN / tD, default tD = D >= 0

float Temp1;
float Temp2;
float Temp3;// Unfortuantely i have no choice but to use this...

//Part 1
Temp1 = (sD < 1e-6f) ? 1.0f : 0.0f;
sN = (1.0f - Temp1) * (b * e - c * d);
sD = ((1.0f - Temp1) * sD) + Temp1;
tN = (Temp1 * e) +  ((1.0f - Temp1) * ((a * e) - (b * d)));
tD = (Temp1 * c) + ((1.0f - Temp1) * tD);

Temp2 = (sN < 0.0f) ? 1.0f : 0.0f;
Temp2 = Temp2 * (1.0f - Temp1);
sN = ((1.0f - Temp2) * sN);
tN = ((1.0f - Temp2) * tN) + (Temp2 * e);
tD = ((1.0f - Temp2) * tD) + (Temp2 * c);

Temp2 = ((sN > sD) ? 1.0f : 0.0f) * (1.0f - Temp2);
Temp2 = Temp2 * (1.0f - Temp1);
sN = ((1.0f - Temp2) * sN) + (Temp2 * sD);
tN = ((1.0f - Temp2) * tN) + (Temp2 * (e + b));
tD = ((1.0f - Temp2) * tD) + (Temp2 * c);

//Part 2.1
Temp1 = (tN < 0.0f) ? 1.0f : 0.0f;
tN = tN * (1.0f - Temp1);

Temp2 = (((-d) < 0.0) ? 1.0f : 0.0f) * Temp1;
sN = (1.0f - Temp2) * sN;//sN = (Temp2 * 0) + ((1.0f - Temp2) * sN);    
Temp3 = ((((-d) > a) ? 1.0f : 0.0f) * (1.0f - Temp2)) * (Temp1);
sN = (Temp3 * sD) + ((1.0f - Temp3) * (sN));

Temp2 = (1.0f - Temp3) * (1.0f - Temp2) * (Temp1);
sN = (Temp2 * (-d)) + ((1.0f - Temp2) * (sN));
sD = (Temp2 * a) + ((1.0f - Temp2) * (sD));

//Part 2.2
Temp1 = ((tN > tD) ? 1.0f : 0.0f) * (1.0f - Temp1);
tN = ((1.0f - Temp1) * tN) + (Temp1 * tD);

Temp2 = (((-d + b) < 0.0) ? 1.0f : 0.0f) * Temp1;
sN = (1.0f - Temp2) * sN;//sN = (Temp2 * 0) + ((1.0f - Temp2) * sN);    
Temp3 = ((((-d + b) > a) ? 1.0f : 0.0f) * (1.0f - Temp2)) * (Temp1);
sN = (Temp3 * sD) + ((1.0f - Temp3) * (sN));

Temp2 = (1.0f - Temp3) * (1.0f - Temp2) * (Temp1);
sN = (Temp2 * (-d)) + ((1.0f - Temp2) * (sN));
sD = (Temp2 * a) + ((1.0f - Temp2) * (sD));

resultSegmentPoint1 = VectorAdd(AB1, VectorMultiply(u, (fabs(sN) < 1e-6f ? 0.0 : sN / sD)));
resultSegmentPoint2 = VectorAdd(CD1, VectorMultiply(v, (fabs(tN) < 1e-6f ? 0.0 : tN / tD)));
}

Upvotes: -1

StefanKssmr
StefanKssmr

Reputation: 1226

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:

  • finding the shortest connection between the lines l1 and l2.
  • finding the shortest connection from either of the points A, B to segment CD (likewise for C,D to segment AB).

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

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

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

You'll find that

df/dsc = [u*sc - v*tc + (A-C)]*u    and    df/dtc = [u*sc - v*tc + (A-C)]*(-v)

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

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

which directly yields

sc = -(A-P)*u/(u*u)    and   rc = -(C-P)*v/(v*v)

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
}

Usage:

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;
}

This prints

P1 = (       -1.24635         1.1212        3.12599)
P2 = (      0.0125125        1.64365        1.49881)

live demo.

Note that this code still requires more testing.

Upvotes: 4

Related Questions