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

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

Опираясь на предыдущее нехитрое упражнение с матрицей Хаусхолдера, попробуем реализовать примитивный алгоритм для трассировки лучей отражённых от неявно заданной поверхности второго порядка, и проведём численный эксперимент. Для этого мне понадобилось значительно модифицировать предыдущую программу (опубликованную ранее одним файлом на gist), и сделать из неё библиотеку -- раскидав декларацию и реализацию, дополнив нейминг, и снабдив сборку функциональным Makefile'ом.

Замечания по поводу конвенций библиотеки

Как это заведено в чистом Си, все типы снабжены именным префиксом библиотеки (sot — simple optics library). Я попытался, так же, употребить некоторое подобие венгерской нотации, из-за чего имена типов и функций приобрели несколько монструозный вид (вроде sot_Real_t вместо Real), однако я буду стараться держать имена переменных по-прежнему лаконичными для схожести с математической записью и сохранения читаемости (так некий целочисленный индекс $n$ будет называться n, а не in, как велит конвенция венгерской нотации).

Решение остаться на чистом Си хотя и полностью отвечает поговорке "работа любит дурака", продиктовано идеей, что этот счётный код затем можно будет в каком-то виде портировать на CUDA. Кроме того, меня давно занимала идея организовать обработку исключений в библиотеке на SJLJ, как это сделано в libpng — из чисто профессионального интереса.

Точка падения луча

Трассировка лучей в общем случае включает в себя задачу отыскания точки падения луча на поверхность в качестве первого шага. Для этого, располагая поверхностью $S(\vec{r})$ и лучом исходящим из точки $\vec{r_0}$ с заданным направлением $\vec{u}$, нужно попытаться отыскать точки пересечения прямой на которой лежит луч с этой поверхностью, что в математическом смысле означает их (прямой и поверхности) общие точки $r_1, r_2, ... r_n$. Объединив параметрическое уравнения прямой и уравнение поверхности в неявном виде, можно получить подстановку:

\begin{equation} \begin{cases} \vec{r} = \vec{r^0} + t \vec{u}, \\ S(\vec{r}) = 0 \end{cases} \label{eq:refCrossSubst} \end{equation}

Корни которого ($t_0, t_1, ... t_n$), подставленные в параметрическое уравнение прямой дадут координаты точек пересечения. Возможный случай отсутствия действительных корней, разумеется, соответствует прямой нигде не пересекающему поверхность. При этом, в силу геометрического смысла параметрического уравнения прямой, взяв за $\vec{r^0}$ точку исхода луча, корнями $t < 0$ можно пренебречь как находящимися "позади" ("за источником").

В то время, как для поверхностей порядков выше второго можно предусмотреть разрешение \eqref{eq:refCrossSubst} средствами того или иного численного метода, есть пока известный соблазн в том чтобы записать точное решение уравнения \eqref{eq:refCrossSubst}, ведь подстановка даёт в лучшем случае квадратный трёхчлен. Такие поверхности обладают известным рядом замечательных оптических свойств, и, как я читал, часто используются в практических инженерно-оптических решениях.

Аналитическое решение для поверхностей второго порядка

В общем случае уравнение поверхности второго порядка удобно записать через симметричную матрицу квадратичной формы $A_{ij}$ и коэффициенты линейной формы $l_{i}$:

\begin{equation} S_2(\vec{r}) = r_i A_{ij} r_j + 2 l_i r_i + d = 0, \end{equation}

где, согласно \eqref{eq:refCrossSubst}, $r_i$ выражаются как $r_i = r^0_i + u_i \cdot t$:

\begin{equation} (r^0_i + u_i \cdot t) A_{ij} (r^0_j + u_j \cdot t) + 2 l_i (r^0_i + u_i \cdot t) + d = 0. \end{equation}

Градиент при такой записи легко выражается через параметры линейной и квадратичной формы: \begin{equation} \nabla S_2 = \frac{\partial S_2(\vec{r})}{\partial r_i} = 2 A_{ij} r_j + 2 l_i, \end{equation} что уже позволяет говорить о вычислительной целесообразности формализации поверхностей второго порядка в виде отдельных структур в коде.

Если привести слагаемые при $t$, легко видеть, что подстановке \eqref{eq:refCrossSubst} удовлетворяют корни квадратного трёхчлена: $t_{1,2} = ( - b \pm \sqrt{ b^2 - 4 a c })/2 a$, где:

\begin{equation} \begin{aligned} a &= u_i A_{ij} u_j, \\ b &= 2 r^0_i A_{ij} u_j + 2 l_i u_j, \\ c &= r^0_i A_{ij} r^0_j + 2 l_i r^0_i + d. \end{aligned} \end{equation}

Примечателен тот алгебраический результат, что коэффициент $c$ выглядит в точности как подстановка вектора $r^0_i$ (отмечающего точку на прямой, которую мы условились выбирать точкой исхода луча) в выражение поверхности $S_2(\vec{r})$. Если точка $r^0_i$ лежит на поверхности, то выражение для корней упрощается до $t_1 = 0, t_2 = - b/a$, что соответствует пересечениям в установленной точке (простое следствие теоремы Виета), и какой-то удалённой от неё. Интересно подумать, есть ли смысле у того факта что $b = \nabla_i S_2(\vec{r^0}) \cdot u_i$?

Счёт и результаты

Прежде всего, добавим в код возможность задавать параметрические поверхности. Это коснётся прежде всего функции nabla_at() из предыдущей заметки, которая будет принимать более сложную структуру данных чем три указателя на функции соответствующие градиенту. Используя объявленные ранее типы данных, любую поверхность второго порядка можно определить через набор из квадратичной и линейной формы, а так же одного скалярного коэффициента:

struct {
    SymMatrix a;
    Vector3 l;
    Real_t d;
};

В действительности же, я сделал эту декларацию вложенной, для того чтобы поддерживать как более простой случай тройки статических функций задающих градиент, так и некоторой произвольной неявной функции задающей поверхность, которую я планирую обсчитывать целиком численными методами (однако пока неясно, как будет удобнее её задавать).

В выражениях для коэффициентов фигурирует произведение квадратичной формы $r_i A_{ij} r_j$, которую может быть эффективно записать в явном виде:

sot_Real_t
sot_quad_prod( const sot_SymMatrix_t * a, const sot_Vector3_t * v ) {
    return a->el[0] * v->x * v->x
         + a->el[1] * v->y * v->y
         + a->el[2] * v->z * v->z
         + 2*( a->el[3] * v->x * v->y
             + a->el[4] * v->y * v->z
             + a->el[5] * v->x * v->z )
         ;
}

Для несимметричного произведения выражение выглядит несколько сложнее:

sot_Real_t
sot_semi_quad_prod( const sot_Vector3_t * u
                  , const sot_SymMatrix_t * a
                  , const sot_Vector3_t * v ) {
    return (a->el[5] * u->r[2] + a->el[3] * u->r[1] + a->el[0]*u->r[0])*v->r[0]
         + (a->el[4] * u->r[2] + a->el[1] * u->r[1] + a->el[3]*u->r[0])*v->r[1]
         + (a->el[2] * u->r[2] + a->el[4] * u->r[1] + a->el[5]*u->r[0])*v->r[2]
         ;
}

Ну и, наконец, функция отыскания корней для решения уравнения на точки пересечения прямой заданной в параметрическом виде и произвольной поверхности второго порядка:

char
sot_quad_roots( const sot_SymMatrix_t * A, const sot_Vector3_t * l, const sot_Real_t d
          , const sot_Vector3_t * u, const sot_Vector3_t * r0
          , const sot_Real_t eps
          , sot_Real_t * t ) {
    const Real_t a = quad_prod( a, u )
               , b = 2*semi_quad_prod( r0, a, u ) + 2*prod( s->l, u )
               , c = quad_prod( a, r0 ) + 2*prod( s->l, r0 ) + s->d
               , D = b*b - 4*a*c
               ;
    if( D < 0. ) return -1;
    if( !a ) return -2; 
    if( D < eps ) {
        t[0] = - b/(2*a);
        return 1;
    }
    const sot_Real_t sqD = sqrt(D)
                   , a2 = 2*a
                   ;
    t[0] = - (b + sqD)/a2;
    t[1] = - (b - sqD)/a2;
    return 2;
}

Я не буду здесь приводить код для выражения затем точек пересечения $x_1, x_2$ из найденных корней, тестовый код и т.д. ввиду виду чрезвычайной тривиальности всей прочей числодробильни. Я сделал небольшую численную проверку чтобы убедиться в работоспособности метода и на более сложных поверхностях, но визуализацию сделал только на сфере: это и является довольно экстремальной проверкой, поскольку на точках экватора сферы возможно падение лучей практически под нулевым углом, а полюс должен давать луч отражённый в строго-обратном направлении.

Few rays on sphere	More rays

Как можно видеть на рисунке слева, отражение, в принципе, ведёт себя вполне ожидаемо и на экваторе, и на полюсах.

Пути дальнейшего развития

Ну, теперь можно либо дальше развивать функциональность счётного ядра, вводя

• обсчёт преломления лучей по закону Снеллиуса (Снелла),
• учёт диффузного рассеяния (касается ли оптических приборов?),
• поддержку сегментированной геометрии, ориентации в пространстве и т.д. — для описания, собственно, оптических приборов, в действительности состоящих из конечных линз и зеркал,
• поддержку вообще любой поверхности заданной функцией в неявном виде,

либо начать писать собственную визуализацию, поскольку для предстоящих задач функциональности Gnuplot'а начинает явно недоставать.