[Closed] Closest point on surface

But we are going into extreme deformations here, the original script was written to find the closest point to relatively regular surfaces, mostly for mapping purposes.

Well, I’m aware of that. But I wouldn’t miss a chance to learn from you, hehe.

Thanks a lot for your interest, guys!

Maybe we first can ignore the “efficient” bit and settle on an algorithm that is guaranteed to find the closest point, even for tricky surfaces?

What would be the theoretical way to solve this?

– MartinB

There can be many points that are closer to our point but whose normal does not hit our point, including the vertices.

Bobo, if I understand you correctly, then I disagree with you.
I think either the closest point is the orthogonal projection of ‘The Point’ on a face or one of the vertices. I mean that if it’s not one of the vertices, then it’s on a face and in this case the closest point is always the orthogonal projection. Looking at it in 2D with a broken line makes it easier to picture.

The fact that ‘The Point’ can be orthogonally projected on one or several faces (e.g. it is in the orthogonally extruded prism of one or several faces) doesn’t guarantee that the projected point is the closest. In the case of a concave surface, ‘The Point’ could be in the prism of a face in the hollow and still be closer to some vertices.

As I see it, without regard for the speed, the script should compute the distance and position of intersection points along normals if they exist and store the pos and dist in arrays. Then add the pos and dist of each vertex in the same arrays, no matter if they’re empty or not after the first part. Then look for the shortest distance, which can bring several equidistant points in some special cases.

Again, the only situation in which it doesn’t work is if there is a hole in the surface and ‘The Point’ is closer to a face on the other side (facing away) than to a point on the surface on this side of the model. And of course, if the point is inside the object, it doesn’t work either.
The only way around I can think of would be to copy the surface and flip the normals, which would simply double the time to compute. Unless there is a way to make the intersectRay() and intersectRayEx() methods work when the normal is facing the other way.

What do you think of this:

Input: Mesh M and query point Q
Output: Point P on M with minimal distance PQ

Pseudocode

0: mindist = MAXFLOAT; closestPoint=undefined
1: for each face F in the mesh M do (
2:   P := orthogonal projection of Q onto the plane that contains F
3:   If P is outside of F, then P := closest point on F to P
4:   d := distance between P and Q
5:   if d < mindist then ( mindist := d; closestPoint := P )
6: )
7: return closestPoint and mindist

Would that work?

I do not know yet how to do step #3 but since we are effectively in 2D (on the plane of the face) it should not be hard. Maybe there is a trick with barycentric coordinates?

Also it is unclear to me whether there is an opportunity to optimize/use some acceleration structure.

Any comments?
– MartinB

@Bobo: I apologize. It seems it came out wrong again. I didn’t mean to criticize your script or anything. I’m very far from having the expertise to do that. I didn’t know my words had such a strong meaning. Please don’t let what I said stop your from helping someone else.
And I was wrong anyway, I missed some other possibilities.

@Martin: in your code, it’s the third point that seems to be hard to figure correctly. So I asked a friend who teaches maths (to kids, but still). He changed the problem to something simpler: finding the closest point on a triangular surface (our faces, including vertices and edges) from a point in 3D space.

First, if the orthogonally projected point is in the triangle, that’s the closest.
This is the part done by the intersectRayEx() method in Bobo’s script.

Second, if the projection misses the triangle, then try an orthogonal projection on each of the three edges.
This is the part I missed. I don’t know if there is a method in MXS to do that. The implementation he explained to me uses trigonometry but there might be some shortcuts with the matrices.

Third, if the second projection misses the edges, then the closest point is one of the three points defining the triangle (our vertices).
This part is the easiest.

The guy knows nothing about programming or 3D modeling. When I told him how 3d models had their surface made of triangles, his opinion was that there were too many possible configurations for the surface and that the best way would be to do the above for each face, then pick the best result (or results if some points are equidistant).

The second part – the orthogonal projection on an edge – is probably the one that would slow down the process the most.
I’ll ask him the next time I see him – that would be wednesday – if there is a general method to find the closest point in a triangle ABC from a point P, no matter the projections. This might be much faster than going through three different processes.

What do you think?

The most-often-described method is pretty much as you quote here, but it starts by projecting the point orthogonally onto the (unbounded) plane defined by the triangle, and then testing the projected point against the extended edges of the triangle, in 2D, in the way you describe. Point>plane and Point>line are both in the sticky at the top of this forum.

These are pretty expensive operations, though, so you would ideally need to cull most of the triangles first, more cheaply, by using something like an octree division of the space enclosing the surface.

getting back to a render based solution is there any mileage in the near/far falloff map or even a z-depth pass to reduce the search parameters before raycasting?

I know how to project a point onto a line, but how do I test whether that point is on the edge of the triangle?

Interesting that you say that. I would expect a full ray-mesh intersection to be slower.

Sounds fairly good to me (and quite similar to what I suggested). I still wonder how I could project a point on a plane onto a triangle on the same plane without testing all the different possibilities explicitly. The idea would be to compute barycentric coordinates for that point and then somehow clamp these to the area of the triangle itself (0 <= barycentric coordinate <= 1) but I have not found a way to do that yet.

– MartinB

Hi Martin.

I’m actually trying to figure how to test the possibility of an orthogonal projection on a face before any calculation of the projected point, by only using distances (or better, squared distances, to avoid square roots).
The idea is that there should be a lot more faces on which the orthogonal projection is not possible than the other way around.

Testing for an edge is simple. Two verts: A and B, a point somewhere: P.
If abs(PAPA – PBPB) > AB*AB, then the orthogonal projection misses the edge [AB]. This is derived from the pythagorean theorem. (Correct me if I’m wrong, it’s getting late here.)
I’m trying to find something similar with a face instead of an edge.

Using squared distances calculated with the XYZ coordinates should be faster, as the length and distance methods calculate the square root. On the other hand, accessing the coordinates with the .pos property through MAXScript might take more time, while the dist and length methods might have some faster ‘internal’ access. I don’t know.
This needs to be tested.

What drdubosc said about projecting P on the plane (ABC) is a very good point. I don’t know if it’s faster to do it first or later. Quite basically, I’m counting additions, multiplications and variable accesses in the formulas. I’ll make some tests when I’m back home (I don’t have max on this computer).

Take a look at
http://www.blackpawn.com/texts/pointinpoly/default.html
the barycentric method. Even though it’s only concerned with finding if a point is inside a triangle, it’s easily extended to the 6 other zones outside it:

u<0 & v<0, closest to A,
u<0 & v<1, on AB,
u<0 & v>1, closest to B,
u<1 & v<0, on AC,
u>1 & v<0, closest to C,
u>0 & v>0 & u+v>1, on BC

furthermore, the dot products for finding the positions on AB, AC or BC if needed, have already been calculated in the process.

BTW, I’ve been having a fair amount of luck with MaxScript: MeshProjIntersect.

1st set up the data structure:

IPM = MeshProjIntersect()
IPM.SetNode $targetSurface
IPM.Build()

Then, having called IPM.ClosestFace [x,y,z], IPM.GetHitPos() returns the closest point on $targetSurface to [x,y,z] … it really does! Don’t know what the hitches are … I guess you have to be careful to call IPM.free()

I did a Bobo-style moving tape test, and one thing that becomes clear, is that (since none of these methods interpolate normals,) unless the point is relatively close to the ‘hard’ faceted surface, most of the closest points will be vertices.