Closest point projection of a 3D point to 3D triangles with numpy/scipy

Question

How would one calculate the closest projections of a point to N triangles using numpy/scipy?

Right now I would make a function to calculate a projection to a single triangle, basically this, then iterate across the entire array of triangles. But before I set off to do this I was wondering if there's already a solution built in scipy. Something like:

# DREAMY PSEUDOCODE
import numpy as np
N_TRIANGLES = 1000

point = np.random.rand(3) * 100 #random 3d point
triangles = np.random.rand(N_TRIANGLES,3,3) * 100 #array of triangles

from scipy.spatial import pointToTriangles
projections = pointToTriangles(point,triangles)

Here's a picture to help you visualize:

In the image above, the red dot in the middle is my query "point", the blue dots are the vertices of each triangles as define in the "triangles" np.array(). The green dots represent the result i want. They're the closest projections of "point" onto the defined triangles, and i wish to get this info back as an array of points.

cheers!

Could you clarify what you're asking for? In your example, what should projections contain? — ali_m
I still don't understand what you mean by projections of "point" to the triangle's surface. If I have a point and a triangle in 3D space then there are infinitely many lines I could draw connecting the point to the surface of the triangle. Could you maybe find/draw a diagram, or at least give an example with real numbers? — ali_m
@ali_m sorry, i'm looking for the closest projection to a point on triangles. — Fnord
You can project the point onto the plane of each triangle, is that what you want? Or do you also want to find the closest point in the triangle to that projected point? — Jaime
Well, projection typically happens along a well defined direction. If it's not the triangle's plane normal, you need to define which one to project along, or your problem is undefined. — Jaime

Fnord Fnord · Accepted Answer · 2015-09-11T18:05:32

Here's the code i came up with. I couldn't find anything in scipy that could help me directly, and this solution is about 2x faster than querying CGAL. It doesn't handle collapsed triangles, but that could be fixed by checking edge lengths and return the closest point on the longest edge instead.

import numpy as np
from numpy.core.umath_tests import inner1d

def pointsToTriangles(points,triangles):

    with np.errstate(all='ignore'):

        # Unpack triangle points
        p0,p1,p2 = np.asarray(triangles).swapaxes(0,1)

        # Calculate triangle edges
        e0 = p1-p0
        e1 = p2-p0
        a = inner1d(e0,e0)
        b = inner1d(e0,e1)
        c = inner1d(e1,e1)

        # Calculate determinant and denominator
        det = a*c - b*b
        invDet = 1. / det
        denom = a-2*b+c

        # Project to the edges
        p  = p0-points[:,np.newaxis]
        d = inner1d(e0,p)
        e = inner1d(e1,p)
        u = b*e - c*d
        v = b*d - a*e

        # Calculate numerators
        bd = b+d
        ce = c+e
        numer0 = (ce - bd) / denom
        numer1 = (c+e-b-d) / denom
        da = -d/a
        ec = -e/c


        # Vectorize test conditions
        m0 = u + v < det
        m1 = u < 0
        m2 = v < 0
        m3 = d < 0
        m4 = (a+d > b+e)
        m5 = ce > bd

        t0 =  m0 &  m1 &  m2 &  m3
        t1 =  m0 &  m1 &  m2 & ~m3
        t2 =  m0 &  m1 & ~m2
        t3 =  m0 & ~m1 &  m2
        t4 =  m0 & ~m1 & ~m2
        t5 = ~m0 &  m1 &  m5
        t6 = ~m0 &  m1 & ~m5
        t7 = ~m0 &  m2 &  m4
        t8 = ~m0 &  m2 & ~m4
        t9 = ~m0 & ~m1 & ~m2

        u = np.where(t0, np.clip(da, 0, 1), u)
        v = np.where(t0, 0, v)
        u = np.where(t1, 0, u)
        v = np.where(t1, 0, v)
        u = np.where(t2, 0, u)
        v = np.where(t2, np.clip(ec, 0, 1), v)
        u = np.where(t3, np.clip(da, 0, 1), u)
        v = np.where(t3, 0, v)
        u *= np.where(t4, invDet, 1)
        v *= np.where(t4, invDet, 1)
        u = np.where(t5, np.clip(numer0, 0, 1), u)
        v = np.where(t5, 1 - u, v)
        u = np.where(t6, 0, u)
        v = np.where(t6, 1, v)
        u = np.where(t7, np.clip(numer1, 0, 1), u)
        v = np.where(t7, 1-u, v)
        u = np.where(t8, 1, u)
        v = np.where(t8, 0, v)
        u = np.where(t9, np.clip(numer1, 0, 1), u)
        v = np.where(t9, 1-u, v)


        # Return closest points
        return (p0.T +  u[:, np.newaxis] * e0.T + v[:, np.newaxis] * e1.T).swapaxes(2,1)

Some test data projecting 100 points to 10k triangles:

    import numpy as np
    import cProfile

    N_TRIANGLES = 10**4 # 10k triangles
    N_POINTS    = 10**2 # 100 points
    points      = np.random.random((N_POINTS,3,)) * 100
    triangles   = np.random.random((N_TRIANGLES,3,3,)) * 100

    cProfile.run("pointsToTriangles(points,triangles)") # 54 function calls in 0.320 seconds

This will quickly become a memory hog, so when dealing with large data sets it's probably better to iterate over either points or triangles one at a time.

Closest point projection of a 3D point to 3D triangles with numpy/scipy

1 Answers