Как использовать np.where в другом np.where (conext: трассировка лучей)

Question

Вопрос в том, как использовать два np.where в одном и том же утверждении, например так (упрощенно):

np.where((ndarr1==ndarr2),np.where((ndarr1+ndarr2==ndarr3),True,False),False)

Чтобы избежать вычисления второго условного выражения, если первое не достигнуто.

Моя первая задача - найти пересечение луча в треугольнике, если он есть.Эта проблема может быть решена с помощью этого алгоритма (находится в stackoverflow):

def intersect_line_triangle(q1,q2,p1,p2,p3):
    def signed_tetra_volume(a,b,c,d):
        return np.sign(np.dot(np.cross(b-a,c-a),d-a)/6.0)

    s1 = signed_tetra_volume(q1,p1,p2,p3)
    s2 = signed_tetra_volume(q2,p1,p2,p3)

    if s1 != s2:
        s3 = signed_tetra_volume(q1,q2,p1,p2)
        s4 = signed_tetra_volume(q1,q2,p2,p3)
        s5 = signed_tetra_volume(q1,q2,p3,p1)
        if s3 == s4 and s4 == s5:
           n = np.cross(p2-p1,p3-p1)
           t = np.dot(p1-q1,n) / np.dot(q2-q1,n)
           return q1 + t * (q2-q1)
    return None

Вот два условных оператора:

1. s1! = S2
2. s3 == s4 & s4 == s5

Теперь, когда у меня есть> 20k треугольников для проверки, я хочу применить эту функцию ко всем треугольникам одновременно.

Первое решение:

s1 = vol(r0,tri[:,0,:],tri[:,1,:],tri[:,2,:])
s2 = vol(r1,tri[:,0,:],tri[:,1,:],tri[:,2,:])

s3 = vol(r1,r2,tri[:,0,:],tri[:,1,:])
s4 = vol(r1,r2,tri[:,1,:],tri[:,2,:])
s5 = vol(r1,r2,tri[:,2,:],tri[:,0,:])

np.where((s1!=s2) & (s3+s4==s4+s5),intersect(),False)

, где s1, s2, s3, s4, s5 - массивы, содержащие значение S для каждого треугольника.Проблема в том, что это означает, что мне нужно вычислить s3, s4 и s5 для всех треугольников.

Теперь идеальным было бы вычислить утверждение 2 (и s3, s4, s5) только тогда, когда утверждение 1 истинно, счто-то вроде этого:

check= np.where((s1!=s2),np.where((compute(s3)==compute(s4)) & (compute(s4)==compute(s5), compute(intersection),False),False)

(чтобы упростить объяснение, я просто указал «вычислять» вместо всего вычислительного процесса. Здесь «вычислять» - это только для соответствующих треугольников).

Теперь, конечно, эта опция не работает (и вычисляет s4 два раза), но я бы с удовольствием дал несколько рекомендаций по аналогичному процессу

Answer 1

Я могу получить небольшое ускорение от короткого замыкания, но я не уверен, что это стоит дополнительного администратора.

full computation 4.463818839867599 ms per iteration (one ray, 20,000 triangles)
short ciruciting 3.0060838296776637 ms per iteration (one ray, 20,000 triangles)

Код:

import numpy as np

def ilt_cut(q1,q2,p1,p2,p3):
    qm = (q1+q2)/2
    qd = qm-q2
    p12 = p1-p2
    aux = np.cross(qd,q2-p2)
    s3 = np.einsum("ij,ij->i",aux,p12)
    s4 = np.einsum("ij,ij->i",aux,p2-p3)
    ge = (s3>=0)&(s4>=0)
    le = (s3<=0)&(s4<=0)
    keep = np.flatnonzero(ge|le)
    aux = p1[keep]
    qpm1 = qm-aux
    p31 = p3[keep]-aux
    s5 = np.einsum("ij,ij->i",np.cross(qpm1,p31),qd)
    ge = ge[keep]&(s5>=0)
    le = le[keep]&(s5<=0)
    flt = np.flatnonzero(ge|le)
    keep = keep[flt]
    n = np.cross(p31[flt], p12[keep])
    s12 = np.einsum("ij,ij->i",n,qpm1[flt])
    flt = np.abs(s12) <= np.abs(s3[keep]+s4[keep]+s5[flt])
    return keep[flt],qm-(s12[flt]/np.einsum("ij,ij->i",qd,n[flt]))[:,None]*qd

def ilt_full(q1,q2,p1,p2,p3):
    qm = (q1+q2)/2
    qd = qm-q2
    p12 = p1-p2
    qpm1 = qm-p1
    p31 = p3-p1
    aux = np.cross(qd,q2-p2)
    s3 = np.einsum("ij,ij->i",aux,p12)
    s4 = np.einsum("ij,ij->i",aux,p2-p3)
    s5 = np.einsum("ij,ij->i",np.cross(qpm1,p31),qd)
    n = np.cross(p31, p12)
    s12 = np.einsum("ij,ij->i",n,qpm1)
    ge = (s3>=0)&(s4>=0)&(s5>=0)
    le = (s3<=0)&(s4<=0)&(s5<=0)
    keep = np.flatnonzero((np.abs(s12) <= np.abs(s3+s4+s5)) & (ge|le))
    return keep,qm-(s12[keep]/np.einsum("ij,ij->i",qd,n[keep]))[:,None]*qd

tri = np.random.uniform(1, 10, (20_000, 3, 3))
p0, p1 = np.random.uniform(1, 10, (2, 3))

from timeit import timeit
A,B,C = tri.transpose(1,0,2)
print('full computation', timeit(lambda: ilt_full(p0[None], p1[None], A, B, C), number=100)*10, 'ms per iteration (one ray, 20,000 triangles)')
print('short ciruciting', timeit(lambda: ilt_cut(p0[None], p1[None], A, B, C), number=100)*10, 'ms per iteration (one ray, 20,000 triangles)')

Обратите внимание, что я немного поиграл с этим алгоритмом, так что это может не во всех случаях получить такой же результат, как у вас.

Что я изменил:

• Я добавил тетра объем, что позволяет сохранить несколько повторных подвычислений
• Я заменяю один из концов луча на среднюю точку M луча. Это позволяет сэкономить на расчете одного тетраобъема (s1 или s2), поскольку можно проверить, пересекает ли луч треугольную плоскость ABC, сравнив объем тетра ABCM с суммой s3, s4, * 1020. * (если они имеют одинаковые знаки).

Answer 2

Вот как я использовал замаскированные массивы для решения этой проблемы:

    loTrue= np.where((s1!=s2),False,True)
    s3=ma.masked_array(np.sign(dot(np.cross(r0r1, r0t0), r0t1)),mask=loTrue)
    s4=ma.masked_array(np.sign(dot(np.cross(r0r1, r0t1), r0t2)),mask=loTrue)
    s5=ma.masked_array(np.sign(dot(np.cross(r0r1, r0t2), r0t0)),mask=loTrue)
    loTrue= ma.masked_array(np.where((abs(s3-s4)<1e-4) & ( abs(s5-s4)<1e-4),True,False),mask=loTrue)

    #also works when computing s3,s4 and s5 inside loTrue, like this:        
    loTrue= np.where((s1!=s2),False,True)
    loTrue= ma.masked_array(np.where(
            (abs(np.sign(dot(np.cross(r0r1, r0t0), r0t1))-np.sign(dot(np.cross(r0r1, r0t1), r0t2)))<1e-4) &
            (abs(np.sign(dot(np.cross(r0r1, r0t2), r0t0))-np.sign(dot(np.cross(r0r1, r0t1), r0t2)))<1e-4),True,False)
            ,mask=loTrue)

Обратите внимание, что тот же процесс, когда не используется такой подход, выполняется следующим образом:

    s3= np.sign(dot(np.cross(r0r1, r0t0), r0t1)  /6.0)
    s4= np.sign(dot(np.cross(r0r1, r0t1), r0t2)  /6.0)
    s5= np.sign(dot(np.cross(r0r1, r0t2), r0t0)  /6.0)
    loTrue= np.where((s1!=s2) & (abs(s3-s4)<1e-4) & ( abs(s5-s4)<1e-4) ,True,False)

ОбаТем не менее, вы можете получить те же результаты, когда зацикливаетесь на этом процессе только для 10 000 итераций, НЕ используя маскированные массивы быстрее!(26 секунд без замаскированных массивов, 31 секунда с замаскированными массивами, 33 с использованием замаскированных массивов только в одну строку (без вычисления s3, s4 и s5 отдельно или до вычисления s4 до этого).

Вывод: использование вложенных массивовздесь решено (обратите внимание, что маска указывает, где она не будет вычислена, поэтому сначала loTri должен поставить значение False (0), когда условие проверено). Однако в этом сценарии это не быстрее.