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

В этой статье я расскажу, как мне удалось портировать алгоритм рендера трехмерных сцен на формулы 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 будет больше <img src=$ для любой точки луча (кроме точек, совсем близких к точке пересечения). Поэтому рано или поздно он приблизится к точке пересечения.

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

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

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

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

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

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

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

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

Нарисуем чертеж. Обозначим $I_n$ и $I_{n+1}$ — точки, полученные на соседних итерациях алгоритма. Предположим, что точки находятся достаточно близко к объекту, что его можно заменить на плоскость, угол с которой мы хотим найти. Пусть $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\\ )$

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

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

$\min(\\ \quad \sqrt{x^2 + (y - 0.27)^2 + z^2} - 0.05\\ \quad \sqrt{x^2 + 2.5\cdot y^2 + z^2} - 0.4,\\ \quad \sqrt{(\sqrt{x^2 + z^2} - 0.3)^2 + (y - 0.18)^2} - 0.02,\\ \quad \max(\\ \qquad x + y - 0.7,\\ \qquad - y + 0.09,\\ \qquad \sqrt{(\sqrt{(x - 0.55)^2 + (y - 0.09)^2) - 0.1}^2 + (z - 0.1)^2} - 0.04,\\ \quad ),\\ \quad \max(\\ \qquad -(- y + 0.09),\\ \qquad \sqrt{(\sqrt{(x - 0.35)^2 + (y - 0.09)^2} - 0.1)^2 + (z - 0.1)^2} - 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)$ . Лучи, соответствующие пикселям, начинаются в точке камеры и проходят через соответствующую точку экрана.

У экрана нет точного местоположения и размера. Точнее, они зависят от расстояния до камеры: если экран отодвинуть дальше, его потребуется сделать больше. Поэтому условимся, что расстояние от камеры до экрана равно 1. Тогда мы можем вычислить значение вектора $(x0, y0, z0)$ , соединяющего точку камеры и центр экрана (он такой же как и у центра камеры, только домноженный не на $dist$ , а на $-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;12 МАКС(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»?
Тот движок подходит только для рендера лабиринтов и эллипсоидов. Используется одномерный ray tracer, из данных которого вверх и вниз рисуются полосы, создающие иллюзию стен. Но зато там реализована полноценная игровая логика, здесь же целью является только трехмерный движок.