How do I get the closest point on mesh of given points? I want the location but not just the distance.

[HelloRicky123]

Thank you for your excellent work! In #193 it shows how to get the closest distance and the index of the triangle of a point. But it did not answer how to get the closest point. Like the function in blender python: new_object.closest_point_on_mesh.

Is there any advice on my title problem?

Thank you for your time!

[gkioxari]

The functionality in PyTorch3D that you are referring given a point p and a mesh M returns the index of the triangular face of M that is closes to p, let's call the closest face F. If I understand correctly, you want to find the projection of p on F, p_proj, which is the closest point of p on F. We don't provide p_proj as the output of our point_to_mesh distance functions but it's simple geometry.

Here is code to get it for a single point and the closest face tri which consists of the three vertices of face F.

def _point_to_tri_distance(point: torch.Tensor, tri: torch.Tensor) -> torch.Tensor:
  """
  Computes the squared euclidean distance of points to edges
  Args:
      point: FloatTensor of shape (3)
      tri: FloatTensor of shape (3, 3)
  Returns:
      dist: FloatTensor of shape (1)
  """
  a, b, c = tri.unbind(0)
  cross = torch.cross(b - a, c - a)
  norm = cross.norm()
  normal = torch.nn.functional.normalize(cross, dim=0)

  # p0 is the projection of p onto the plane spanned by (a, b, c)
  # p0 = p + tt * normal, s.t. (p0 - a) is orthogonal to normal
  # => tt = dot(a - p, n)
  tt = normal.dot(a) - normal.dot(point)
  p0 = point + tt * normal
  return p0

[HelloRicky123]

However, I think simply using normal to project the point onto the triangle is not enough. When the closest point is on the edges or the vertices, the normal will project the point out of the range of the triangle. I have an example code below.

import torch
import trimesh
import numpy as np

def _point_to_tri_distance(point: torch.Tensor, tri: torch.Tensor) -> torch.Tensor:
    """
    Computes the squared euclidean distance of points to edges
    Args:
      point: FloatTensor of shape (3)
      tri: FloatTensor of shape (3, 3)
    Returns:
      dist: FloatTensor of shape (1)
    """
    a, b, c = tri.unbind(0)
    cross = torch.cross(b - a, c - a)
    norm = cross.norm()
    normal = torch.nn.functional.normalize(cross, dim=0)

    # p0 is the projection of p onto the plane spanned by (a, b, c)
    # p0 = p + tt * normal, s.t. (p0 - a) is orthogonal to normal
    # => tt = dot(a - p, n)
    tt = normal.dot(a) - normal.dot(point)
    p0 = point + tt * normal
    return p0

vertices = np.array([[[0,0,0],[1,1,0],[-1,1,0]]],np.float32)
tris = torch.from_numpy(vertices).cuda()
points = np.array([[0,0.5,2],[0,0,2],[0,-0.5,2]])
pts = torch.from_numpy(points).float().cuda()
faces = np.array([[0,1,2]])

# mesh = trimesh.Trimesh(vertices, faces)

# scene = trimesh.Scene([
#     mesh,
#     trimesh.points.PointCloud(vertices),
#     trimesh.points.PointCloud(points),
# ])
# scene.show()
# (closest_points, distances, triangle_id) = mesh.nearest.on_surface(points)
# print(closest_points, distances, triangle_id)
# print(points)
# [[0.  0.5 0. ]
#  [0.  0.  0. ]
#  [0.  0.  0. ]] [1.         1.         1.11803399] [0 0 0]
# [[ 0.   0.5  1. ]
#  [ 0.   0.   1. ]
#  [ 0.  -0.5  1. ]]

import pytorch3d
from pytorch3d import _C

nrPts = len(points)
startIndex = torch.LongTensor([0]).cuda()
endIndex = torch.LongTensor([nrPts]).cuda()

dists, idxs = _C.point_face_dist_forward(pts, startIndex, tris, startIndex, endIndex)
print(dists, idxs)
print(_point_to_tri_distance(pts[2], tris[0]))
# should be [0.  0.  0. ] as supposed, however get tensor([ 0.0000, -0.5000,  0.0000], device='cuda:0')

[gkioxari]

The code I gave you will project the point on the plane spanned by the vertices of the closest face. What you describe is absolutely correct and reasoning about whether the point is inside or outside is what we do to compute the distances (and yes it is outside the triangle then we compute distances to edges). If you want to follow the full algorithm and computations we do for point to mesh distances check out our source code and the tests.

In particular, this test function might be of use

https://github.com/facebookresearch/pytorch3d/blob/fc4dd80208bbcf6f834e7e1594db0e106961fd60/tests/test_point_mesh_distance.py#L172

[gaoyuankidult]

It seems that finding the projected point on a triangle is not that trivial. One may check this thread for an answer: https://stackoverflow.com/questions/32342620/closest-point-projection-of-a-3d-point-to-3d-triangles-with-numpy-scipy Another idea: https://gdbooks.gitbooks.io/3dcollisions/content/Chapter4/closest_point_to_triangle.html

[SuwoongHeo]

@gaoyuankidult For someone who is interested, I had implemented the method in PyTorch. By referring to this cool article, it might be a bit different from the code in that of the StackOverflow link tho. I got ~0.25sec elapsed time for 256*30 query points. I attach my implementation below. Note that mesh is not the class of pytorch3d.

pmkeys = ['a','b','c','h', 'delta'] # values that can be precomputed
            # values that should be computed on the fly
            # 'd','e','f','sbar','tbar','g','k','bd','ce','be','ad']
class ProjectMesh2Point(object):
    def __init__(self, mesh :Mesh):
        # Note. edge is oriented in right-handed direction,
        self.E0 = mesh.edges[:,0].squeeze()
        self.E1 = -mesh.edges[:,1].squeeze()
        self.B = mesh.triangles[:,0].squeeze()
        for key in pmkeys:
            setattr(self, key, [])

        # Values that can be precomputed
        self.a = T.sum(self.E0 * self.E0, dim=-1) # inner product
        self.b = T.sum(self.E0 * self.E1, dim=-1)
        self.c = T.sum(self.E1 * self.E1, dim=-1)
        self.delta = self.a * self.c - self.b**2
        self.h = self.a - 2*self.b + self.c

    def __call__(self, query, index):
        """
        The projected point is computed by
        Pout = B + sout*E0 + tout*E1
        query : Px3 vector
        index : Px1 indices for updateing vars
        """
        # Update projection variables according to input queries
        # variables of "d,e,f,sbar,tbar,g,k,bd,ce,be,ad" are updated

        # Get vars in indices
        B, E0, E1 = self.B[index,:], self.E0[index,:], self.E1[index,:]
        a, b, c, h = self.a[index], self.b[index], self.c[index], self.h[index]
        delta = self.delta[index]
        D = B - query
        d = T.sum(E0*D, dim=-1)
        e = T.sum(E1*D, dim=-1)
        # f = T.sum(D*D, dim=-1)

        sbar = b * e - c * d
        tbar = b * d - a * e

        bd = b + d
        ce = c + e

        ad = a + d
        be = b + e

        g = ce - bd
        k = ad - be

        # output
        sout = T.zeros_like(sbar, device=sbar.device)
        tout = T.zeros_like(tbar, device=tbar.device)

        # Region classification
        r_conds = T.stack([(sbar+tbar)<=delta, sbar>=0., tbar>=0., bd>ce, d<0, be>ad])

        # Inside triangle
        r_0 = r_conds[0] & r_conds[1] & r_conds[2] # region 0
        sout[r_0] = sbar[r_0]/delta[r_0]
        tout[r_0] = tbar[r_0]/delta[r_0]

        # region 1
        r_1 = ~r_conds[0] & r_conds[1] & r_conds[2]
        sout[r_1] = T.clip(g[r_1]/h[r_1], 0., 1.)
        tout[r_1] = 1 - sout[r_1]

        # region 2
        r_2 = ~r_conds[0] & ~r_conds[1] & r_conds[2]
        r_2a = r_2 & r_conds[3] # region 2-a
        sout[r_2a] = T.clip(g[r_2a] / h[r_2a], 0., 1.)
        tout[r_2a] = 1 - sout[r_2a]
        r_2b = r_2 & ~r_conds[3] # region 2-b
        tout[r_2b] = T.clip(-e[r_2b]/c[r_2b], 0., 1.) # Note. sout=0 in r_2b

        # region 3
        r_3 = r_conds[0] & ~r_conds[1] & r_conds[2]
        tout[r_3] = T.clip(-e[r_3] / c[r_3], 0., 1.)  # Note. sout=0 in r_3

        # region 4
        r_4 = r_conds[0] & ~r_conds[1] & ~r_conds[2]
        r_4a = r_4 & r_conds[4] # region 4-a
        sout[r_4a] = T.clip(-d[r_4a]/a[r_4a], 0., 1.) # Note tout=0 in r_4a
        r_4b = r_4 & ~r_conds[4] # region 4-b
        tout[r_4b] = T.clip(-e[r_4b] / c[r_4b], 0., 1.)  # Note sout=0 in r_4b

        # region 5
        r_5 = r_conds[0] & r_conds[1] & ~r_conds[2]
        sout[r_5] = T.clip(-d[r_5]/a[r_5], 0., 1.) # Note tout=0 in r_5

        # region 6
        r_6 = ~r_conds[0] & r_conds[1] & ~r_conds[2]
        r_6a = r_6 & r_conds[5]
        tout[r_6a] = T.clip(k[r_6a]/h[r_6a], 0., 1.)
        sout[r_6a] = 1 - tout[r_6a]
        r_6b = r_6 & ~r_conds[5]
        tout[r_6b] = T.clip(-d[r_6b]/a[r_6b], 0., 1.) # Note sout=0 in r_6b

        # Sanity check
        # Should be false all
        # r_1 & r_2 & r_3 & r_4 & r_5 & r_6
        # (r_2a & r_2b) | (r_4a & r_4b) | (r_6a & r_6b)

        Pout = B + sout[...,None]*E0 + tout[...,None]*E1

        return Pout, (sout, tout)

[gaoyuankidult]

@SuwoongHeo Thanks for your information. Do you think it is faster than this idea ? The PyTorch3D implementation also used it in the following function. https://github.com/facebookresearch/pytorch3d/blob/fc4dd80208bbcf6f834e7e1594db0e106961fd60/tests/test_point_mesh_distance.py#L172 I think the p0 in the function is the closest point on the triangle.

It seems to be easier than the idea mentioned in the article

[SuwoongHeo]

@gaoyuankidult I'm not good at C# so it is hard to guess. But I think for that idea, it may depend on how fast the ClosetPoint that I cannot find in the article is. And for Pytorch3D implementation, it seems not to return where the closest point is but distance value. My implementation is intended to get texture coordinate information. For getting point to surface distance only, it might be better to use Pytorch3D's implementation.

+) Since this article provides a closed-form solution. I guess it may faster than the others. When it comes to parallel implementation. I`m not sure again :(

[JohnYeung-dojjy]

I have made an implementation here following the discussion at 193 and this thread.

The _point_to_tri_distance() method provided by @gkioxari indeed gives the coordinates of the projected point onto the closest face on the mesh, although the docstring hinted otherwise. I have made the comments clearer in my repo.