Отражение света векторами Numpy

[Vitalij Lysanov]

Фотореалистику выделили а отдельную программу, и это не случайно. Чтобы качественно отобразить освещенную поверхность нужны огромные затраты машинного времени. Сколько ни есть времени, можно все задействовать, особенно если вспомнить что лучи разного цвета по разному преломляются, а радуга вообще красота.


Немного про вектора

    import numpy as np
    from numpy.linalg import norm

>>> a=np.array([0,1,0])       (X,Y,Z)   Вертикальный вектор длиной 1
>>> b=np.array([1,0,0])                     Горизонтальный вектор длиной 1

>>> a+b
array([1, 1, 0])                                    Сложение векторов, (1,1) точка гипотенузы

>>> norm(a+b)               
1.4142135623730951                        Длина гипотенузы

>>> norm([1,1,1])
1.7320508075688772                        Длина вектора образованного единичными
                                                           векторами по осям


Точки x,y,z поверхности получены математически.
Позаимствовал на время часть поверхности модели:
http://forum.ascon.ru/index.php/topic,10204.msg223224.html#msg223224   Ответ №34


Для построения нормали к поверхности нужна следующая точка  по образующей поверхности, получаем математически. x2,y2,z2.
И точку вдоль другой образующей. x3,y3,z3.
По трем точкам строим вектора и находим нормаль к поверхности в каждой точке.     
Точки поверхности можно получить из формата stl. Вектора это две стороны треугольника.


            v2=np.array(    [(x2-x),(y2-y),(z2-z)]    )        Первый вектор
            v3=np.array(    [(x3-x),(y3-y),(z3-z)]    )       Второй вектор

            napr=np.cross(v2,v3)                             Нормаль это векторное произведение
            napr=edv(napr)                                      Сразу ставим единичную длину вектора


edv() -  функция вычисляет длину вектора и делит вектор на его длину

   def edv(tmp):
       dlina_tmp= norm(tmp, ord=2)
       ed_vektor=tmp/dlina_tmp
       return ed_vektor




Отраженный вектор считаем функцией otr(), рассчитывая равнобедренный треугольник.

def otr(t3,tn):

    t3_ed=edv(t3)
    tn_ed=edv(tn)

    dlina_do_tochki_normali=np.dot(t3,tn_ed)
    tp=dlina_do_tochki_normali*tn_ed

    T_otr=t3+2*(tp-t3)
    return T_otr



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

             T_otr_z=T_otr[2]                          Координата Z точки
            dz=50-z
            m=dz/T_otr_z                                Множитель
            t_vstr=[x,y,z]+m*T_otr                  Умножаем единичный отраженный вектор на
                                                                 множитель



Точки это как  свет лазерной указки. В каждой точке при диффузии света нужно нарисовать световое пятно с убывающей яркостью к краям пятна.

>>> cv
array([[0, 0, 1, 1, 1, 0, 0],         Массив для отображения пятна диффузии
       [0, 1, 1, 2, 1, 1, 0],
       [1, 1, 2, 3, 2, 1, 1],
       [1, 2, 3, 3, 3, 2, 1],
       [1, 1, 2, 3, 2, 1, 1],
       [0, 1, 1, 2, 1, 1, 0],
       [0, 0, 1, 1, 1, 0, 0]])
Размер массива 7*7.
Размер массива для накопления яркостей 100*100
Координаты отраженного луча дают смещение массива 7*7 для накопления яркостей.

k[pr:pr+7,ver:ver+7]=k[pr:pr+7,ver:ver+7]+cv
В массиве 100*100 берем срез чтобы выровнять размеры, тогда матрицы можно сложить.

Тут еще ошибка округления. Смещения pr и ver должны быть целыми. Можно и поточнее сложить, но медленней получится? 

          pr=50+int(t_vstr[0])
           ver=50+int(t_vstr[1])


   


Р.S.
____
Это не все. Это только начало.

При формате STL треугольники могут быть разного размера, а отраженный вектор один. Наверное нужно умножить на площадь треугольника. Как раз модуль векторного произведения векторов дает площадь.

И лучи заходит на плоскость 50 под разными углами, тоже нужно учесть.

И еще белый свет состоит из нескольких цветов.

И циклы нужно убрать, перевести на матрицы, тогда быстрее будет, и в динамике можно все отобразить
https://drive.google.com/file/d/1x-L5bCekFgXkohh_X1nPPQtJWu_nZuds/view?usp=sharing
Это тоже сложение волн только с учетом фазы волны

И ...

И зачем это нужно?,Колёсико в Компасе повернул и готово... 
Вспомнили школьные вектора, и хорошо,       может  пригодится.



___
PS PS

Вообще то видится задача не для теории, или нелинейной теории.
Допустим в Компасе получили нужное распределение яркости поверхности.
Это для параллельного пучка света.
И как форму этой поверхности преобразовать чтобы получить такой же пучек от точесного источника света.
Видится только проверять лучами.

Пластина в паралельном пучке равномерно освещена.
Точечный источник и парабола тоже дают параллельный пучек

Точечный и сфера соединяют в точку.
Парабола и паралельный пучек тоже соединяют в точку.

Закономерности?

Я про то что можно получить картину в Компасе и преобразовать в отражатель для точечного источника.
 

Код Выделить Развернуть

def edv(tmp):
    dlina_tmp= norm(tmp, ord=2)
    ed_vektor=tmp/dlina_tmp
    return ed_vektor


def otr(t3,tn):

    t3_ed=edv(t3)
    tn_ed=edv(tn)

    dlina_do_tochki_normali=np.dot(t3,tn_ed)
    tp=dlina_do_tochki_normali*tn_ed

    T_otr=t3+2*(tp-t3)
    return T_otr



if __name__ == "__main__":

    import numpy as np
    import matplotlib.cm as cm
    import matplotlib.pyplot as plt
    import math
    import cmath
    iDocument3D = iKompasObject.ActiveDocument3D()

    from numpy.linalg import norm

    import math

    no=0


    новый_деталь()                                      # Открываем новую деталь


    import math
    f2 = open("D:/tochki41.txt", "w")
    f3 = open("D:/50.txt", "w")



    cv=np.array([[0,0,1,1,1,0,0],[0,1,1,2,1,1,0],[1,1,2,3,2,1,1],[1,2,3,3,3,2,1],[1,1,2,3,2,1,1],[0,1,1,2,1,1,0],[0,0,1,1,1,0,0]])  #  Это массив для диффузии
# Тока встречи луча с плоскостью не целое число. Массив для накопления света от диффузии 100 на 100 и конечно целые числа. округлили. Можно точнее считать, но медленнее получится.



    cvet=np.array([5,5,33]    )
    точка_3д(cvet,6)

    k=np.zeros((100, 100), dtype=np.float)
#    k  Массив для диффузии


# Дальше идут циклы. желательно без циклов на матрицах посчитать.

    for nn in range (170,330,2):
        for i in range(0,150,2):


            ii=i/100.0
            n=nn*0.1

            x=n*n/486.0+n*math.sin(ii)
            f2.write(str(n*n/486.0+n*math.sin(ii)))
            f2.write("\t")

            y= n*math.cos(ii)
            f2.write(str(n*math.cos(ii)))
            f2.write("\t")


            z=n*1.168
            f2.write(str(n*1.168))
            f2.write("\n")



#____________________________________________

            ii=(i+1.0)/100.0
            n=n*1.0

            x2=n*n/486.0+n*math.sin(ii)
            y2= n*math.cos(ii)
            z2=n*1.168


#____________________________________________
            ii=(i)/100.0
            n=(n+1.0)*1.0

            x3=n*n/486.0+n*math.sin(ii)
            y3= n*math.cos(ii)
            z3=n*1.168
#______________________________________________



            v2=np.array(    [(x2-x),(y2-y),(z2-z)]    )
            v3=np.array(    [(x3-x),(y3-y),(z3-z)]    )

            napr=np.cross(v2,v3)
            napr=edv(napr)

            rast=[x,y,z]- cvet

            T_otr=otr(rast,napr)
#_________________________________________________________________
            T_otr_z=T_otr[2]




#________________________________________________________________________                 
            dz=50-z
            m=dz/T_otr_z
            t_vstr=[x,y,z]+m*T_otr


            pr=50+int(t_vstr[0])
            ver=50+int(t_vstr[1])



            k[pr:pr+7,ver:ver+7]=k[pr:pr+7,ver:ver+7]+cv

            f3.write(str(t_vstr[0]))
            f3.write("\t")
            f3.write(str(t_vstr[1]))
            f3.write("\t")
            f3.write(str(t_vstr[2]))
            f3.write("\n")





        f2.write("\n")
        f3.write("\n")

    f2.close()
    f3.close()


    im = plt.imshow(k, interpolation='bilinear', cmap=cm.RdYlGn,
                  origin='lower',
                  vmax=abs(k).max(), vmin=-abs(k).max())

    no=no+1
    plt.savefig("d:/1717/V"+str(no)+".png")






[Vitalij Lysanov]

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

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

За основу можно взять сами лучи света и для тренировки построить в 2D.
Для точного построения лучей нужно много, поэтому можно задать две линии в сечении пучка.  На линиях равномерно расположим точки и проведем лучи.

Также есть точка с источником света.

Берем последовательно два луча. От точечного источника света проводим отрезок к любой точке первого луча. Полученный угол делим пополам и в точке проводим нормаль к биссектриссе. Получили участок поверхности который от точечного источника света направит луч вдоль первой прямой.

Ищем точку встречи нормали м второго отрезка. Эта точка будет исходной для следующей пары лучей. И т.д. Лучей ложно взять хоть тысячу.


Полученная кривая профиля отражателя примерно аппроксимируется третьей степенью многочлена.



https://drive.google.com/open?id=1ukbWIgPjqwhCQO-9UwPIijbLcjOJkBo2aCjUaAszF0Q
Этот текст с картинками



   

Код Выделить Развернуть

f2 = open("D:/otr.txt", "w")


    toch=complex(20,-5)                        #   Точка источника света
    print (toch)
    точка(toch,5)


    shag=1

    n=0
    for y in range(1,120,shag):
        n=n+1
        t1=0+y*1j
        y2=100.0*math.sin(y/100.0)
        t2=100+y2*1j

        t3=0+(y+shag)*1j
        y2=100.0*math.sin((y+shag)/100.0)
        t4=100+y2*1j

# t1 t2     первый луч
# t3 t4     второй луч


        отрезок(t1,t2)
        отрезок(t3,t4)
        ug34=угол_отрезка(t3,t4)



        if n>1:
            tmp=pe                                # Для следующей пары берем координаты из
                                                        # предыдущего цикла
        else:
            tmp=t1+(t2-t1)*0.9



        точка(tmp,0)
        u1=угол_отрезка(tmp,t1)                    # Угол первого луча
        u2=угол_отрезка(tmp,toch)                # Угол первого луча

        u=0.5*math.pi+(u1+u2)/2.0                # Угол нормали к биссектриссе




        pr1=прямая(tmp,math.degrees(u))
        pr2=прямая(t3,math.degrees(ug34))

        pe=пересечение_прямых(pr1, pr2)
        точка(pe,1)

        f2.write(str(pe.real))
        f2.write("\t")

        f2.write(str(pe.imag))
        f2.write("\t")

        f2.write("\n")


    f2.write("\n")

    f2.close()

/code]



https://drive.google.com/open?id=1puM9SxDXwv02DadG2iV78mxSZFUWOZtVl2Au4_mwHK8
https://drive.google.com/open?id=165jHZZuPNVuarNWW76fVL1viZ37gr-Q1jvR0zzPT2Hs
Здесь подробней о примененных функциях Python.