Как добраться ортогональных векторов расстояния от самолета в пакете numpy/составляющей?

У меня есть набор векторов в виде массива numpy. Мне нужно получить ортогональные расстояния каждого из них от плоскости, определяемой 2 векторами v1 и v2. Я могу легко получить это для одного вектора, используя процесс грамма-Шмидта. Есть ли способ сделать это очень быстро для многих векторов, не проходя через каждый из них, или используя np?векторизовать?



Спасибо!

611   3  

3 ответов:

Более явный способ получить ответ @Jaime, вы можете быть явным относительно построения оператора проекции: 

def normalized(v):
    return v/np.sqrt( np.dot(v,v))
def ortho_proj_vec(v, uhat ):
    '''Returns the projection of the vector v  on the subspace
    orthogonal to uhat (which must be a unit vector) by subtracting off
    the appropriate multiple of uhat.
    i.e. dot( retval, uhat )==0
    '''
    return v-np.dot(v,uhat)*uhat

def ortho_proj_array( Varray, uhat ):
     ''' Compute the orhogonal projection for an entire array of vectors.
     @arg Varray:  is an array of vectors, each row is one vector
          (i.e. Varray.shape[1]==len(uhat)).
     @arg uhat: a unit vector
     @retval : an array (same shape as Varray), where each vector
               has had the component parallel to uhat removed.
               postcondition: np.dot( retval[i,:], uhat) ==0.0
               for all i. 
    ''' 
    return Varray-np.outer( np.dot( Varray, uhat), uhat )




# We need to ensure that the vectors defining the subspace are unit norm
v1hat=normalized( v1 )

# now to deal with v2, we need to project it into the subspace
# orhogonal to v1, and normalize it
v2hat=normalized( ortho_proj(v2, v1hat ) )
# TODO: if np.dot( normalized(v2), v1hat) ) is close to 1.0, we probably
# have an ill-conditioned system (the vectors are close to parallel)



# Act on each of your data vectors with the projection matrix,
# take the norm of the resulting vector.
result=np.zeros( M.shape[0], dtype=M.dtype )
for idx in xrange( M.shape[0] ):
    tmp=ortho_proj_vec( ortho_proj_vec(M[idx,:], v1hat), v2hat )             
    result[idx]=np.sqrt(np.dot(tmp,tmp))

 # The preceeding loop could be avoided via
 #tmp=orhto_proj_array( ortho_proj_array( M, v1hat), v2hat )
 #result=np.sum( tmp**2, axis=-1 )
 # but this results in many copies of matrices that are the same 
 # size as M, so, initially, I prefer the loop approach just on
 # a memory usage basis.
Это действительно всего лишь обобщение процедуры ортогонализации Грама-Шмидта. заметим, что в конце этого процесса мы имеем np.dot(v1hat, v1hat.T)==1, np.dot(v2hat,v2hat.T)==1, np.dot(v1hat, v2hat)==0 (в пределах числовой точности)
4  
2013-02-05 16:03:03

Вам нужно построить единицу-нормаль к плоскости:

В трех измерениях это можно легко сделать: 

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

И затем расставьте точки над каждым из ваших векторов; я предполагаю, что это матрица Nx3 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,:] ))

Так что result[idx] - это расстояние вектора idx'th от плоскости.

4  
2013-02-01 18:48:06

Исходный код, который я написал, не работает должным образом, поэтому я удалил его. Но следуя той же идее, объясненной ниже, если вы потратите некоторое время на обдумывание, нет необходимости в правиле Крамера, и код можно немного упростить следующим образом: 

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,)
Чтобы убедиться, что он работает правильно, я упаковал код Дейва в функцию как distance_3d и попробовал следующее: 
>>> 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))

Но, конечно, теперь это работает для любого 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])

Ты придется разложить ваш вектор, назовем его u, в сумму двух векторов, u = v + w, v находится в плоскости, и поэтому может быть разложен как v = a * v1 + b * v2, в то время как w перпендикулярен плоскости, и таким образом np.dot(w, v1) = np.dot(w, v2) = 0.

Если написать u = a * v1 + b * v2 + w и взять точечное произведение этого выражения с v1 и v2, то получим два уравнения с двумя неизвестными: 

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)
Поскольку это всего лишь система 2x2, мы можем решить ее, используя правило Крамера Как: 
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

Отсюда вы можете получить: 

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

А расстояние до плоскости-это модуль w.

2  
2013-02-01 19:18:40

Comments

    Ничего не найдено.
Click here to cancel reply.