How to get orthogonal distances of vectors from plane in Numpy/Scipy?

Question

I have a set of vectors as a numpy array. I need to get the orthogonal distances of each of those from a plane defined by 2 vectors v1 and v2. I can get this easily for a single vector using the Gram-Schmidt process. Is there a way to do it very fast in for many vectors without having to loop through each of them, or using np.vectorize?

Thanks!

Answer 1

You need to construct the unit-normal to the plane:

In three dimensions this can be easily done:

nhat=np.cross( v1, v2 )
nhat=nhat/np.sqrt( np.dot( nhat,nhat) )

and then dot this with each of your vectors; which I assume is an Nx3 matrix M

result=np.zeros( M.shape[0], dtype=M.dtype )
for idx in xrange( M.shape[0] ):
    result[idx]=np.abs(np.dot( nhat, M[idx,:] ))

so that result[idx] is the distance of the idx'th vector from the plane.

Answer 2

EDIT The original code I wrote does not work properly, so I have removed it. But following the same idea, explained below, if you spend some time thinking at it, there is no need for Cramer's rule, and the code can be streamlined quite a bit as follows :

def distance(v1, v2, u) :
    u = np.array(u, ndmin=2)
    v = np.vstack((v1, v2))
    vv = np.dot(v, v.T) # shape (2, 2)
    uv = np.dot(u, v.T) # shape (n ,2)
    ab = np.dot(np.linalg.inv(vv), uv.T) # shape(2, n)
    w = u - np.dot(ab.T, v)
    return np.sqrt(np.sum(w**2, axis=1)) # shape (n,)

To make sure it works properly, I have packed Dave's code into a function as distance_3d and tried the following:

>>> d, n = 3, 1000
>>> v1, v2, u = np.random.rand(d), np.random.rand(d), np.random.rand(n, d)
>>> np.testing.assert_almost_equal(distance_3d(v1, v2, u), distance(v1, v2, u))

But of course it now works for any d:

>>> d, n = 1000, 3
>>> v1, v2, u = np.random.rand(d), np.random.rand(d), np.random.rand(n, d)
>>> distance(v1, v2, u)
array([ 10.57891286,  10.89765779,  10.75935644])

---

You have to decompose your vector, lets call it u, in the sum of two vectors, u = v + w, v is in the plane, and so can be decomposed as v = a * v1 + b * v2, while w is perpendicular to the plane, and thus np.dot(w, v1) = np.dot(w, v2) = 0.

If you write u = a * v1 + b * v2 + w and take the dot product of this expression with v1 and v2, you get two equations with two unknowns:

np.dot(u, v1) = a * np.dot(v1, v1) + b * np.dot(v2, v1)
np.dot(u, v2) = a * np.dot(v1, v2) + b * np.dot(v2, v2)

Since it is only a 2x2 system, we can solve it using Cramer's rule as:

uv1 = np.dot(u, v1)
uv2 = np.dot(u, v2)
v11 = np.dot(v1, v2)
v22 = np.dot(v2, v2)
v12 = np.dot(v1, v2)
v21 = np.dot(v2, v1)
det = v11 * v22 - v21 * v12
a = (uv1 * v22 - v21 * uv2) / det
b = (v11 * uv2 - uv1 * v12) / det

From here, you can get:

w = u - v = u - a * v1 - b * v2

and the distance to the plane is the modulus of w.