0

I have a case where I need to find the closest point to a sphere's center on a LOT of triangles. However, the accuracy isn't quite so important... it needs to be vaguely accurate, but an error rate of like 5% would be acceptable.

Currently I'm using two tests-- IsPointInTriangle and closest point on line for all three triangle edges. Specifically, this code:

Vector Math::GetClosestPointOnLine(Vector thePos, Vector theL1, Vector theL2)
{
    //
    // Do NOT try to optimize this like the 2D version... when I did it screwed up all the Hamsterball physics!
    //
    Vector aC=thePos-theL1;
    Vector aV=theL2-theL1; 

    float aD=aV.Length();
    aV/=aD;

    float aT=Dot(aV,aC);

    if (aT<0.0f) return (theL1);
    if (aT>aD) return (theL2);


    // Return the point between ‘a’ and ‘b’
    //set length of V to t. V is normalized so this is easy
    aV*=aT;
    return (theL1+aV);  
}

bool Math::IsPointInTriangle(Vector thePos, Vector theT1, Vector theT2, Vector theT3)

//bool IsPointInTriangle(const VECTOR& point,
//  const VECTOR& pa,const VECTOR& pb, const VECTOR& pc)
{
    Vector e10=theT2-theT1;
    Vector e20=theT3-theT1;
    float a = e10.Dot(e10);
    float b = e10.Dot(e20);
    float c = e20.Dot(e20);
    float ac_bb=(a*c)-(b*b);
    Vector vp(thePos.mX-theT1.mX, thePos.mY-theT1.mY, thePos.mZ-theT1.mZ);
    float d = vp.Dot(e10);
    float e = vp.Dot(e20);
    float x = (d*c)-(e*b);
    float y = (e*a)-(d*b);
    float z = x+y-ac_bb;
    return (( in(z)& ~(in(x)|in(y)) ) & 0x80000000)!=0;
}

Both of these are pretty fast, aside from a sqrt in closestpointonline. But it seems like there is probably a faster way. Any math wizards who can tell me a speedier-- at the cost of accuracy-- way to roughly say "closestpointontriangle, with maybe a 5% error threshold?"

Edit to add: The only alternative way I can think to do this would be to make four planes (plane of triangle + plane of edges in the direction of triangle normal) and fall to the triangle's plane then clip to all edge planes. But I'm not sure if that would be faster or not.

KiraHoneybee
  • 495
  • 3
  • 12

2 Answers2

1

Not an algorithmic optimization, but you can get rid of the sqrt entirely in GetClosestPointOnLine.

This is what you are computing in the first part:

aC = thePos - theL1
aV1 = theL2 - theL1 // first value of aV
aD = sqrt(aV1x*aV1x + aV1y*aV1y + aV1z*aV1z)
aV2 = aV1/aD // second value of aV
aT = aV2x*aCx + aV2y*aCy + aV2z*aCz = (aV1x*aCx + aV1y*aCy) / aD
// comparisons aT < 0.0f and aT > aD
aV3 = aV2*aT = aV1 * (aV1x*aCx + aV1y*aCy) / (aD*aD) // third value of aV

We will assume aD > 0, otherwise you already have UB in aV/=aD.
For aT < 0.0f, you don't need the 1/aD scale factor.
For aT > aD, you can instead write aT*aD > aD*aD, which is:

aV1x*aCx + aV1y*aCy > aV1x*aV1x + aV1y*aV1y

If you drop the sqrt in aD, remove the aV/=aD line, and change the last return value (aV3) accordingly, then the computation becomes:

aC = thePos - theL1
aV1 = theL2 - theL1
aD = aV1x*aV1x + aV1y*aV1y
aT = aV1x*aCx + aV1y*aCy
// comparisons aT < 0.0f and aT > aD
aV3 = aV1*aT/aD

and the code:

Vector Math::GetClosestPointOnLine(Vector thePos, Vector theL1, Vector theL2)
{
    Vector aC=thePos-theL1;
    Vector aV=theL2-theL1; 

    float aD = aV.LengthSquared(); // or Dot(aV, aV)
    float aT = Dot(aV, aC);

    if (aT<0.0f) return (theL1);
    if (aT>aD) return (theL2);

    aV*=aT/aD;
    return (theL1+aV);  
}

I did a quick test with a few random vectors, taking all three paths, to make sure both versions return the same vector.

Nelfeal
  • 12,593
  • 1
  • 20
  • 39
  • This works perfectly in 2D space-- in 3D space it starts producing acceleration of the "closest point" toward the endpoints as you get further from the center of the line. If you convert it to 3D and make one line {-1,-1,-1} - {1,1,1} you'll see that the result of closestpoint2 vs closestpoint begins to accelerate toward the endpoints. This is the kind of thing I'm looking for though, if there's a way to make this work in 3D! Reducing all the equations down to minimum. – KiraHoneybee Jul 30 '23 at 14:19
  • I'm not sure what you mean. [It seems to work in 3D as well.](https://godbolt.org/z/rsGnhYWdj) – Nelfeal Jul 30 '23 at 14:37
  • When I implement it (pasting your code) I get these results: https://imgur.com/a/vmqSflK red = the line, blue is the test point, green = closest point via sqrt algo, yellow=closest point via your algo. The line is {-1,-1,-1} to {1,1,1}. The blue dot starts at {0,0,0} and I am moving it only in the x and y axis, leaving z zero. – KiraHoneybee Jul 30 '23 at 14:45
  • Ah, never mind, I kept the original code, so the lengthsquared and dot weren't working incorporating z. You are correct, this works, and is faster! Before I check yours as the answer, do you see any better way to crunch down closest point on triangle rather than doing IsPointInTriangle() and then ClosestPoint to each line? Is there some nice way to crunch that down? – KiraHoneybee Jul 30 '23 at 14:51
  • Honestly I don't quite understand to problem statement. What I got is that you want to find the point in a triangle that is closest to another point in 3D space, and you want to do that for many different triangles. But I don't understand what `IsPointInTriangle` does exactly or why you use it. It would help to have a minimal reproducible example. – Nelfeal Jul 30 '23 at 15:37
0

Let the vertices of your triangle be A, B and C.

Find the centroid of the triangle P (we need an arbitrary point inside the triangle, and this is an easy one to calculate as the mean of the vertices).

Now say you want to find the point in the triangle closest to point X. First project X into the plane of the triangle (call this point Y).

Now calculate norm(P-A).(P-Y), norm(P-B).(P-Y), norm(P-C).(P-Y). Find for which two points this value is the greatest. These are the two vertices with the smallest angle triangle_vertex-centroid-Y, and will be the two points of the line the cloest point on the triangle will lie on.

Now find the closest point from Y to the line defined by these two points. It will necessarily lie in between the two points. Call this point Z.

Now return abs(P-Z)<abs(P-Y) ? Z : Y (This will return Y if it was already on the triangle and Z is further away).

This should be faster as it only calculates the closest point on a line once, and you skip the explicit check for if the point is in a triangle.

(norm(vector) is the normalized vector, or vector/abs(vector), abs(vector) is the magnitude of the vector, which seems to be vector.Length() in your code)

Artyer
  • 31,034
  • 3
  • 47
  • 75