Let's say A1&A2 are Two 3D points that make a line Segment.
T1,T2,T3 are Three 3D points that make a Triangle Polygon in 3D space.
Let P1 be a point on the Line segment, Let P2 be a point on the Triangle Polygon
P1 And P2 Are Closest To each other
Now, how can I Calculate P1 and P2 which method shall I use?
Right now I know how to find the closest point on a line segment from a point.
and Closest two Points between Two Line Segment.
I use this Below Function to find the closest line segment between two line segment
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
}
This function Below Gets the Closest 2 3D points(line segment) between a Triangle(defined by 3 points) and a Line Segment(Defined by 2 points). These math functions are based on @StefanKssmr 's Answer And @Spektre 's line closest(line l0,triangle t0). Here is the Entire Code Online Compiler with Example. Here is the Visual Representation for the examples
void ClosestPointBetweenTriangleAndLineSegment(Vector3D LineStart, Vector3D LineEnd, Vector3D Triangle0, Vector3D Triangle1, Vector3D Triangle2, Vector3D& ClosestPointOnLine, Vector3D& ClosestPointOnTriangle)
{
ClosestPointOnLine = LineStart;//For FirstPart Calculation Only
ClosestPointOnTriangle;
Vector3D Return1 = LineEnd;//For FirstPart Calculation Only
Vector3D Return2;
NearestPointBetweenPointAndPlane(ClosestPointOnLine, Triangle0, Triangle1, Triangle2, ClosestPointOnTriangle);
NearestPointBetweenPointAndPlane(Return1, Triangle0, Triangle1, Triangle2, Return2);
if (Magnitude(VectorMinus(Return1, Return2)) < Magnitude(VectorMinus(ClosestPointOnLine, ClosestPointOnTriangle)))
{
ClosestPointOnLine = Return1;
ClosestPointOnTriangle = Return2;
}
if (FastNoCheckIsPointWithinTraingle(ClosestPointOnTriangle, Triangle0, Triangle1, Triangle2))
{
return;
}
else
{
NearestPointBetweenTwoLineSegment(LineStart, LineEnd, Triangle0, Triangle1, ClosestPointOnLine, ClosestPointOnTriangle);// Only For First
NearestPointBetweenTwoLineSegment(LineStart, LineEnd, Triangle0, Triangle2, Return1, Return2);
if (Magnitude(VectorMinus(Return1, Return2)) < Magnitude(VectorMinus(ClosestPointOnLine, ClosestPointOnTriangle)))
{
ClosestPointOnLine = Return1;
ClosestPointOnTriangle = Return2;
}
NearestPointBetweenTwoLineSegment(LineStart, LineEnd, Triangle1, Triangle2, Return1, Return2);
if (Magnitude(VectorMinus(Return1, Return2)) < Magnitude(VectorMinus(ClosestPointOnLine, ClosestPointOnTriangle)))
{
ClosestPointOnLine = Return1;
ClosestPointOnTriangle = Return2;
}
}
}
Below are the Helper Function for the above function All are based on common math
#include <iostream>
#include <cmath>
struct Vector3D
{
float x = 0;
float y = 0;
float z = 0;
};
struct Vector4D
{
float x = 0;
float y = 0;
float z = 0;
float d = 0;
};
Vector3D VectorAdd(Vector3D V1, Vector3D V2)
{
return { V1.x + V2.x, V1.y + V2.y, V1.z + V2.z };
}
Vector3D VectorMinus(Vector3D V1, Vector3D V2)
{
return { V1.x - V2.x, V1.y - V2.y, V1.z - V2.z };
}
Vector3D VectorMultiply(Vector3D V1, float Multiplier)
{
return { (V1.x * Multiplier), (V1.y * Multiplier), (V1.z * Multiplier) };
}
Vector3D VectorMultiplyVector(Vector3D V1, Vector3D V2)
{
return { (V1.x * V2.x), (V1.y * V2.y), (V1.z * V2.z) };
}
Vector3D VectorDivide(Vector3D V1, float Divider)
{
return{ (V1.x / Divider), (V1.y / Divider), (V1.z / Divider) };
}
float DotProduct(Vector3D V1, Vector3D V2)
{
return (float)(V1.x * V2.x) + (V1.y * V2.y) + (V1.z * V2.z);
}
Vector3D CrossProduct(Vector3D V1, Vector3D V2)
{
return { (V1.y * V2.z) - (V1.z * V2.y),
(V1.z * V2.x) - (V1.x * V2.z),
(V1.x * V2.y) - (V1.y * V2.x) };
}
Vector3D VectorPower(Vector3D V1, float Power)
{
return { pow(V1.x, Power), pow(V1.y, Power), pow(V1.z, Power) };
}
float VectorTotal(Vector3D V1)
{
return DotProduct(V1, { 1,1,1 });
}
float Magnitude(Vector3D Vector)
{
float CalculationFloat = (float)(sqrt(pow(Vector.x, 2.0) + pow(Vector.y, 2.0) + pow(Vector.z, 2.0)));
if (std::isnan(CalculationFloat) || std::isinf(CalculationFloat))
{
return 0;
}
return CalculationFloat;
}
float AreaOfTriangle3D(Vector3D VA, Vector3D VB, Vector3D VC)
{
return Magnitude(CrossProduct(VectorMinus(VB, VA), VectorMinus(VC, VA))) / 2;
}
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)));
}
bool FastNoCheckIsPointWithinTraingle(Vector3D Point, Vector3D Triangle0, Vector3D Triangle1, Vector3D Triangle2)
{
return fabsf((AreaOfTriangle3D(Point, Triangle1, Triangle2) + AreaOfTriangle3D(Triangle0, Point, Triangle2) + AreaOfTriangle3D(Triangle0, Triangle1, Point)) - AreaOfTriangle3D(Triangle0, Triangle1, Triangle2)) < 1.0f;
}
Vector4D equation_plane(float x1, float y1,
float z1, float x2,
float y2, float z2,
float x3, float y3, float z3)
{
float a1 = x2 - x1;
float b1 = y2 - y1;
float c1 = z2 - z1;
float a2 = x3 - x1;
float b2 = y3 - y1;
float c2 = z3 - z1;
float a = b1 * c2 - b2 * c1;
float b = a2 * c1 - a1 * c2;
float c = a1 * b2 - b1 * a2;
float d = (-a * x1 - b * y1 - c * z1);
//std::cout << std::fixed;
//std::cout << std::setprecision(2);
std::cout << "\nequation of plane is " << a << " x + " << b << " y + " << c << " z + ";
printf("%f", d);
std::cout << " = 0.";
return { a,b,c,d };
}
void NearestPointBetweenPointAndPlane(Vector3D Point, Vector3D Triangle0, Vector3D Triangle1, Vector3D Triangle2, Vector3D& resultSegmentPoint)
{
//Note: Plane equation is x + y + z + d = 0 format//
Vector4D Plane = equation_plane(Triangle0.x, Triangle0.y, Triangle0.z, Triangle1.x, Triangle1.y, Triangle1.z, Triangle2.x, Triangle2.y, Triangle2.z);
Vector3D PlanePartial;
PlanePartial.x = Plane.x;
PlanePartial.y = Plane.y;
PlanePartial.z = Plane.z;
resultSegmentPoint = VectorAdd(Point, VectorMultiply(PlanePartial, -((VectorTotal(VectorMultiplyVector(Point, PlanePartial)) + Plane.d) / VectorTotal(VectorPower(PlanePartial, 2)))));
}
Optional Debug Function
int Clamp(double Number, double Min, double Max)
{
if (Number > Max)
{
return Max;
}
if (Number < Min)
{
return Min;
}
return Number;
}
void PrintVector3Dfor(Vector3D Point, std::string Name, bool InvertXY)
{
int SpaceBar = 50;
SpaceBar = Clamp(SpaceBar - Name.length(), 0, 50);
std::cout << "\n" << Name << " =";
for (int i = 0; i < SpaceBar; ++i)
{
std::cout << " ";
}
std::cout << "(";
if (InvertXY)
{
printf("%f", Point.y);
printf(", %f", Point.x);
}
else
{
printf("%f", Point.x);
printf(", %f", Point.y);
}
printf(", %f )", Point.z);
}
void Printfloatfor(float val, std::string Name)
{
int SpaceBar = 50;
SpaceBar = Clamp(SpaceBar - Name.length(), 0, 50);
std::cout << "\n" << Name << " =";
for (int i = 0; i < SpaceBar; ++i)
{
std::cout << " ";
}
printf(" %f", val);
}