Хабрахабр

Трехмерный движок на формулах Excel для чайников

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

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

Также я сделал небольшое веб-приложение, где можно потренироваться в создании формул для произвольных фигур и сгенерировать свой файл Excel.

Осторожно: 19 картинок и 3 анимации под катом.

В Excel-файлах вверху можно менять следующие ячейки:

  • alpha: горизонтальное вращение камеры в градусах,
  • beta: вертикальное вращение камеры в градусах,
  • dist: расстояние от камеры до начала координат,
  • fov: «зум» камеры.

Чтобы вращать камеру, необходимо дать Excel права на изменение файла.

Нужно выделить все маленькие квадратные ячейки и выбрать «Условное форматирование» → «Цветовые шкалы». В случае скачивания своего Excel-файла с ObservableHQ, ячейки необходимо раскрасить вручную.

Алгоритм Ray Marching был популяризован Iñigo Quilez в 2008 году после презентации «Rendering Worlds with Two Triangles» о компактном трехмерном движке для демосцены.

Оказалось, что при помощи нее можно создавать фотореалистичные сцены, посмотрите, например эту галерею. После этого Iñigo Quilez создал сайт Shadertoy, и большая часть отправляемых туда работ использует описанную методику.

Почему тогда он не используется на видеокартах? Ray Marching — это очень простой в реализации и эффективный алгоритм, позволяющий рендерить довольно сложные сцены. Дело в том, что видеокарты заточены под работу с трехмерными фигурами, состоящими из треугольников, для чего Ray Marching не так эффективен.

  • Объект описывается в виде функции.
  • Сопоставить каждому пикселю результата луч, выходящий из камеры.
  • Для каждого луча найти координаты пересечения с объектом.
  • По точке пересечения определить цвет пикселя.

Думаю, что в том или в ином виде все эти подзадачи решаются в любом трехмерном движке.
Можно придумать три способа описания трехмерного объекта:

  • Явная формула, которая по лучу, выходящему из камеры, возвращает точку пересечения с объектом (или, например, направление касательной в ней).
  • Функция принадлежности. По точке в пространстве сообщает, принадлежит она объекту или нет.
  • Вещественная функция. Для точек внутри объекта она отрицательная, снаружи — положительная.

В этой статье будет использоваться частный случай третьего варианта, называемый signed distance function (SDF).

Для точек внутри она равна расстоянию до границы объекта со знаком «минус» SDF получает на вход точку пространства и возвращает расстояние до ближайшей точки объекта.

Допустим, у нас есть луч, исходящий из камеры, и мы хотим найти точку его пересечения с объектом.

В голову приходит несколько способов найти эту точку:

  • Перебрать точки на луче с некоторым шагом.
  • Имея точки внутри объекта и снаружи, запустить двоичный поиск для уточнения места пересечения.
  • Либо запустить любой метод нахождения корней функции SDF (которая внутри объекта отрицательна, а снаружи положительна), ограниченной на луче.

Но здесь мы будем использовать другой метод: Ray Marching.

Это простой алгоритм, который работает только с SDF.

Вот алгоритм: Давайте параметризуем луч по расстоянию от его начала функцией $ray(t) = (x, y, z)$.
Тогда $ray(0)$ — это сама камера.

Здесь $N$ число итераций, которое мы можем себе позволить.
Чем больше $N$, тем точнее алгоритм.

Число $TMax$ это расстояние от камеры, в пределах которого мы ожидаем найти объект.

Всегда ли алгоритм дойдет до точки пересечения с объектом (точнее, будет стремиться к ней)? Объект может иметь сложную форму и другие ее части могут близко приближаться к лучу, не давая алгоритму сойтись к точке пересечения. На самом деле такого быть не может: другие части объекта обязательно будут на некотором ненулевом расстоянии от луча (иначе он бы их пересекал), которое мы обозначим $\delta>0$. Тогда функция SDF будет больше $\delta$ для любой точки луча (кроме точек, совсем близких к точке пересечения). Поэтому рано или поздно он приблизится к точке пересечения.

С другой стороны, если алгоритм сошелся к некоторой точке, то расстояние между соседними итерациями стремится к нулю $\Rightarrow$ SDF (расстояние до ближайшей точки объекта) тоже стремится к нулю $\Rightarrow$ в пределе $SDF = 0$$\Rightarrow$ алгоритм сходится к точке пересечения луча с объектом.
Представим, что для каждого пикселя мы нашли точку пересечения с объектом. Мы можем нарисовать эти значения (то есть расстояние от камеры до точки пересечения) напрямую и получить такие картинки:

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

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

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

Искомый угол $\varphi$ я обозначил одной дугой на чертеже. Значение его синуса обозначено как $k$.

Мне показалось, что это слишком много вычислений для Excel, и хотелось бы находить искомый угол проще. Обычно, когда используется метод Ray Marching, для нахождения нормали к объекту делают следующее: считают значение SDF в шести точках (несложно уменьшить число точек до четырех) рядом с точкой пересечения и находят градиент этой функции, который должен совпадать с нормалью.

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

Нельзя ли получить информацию об угле к поверхности исходя из скорости сходимости алгоритма? Если это так, что вообще никаких дополнительных вычислений не потребуется: достаточно проанализировать значения предыдущих итераций.

Обозначим $I_n$ и $I_$ — точки, полученные на соседних итерациях алгоритма. Нарисуем чертеж. Пусть $R = SDF(I_n)$, $r = SDF(I_{n+1})$ — расстояния от точек $I_n$ и $I_{n+1}$ до плоскости. Предположим, что точки находятся достаточно близко к объекту, что его можно заменить на плоскость, угол с которой мы хотим найти. Тогда, согласно алгоритму, расстояние между $I_n$ и $I_{n+1}$ равно $R$.

С другой стороны, если обозначить $X$ — точку пересечения луча и фигуры, то
$I_nX = R / k$, а $I_{n + 1}X = r / k$, следовательно

$R = I_nI_{n+1} = I_nX - I_{n+1}X = \frac{R - r}{k}$

$k = \frac{R - r}{R} = 1 - \frac{r}{R} = 1 - \frac{SDF(I_{n+1})}{SDF(I_{n})}$

То есть интенсивность пикселя равна единица минус отношение функций SDF соседних итераций!
Примеры SDF для сферы, куба и тора:
Фигуры можно перемещать и сжимать вдоль осей:
Фигуры также можно объединять, пересекать и вычитать. То есть SDF поддерживает основные операции конструктивной сплошной геометрии:
Внимательный читатель мог заметить, что у некоторых фигур (например, куба), а также для некоторых операций (например, сжатия), приведенные формулы не всегда возвращают расстояние до ближайшей точки объекта, то есть формально не являются SDF. Тем не менее, выяснилось, что их все равно можно подавать на вход движку SDF и он будет их корректно отображать.

Это малая часть того, что можно делать при помощи SDF, полный список фигур и операций можно найти на странице www.iquilezles.org/www/articles/distfunctions/distfunctions.htm

Совмещая все приведенные формулы, получаем первую фигуру:

$\min(\\ \quad\max( |x| - 0.3, |y| - 0.3, |z| - 0.3, -\sqrt{x^2 + y^2 + z^2} + 0.375 ),\\ \quad\sqrt{(\sqrt{(x - 0.25)^2 + (z - 0.25)^2} - 0.25)^2 + y^ 2} - 0.05\\ )$

С чайником немного сложнее: нужен план

Аккуратно записываем формулу:

27)^2 + z^2} - 0. $\min(\\ \quad \sqrt{x^2 + (y - 0. 5\cdot y^2 + z^2} - 0. 05\\ \quad \sqrt{x^2 + 2. 3)^2 + (y - 0. 4,\\ \quad \sqrt{(\sqrt{x^2 + z^2} - 0. 02,\\ \quad \max(\\ \qquad x + y - 0. 18)^2} - 0. 09,\\ \qquad \sqrt{(\sqrt{(x - 0. 7,\\ \qquad - y + 0. 09)^2) - 0. 55)^2 + (y - 0. 1)^2} - 0. 1}^2 + (z - 0. 09),\\ \qquad \sqrt{(\sqrt{(x - 0. 04,\\ \quad ),\\ \quad \max(\\ \qquad -(- y + 0. 09)^2} - 0. 35)^2 + (y - 0. 1)^2} - 0. 1)^2 + (z - 0. 04,\\ \quad ),\\ ) $

Готово: И подбираем нужный ракурс.

Все, что нам осталось — для данного пикселя на экране найти соответствующий ему луч в пространстве, выходящий из камеры. Если точнее, то надо уметь находить точку этого луча по расстоянию до камеры. То есть нам нужна функция $ray(s, t, d)$, где $s, t$ — координаты пикселя, а $d$ — расстояние от начала луча (камеры).

С учетом того, что экран состоит из $rows$ строк и $cols$ столбцов, мы ожидаем координаты в пределах $s\in\left[-\frac{cols}{2}, \frac{cols}{2}\right], t\in\left[-\frac{rows}{2}, \frac{rows}{2}\right]$ Для удобства вычислений координаты пикселей будут задаваться относительно центра экрана.

По сложности реализации они примерно одинаковы, но на практике чаще всего используется перспективная, поэтому мы будем использовать ее. Далее необходимо определиться с типом проекции камеры: ортогональная, перспективная или «рыбий глаз».

Но с учетом того, что фигуры также выровнены вдоль осей, в итоге получился бы не очень интересный ракурс. Большей части вычислений в этой главе можно было бы избежать, например, поместив камеру в точку $(1, 0, 0)$ и направив вдоль оси $X$. Так в случае куба мы бы видели его как квадрат.

Таким образом мы получаем три переменные на входе: угол $\alpha$, угол $\beta$ и расстояние $dist$. Чтобы обеспечить возможность поворота камеры, нужно аккуратно провести некоторое количество вычислений с Эйлеровыми углами. Они определяют как положение, так и направление камеры (камера всегда смотрит в начало координат).

При помощи WolframAlpha находим матрицу поворота:

$\begin{pmatrix}x\\y\\z\end{pmatrix} = \underbrace{\begin{pmatrix} \cos(\alpha) \cos(\beta) & -\cos(\alpha) \sin(\beta) & \sin(a) \\ \sin(\beta) & \cos(\beta) & 0 \\ -\cos(\beta) \sin(\alpha) & \sin(\alpha) \sin(\beta) & \cos(\alpha) \end{pmatrix}}_{M(\alpha, \beta)} \begin{pmatrix}x'\\y'\\z'\end{pmatrix}$

Если применить ее к вектору $(dist, 0, 0)$, получим координаты камеры (не спрашивайте меня, куда пропал минус):

$\begin{pmatrix} camX \\ camY \\ camZ \end{pmatrix} = M(\alpha, \beta) \begin{pmatrix}dist \\ 0 \\ 0 \end{pmatrix} = dist\begin{pmatrix} \cos(\alpha) \cdot \cos(\beta) \\ \sin(\beta)\\ \sin(\alpha) \cdot \cos(\beta) \end{pmatrix}$

Основной объект — экран (на рисунке он красного цвета, в тексте — курсивом). Последующие вычисления будут специфичны для перспективной проекции. Камера фактически просто является точкой с координатами $(camX, camY, camZ)$. Это воображаемый прямоугольник на некотором расстоянии перед камерой, имеющий, как можно догадаться, однозначное соответствие с пикселями обычного экрана. Лучи, соответствующие пикселям, начинаются в точке камеры и проходят через соответствующую точку экрана.

Точнее, они зависят от расстояния до камеры: если экран отодвинуть дальше, его потребуется сделать больше. У экрана нет точного местоположения и размера. Тогда мы можем вычислить значение вектора $(x0, y0, z0)$, соединяющего точку камеры и центр экрана (он такой же как и у центра камеры, только домноженный не на $dist$, а на $-1$): Поэтому условимся, что расстояние от камеры до экрана равно 1.

$\begin{pmatrix} x0\\y0\\z0 \end{pmatrix} = M(\alpha, \beta) \begin{pmatrix}-1\\ 0 \\ 0 \end{pmatrix} = -\begin{pmatrix} \cos(\alpha) \cdot \cos(\beta) \\ \sin(\beta)\\ \sin(\alpha) \cdot \cos(\beta) \end{pmatrix}$

Он зависит от угла обзора камеры, который измеряется в градусах и соответствует тому, что в видеокамерах называется «zoom». Теперь нам надо определить размер экрана. Так как экран не квадратный, необходимо уточнить, что имеется в виду вертикальный угол обзора. Пользователь задает его при помощи переменной $fov$ (field of view).

На основании этого можно определить размер одного пикселя на экране:
Итак, чтобы определить высоту экрана, нужно найти основание равнобедренного треугольника с вершинным углом $fov$ и высотой 1: вспоминая тригонометрию, получаем $2\tan\left(fov / 2\right)$.

$pixelSize = \frac{2\tan\left(fov / 2\right)}{rows}$

Далее, применяя матрицу поворота к векторам $(0, 0, 1)$ и $(0, 1, 0)$, получаем вектора $\overline{u}$ и $\overline{v}$, задающие горизонтальные и вертикальные направления экрана (для упрощения расчетов они также предварительно домножены на $pixelSize$):

$\begin{pmatrix} ux\\uy\\uz \end{pmatrix} = pixelSize\cdot M(\alpha, \beta) \begin{pmatrix}0\\ 0 \\ 1 \end{pmatrix} = pixelSize\begin{pmatrix} \sin(\alpha)\\ 0\\ \cos(\alpha) \end{pmatrix}$

$\begin{pmatrix} vx\\vy\\vz \end{pmatrix} = pixelSize\cdot M(\alpha, \beta) \begin{pmatrix}0\\ 1 \\ 0 \end{pmatrix} = pixelSize\begin{pmatrix} \cos(\alpha)\cdot\sin(\beta)\\ \cos(\beta)\\ \sin(\alpha)\cdot\cos(\beta) \end{pmatrix}$

Таким образом у нас теперь есть все компоненты, чтобы найти направление луча, выходящего из камеры и соответствующего пикселю с координатами $s, t$

$raydir(s, t) = \begin{pmatrix}x0\\y0\\z0\end{pmatrix} + s\begin{pmatrix}ux\\uy\\uz\end{pmatrix} + t \begin{pmatrix}vx\\vy\\vz\end{pmatrix}$

Только надо учесть, что начало луча находится в точке камеры и что вектор направления нужно нормировать:
Это почти то, что нужно.

$ray(s, t, d) = \begin{pmatrix}camX\\camY\\camZ\end{pmatrix} + d\cdot \frac{raydir(s, t)}{\lVert raydir(s, t)\rVert}$

Так мы получили искомую функцию $ray(s, t, d)$, которая возвращает точку луча на расстоянии $d$ от его начала, соответствую пикселю с координатами $s, t$.
Файл Excel, который получается в результате, является книгой, состоящей из >6 листов:
Все листы построены по одинаковой схеме:

Лист R

I1

rows:

50

V1

cols:

77

AI1

fov:

39

AV1

dist:

1,4

BI1

alpha:

35

BV1

beta:

20

**

=ЕСЛИ(i14!XN - i13!XN < 0,00000000000001; 0,09; МАКС(0; МИН(1; (i15!XN - i14!XN) / (i14!XN - i13!XN))))

Лист N

I1

pixelSize:

=TAN(R!AI1 / 2) / (R!I1 / 2)

**

=1 / КОРЕНЬ(СТЕПЕНЬ(X!AN + X!X2; 2) + СТЕПЕНЬ(Y!AN; 2) + СТЕПЕНЬ(Z!AN + Z!X2; 2))

Лист X

I1

camX:

=R!AV1 * COS(R!BV1) * COS(R!BI1)

V1

ux:

=-N!I1 * SIN(R!BI1)

AI1

vx:

=N!I1 * SIN(R!BV1) * COS(R!BI1)

AV1

x0:

=-COS(R!BV1) * COS(R!BI1)

A*

=AI1 * (СТРОКА() - 2 - (R!I1 + 1) / 2)

*2

=AV1 + V1 * (СТОЛБЕЦ() - 1 - (R!V1 + 1) / 2)

**

=(Z2 + AN) * N!ZN

Лист Y

I1

camY:

=R!AV1 * SIN(R!BV1)

V1

vy:

=-N!I1 * COS(R!BV1)

AI1

y0:

=-SIN(R!BV1)

A*

=AI1 + V1 * (СТРОКА() - 2 - (R!I1 + 1) / 2)

**

=AN * N!ZN

Лист Z

I1

camZ:

=R!AV1 * COS(R!BV1)) * SIN(R!BI1)

V1

uz:

= N!I1 * COS(R!BI1)

AI1

vz:

= N!I1 * SIN(R!BV1) * SIN(R!BI1)

AV1

z0:

=-COS(R!BV1) * SIN(R!BI1)

A*

=AI1 * (СТРОКА() - 2 - (R!I1 + 1) / 2)

*2

=AV1 + V1 * (СТОЛБЕЦ() - 1 - (R!V1 + 1) / 2)

**

=(Z2 + AN) * N!ZN

Лист i1

I1

dist0:

=formula(X!I1, Y!I1, Z!I1)

**

=formula(X!I1 + X!XN * I1, Y!I1 + Y!XN * I1, Z!I1 + Z!XN * I1)

Листы i2, i3, ...

**

=formula(X!I1 + X!XN * i*!XN, Y!I1 + Y!XN * i*!XN, Z!I1 + Z!XN * i*!XN)

Примечания:

  • Так как в Excel вычисления производятся в радианах, аргументы всех тригонометрических функций домножаются на $\frac{\pi}{180}$ (в Excel для этого есть функция RADIANS). Чтобы не запутывать формулы, я убрал эти умножения в таблицах выше.
  • Там, где написано formula, необходимо вставить одну из этих формул:

    Формула кубика с тором на Excel

    MIN( MAX( ABS(x) - 0.3, ABS(y) - 0.3, ABS(z) - 0.3, -SQRT(POWER(x, 2) + POWER(y, 2) + POWER(z, 2)) + 0.375 ), SQRT(POWER(SQRT(POWER(x - 0.25, 2) + POWER(z - 0.25, 2)) - 0.25, 2) + POWER(y, 2)) - 0.05
    )

    Формула чайника на Excel

    MIN( SQRT(POWER(SQRT(POWER(x, 2) + POWER(z, 2)) - 0.3, 2) + POWER(y - 0.18, 2)) - 0.02, SQRT(POWER(x, 2) + POWER(y, 2) * 2.5 + POWER(z, 2)) - 0.4, MAX( x + y - 0.15 - 0.05 - 0.5, - (y) + 0.19 - 0.1, SQRT(POWER(SQRT(POWER(x - 0.55, 2) + POWER(y - 0.09, 2)) - 0.1, 2) + POWER(z - 0.1, 2)) - 0.04 ), MAX( -(- (y) + 0.19 - 0.1), SQRT(POWER(SQRT(POWER(x - 0.35, 2) + POWER(y - 0.09, 2)) - 0.1, 2) + POWER(z - 0.1, 2)) - 0.04 ), SQRT(POWER(x, 2) + POWER(y - 0.27, 2) + POWER(z, 2)) - 0.05
    )

Как эта статья соотносится cо статьей «3D-движок, написанный на формулах MS Excel»?

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

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть