Трехмерный движок внутри запроса SQL / Хабр

2019-01-11 19:18:13

Трехмерный движок внутри запроса SQL

Несколько лет назад на форуме SQL.ru решили провести сравнение реализаций трассировщиков лучей на разных языках программирования. К сожалению, моя заявка не может участвовать т.к. она не выводит надпись «PIXAR», поэтому публикую ее здесь.

Для чистоты эксперимента я использовал SQLite без расширений. Оказалось, что там нет даже функции SQRT.

WITH RECURSIVE numbers AS (SELECT 0 AS n UNION ALL SELECT n+1 FROM numbers WHERE n<89),
pixels AS (SELECT rows.n as row, cols.n as col FROM numbers as rows CROSS JOIN
numbers as cols WHERE rows.n > 4 AND rows.n < 38 AND cols.n > 9 AND cols.n < 89),
rawRays AS (SELECT row, col, -0.9049 + col * 0.0065 + row * 0.0057 as x,
-0.1487 + row * -0.0171 as y, 0.6713 + col * 0.0045 + row * -0.0081 as z FROM pixels),
norms AS (SELECT row, col, x, y, z, (1 + x * x + y * y + z * z) / 2 as n FROM rawRays),
rays AS (SELECT row, col, x / n AS x, y / n AS y, z / n AS z FROM norms),
iters AS (SELECT row, col, 0 as it, 0 as v FROM rays UNION ALL
SELECT rays.row, rays.col, it + 1 AS it, v + MAX(ABS(0.7+v*x) - 0.3,
ABS(0.7+v*y) - 0.3, ABS(-1.1+v*z) - 0.3, -((0.7+v*x) * (0.7+v*x) +
(0.7+v*y) * (0.7+v*y) + (-1.1+v*z) * (-1.1+v*z)) * 1.78 + 0.28) AS v
FROM iters JOIN rays ON rays.row = iters.row AND rays.col = iters.col WHERE it < 15),
lastIters AS (SELECT it0.row, it0.col, it0.v AS v0, it1.v AS v1, it2.v AS v2
FROM iters as it0 JOIN iters AS it1 ON it0.col = it1.col AND it0.row = it1.row
JOIN iters AS it2 ON it0.col = it2.col AND it0.row = it2.row
WHERE it0.it = 15 AND it1.it = 14 AND it2.it = 13),
res AS (SELECT col, (v0 - v1) / (v1 - v2) as v FROM lastIters)
SELECT group_concat(
substr('$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/|()1{}[]?-_+~<>i!lI;:,"^. ',
round(1 + max(0, min(66, v * 67))), 1) || CASE WHEN col=88 THEN X'0A' ELSE '' END, '')
FROM res;



Здесь можно покрутить кубик

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


Disclaimer: я далек от мира БД, поэтому буду раз замечаниям в личку.

Версия для Postgres (UPD: благодаря флагам стало работать на порядок быстрее, UPD2: еще ряд улучшений, теперь время выполнения 150мс)

Для понимания терминологии и принципа работы алгоритма рекомендуется ознакомиться со статьей про ray marching для Excel.

Общая структура


Список промежуточных таблиц:



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

Все таблицы (за исключением numbers и iters) содержат по 81 x 29 строк (по одной на каждый «пиксель»), колонки row и col индексируют их координаты. Таблица iters содержит 81 x 29 x 15 строк (по одной на каждую итерацию ray marching для каждого «пикселя»). Номер итерации содержится в колонке it.

Финальный запрос выдает таблицу из одной строки и колонки с текстом, все остальные таблицы содержат только вещественные числа (при этом колонки row, col и it – целые неотрицательные).

Получается, если опустить вспомогательные таблицы, очень простая структура запроса:

WITH RECURSIVE
 numbers AS (SELECT ...),
 pixels AS (SELECT ...),
 rawRays AS (SELECT ...),
 normsSq AS (SELECT ...),
 norms AS (SELECT ...),
 rays AS (SELECT ...),
 iters AS (SELECT ...),
 lastIters AS (SELECT ...),
 res AS (SELECT ...)

SELECT group_concat(substr('$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/|()1{}[]?-_+~<>i!lI;:,"^. ', round(1 + max(0, min(66, v * 67))), 1) || CASE WHEN col=88 THEN X'0A' ELSE '' END, '') FROM res;

Рекурсивные запросы


Вот стандартный способ получить таблицу, содержащую числа от 0 до 89:

WITH RECURSIVE numbers AS (
 SELECT 0 AS n
 UNION ALL
 SELECT n+1
  FROM numbers
  WHERE n<89
) ...


Извлекаем квадратный корень


Мы будем использовать метод Герона (вавилонский метод) вычисления корня. Допустим, мы хотим вычислить $\sqrt{S}$. Мы строим ряд $x_0, x_1, \ldots$ по следующей формуле:

$x_{n+1} = \frac{x_n + \frac{S}{x_n}}{2}$



Логика метода очень простая: $\sqrt{S}$ всегда лежит между $x$ и $\frac{S}{x}$ для любого числа $x$. Поэтому естественно взять середину отрезка между этими числами как приближение.

Геометрически это можно изобразить так:



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

Начальное значение $x_0$ может быть любым положительным числом, например 1. В игре Quake для этого использовалась магическая константа 0x5f3759df (точнее, она использовалась для инвертированного квадратного корня, тем не менее аналогичный метод может быть придуман и для обычного корня), но, к сожалению, в SQL нет доступа к двоичному представлению чисел с плавающей запятой.

В этой статье корень нужен в двух местах:

  • при нормировании векторов, выходящих из камеры: ray marching сильно зависит от расстояний, а для того, чтобы их откладывать, нужен вектор длины 1.
  • при вычислении расстояния до границы сферы, которая вырезается из квадрата.


В первом случае значения находятся в узком диапазоне $[0.7, 1.5]$ и начальное приближение 1 идеально подходит. Во втором случае, собрав статистику по вызовам, получил среднее значение $0.28$, которое было взято в качестве начального.

Оказалось, что при правильном начальном значении достаточно одной итерации! То есть в нашем случае корень приближается линейной функцией:

$sqrt_1(x) = \frac{1 + x}{2}$


$sqrt_2(x) = \frac{0.28 + \frac{x}{0.28}}{2} = 0.14 + 1.78 x$



Вычисляем лучи из камеры


Задача первых четырех таблиц — каждому «пикселю» сопоставить трехмерный вектор длины 1, выходящий из камеры и проходящий через соответствующую точку экрана.

Сначала необходимо получить таблицу с нужной структурой, то есть с ячейками, для которых указаны номер строки и номер столбца. Для этого берется декартово произведение набора чисел от 0 до 89 и из него вырезаются пустые строки и столбцы:

...
pixels AS (
 SELECT rows.n as row, cols.n as col
 FROM numbers as rows
  CROSS JOIN numbers as cols
 WHERE rows.n >= 5 AND rows.n < 38 AND cols.n >= 10 AND cols.n < 89
),
...

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

...
rawRays AS (
 SELECT
  row, col,
  -0.9049 + col * 0.0065 + row * 0.0057 as x,
  -0.1487 + row * -0.0171 as y,
  0.6713 + col * 0.0045 + row * -0.0081 as z
 FROM pixels
),
...

После этого мы должны посчитать (приблизительные) длины этих векторов по формуле $sqrt_1$:

...
norms AS (
 SELECT
  row, col, x, y, z,
  (1 + x * x + y * y + z * z) / 2.0 AS n
 FROM rawRays
),
...

Осталось разделить координаты векторов на их длину для получения векторов длины 1:

...
rays AS (SELECT row, col, x / n AS x, y / n AS y, z / n AS z FROM norms),
...

Итерации ray marching


Здесь используется чуть более сложная конструкция с рекурсивными запросами, содержащая JOIN. Мы хотим произвести 15 итераций алгоритма ray marching для каждого пикселя. Если при рекурсивном вычислении таблицы numbers каждый раз таблица содержала по одной строке, которые потом объединялись, здесь промежуточные таблицы содержат по 81 x 29 строк и вычисляются 15 раз.

Вся трехмерная геометрия содержится в формуле

$SDF\begin{pmatrix}x\\ y \\ z\end{pmatrix} = \max\left(|x| - 0.3, |y| - 0.3, |z| - 0.3, -\left(sqrt_2(x^2 + y^2 + z^2) - 0.42\right)\right)$


  • функция $\max$ означает пересечение
  • $|x| - 0.3, |y| - 0.3, |z| - 0.3$ задают три пары полуплоскостей, образующих куб со стороной $0.3 \cdot 2$
  • $-\left(sqrt_2(x^2 + y^2 + z^2) - 0.42\right)$ — наружная часть сферы радиуса $0.42$. Радиус взят больше видимого, чтобы компенсировать неточность приближения квадратного корня.


Далее нам надо просто вычислить последовательность $0 = v_0, v_1, v_2 \ldots, v_{15}$ для каждого пикселя по формуле:

$v_{n+1} = v_n + SDF\left( \begin{pmatrix}camX\\camY\\camZ\end{pmatrix} + v_n \begin{pmatrix}x\\y\\z\end{pmatrix} \right)$


где $x, y, z$ — координаты нормированного вектора. Так как координаты камеры повторяются несколько раз, я округлил их до одного десятичного знака.

...
iters AS (
 SELECT
  row, col,
  0 as it,
  0 as v
 FROM rays
 UNION ALL
 SELECT
  rays.row,
  rays.col,
  it + 1 AS it,
  v + MAX(
   ABS(0.7+v*x) - 0.3,
   ABS(0.7+v*y) - 0.3,
   ABS(-1.1+v*z) - 0.3,
   -(
    (0.7+v*x) * (0.7+v*x) + (0.7+v*y) * (0.7+v*y) + (-1.1+v*z) * (-1.1+v*z)
   ) * 1.78 + 0.28
  ) AS v
 FROM iters
  JOIN rays
   ON rays.row = iters.row AND rays.col = iters.col
 WHERE it < 15
),
...

Получаем интенсивности «пикселей»


Здесь используется та же формула, что и в Excel, которая аппроксимирует компоненту diffuse из затенения по Фонгу:

$intensity = \max\left(0, \min\left(1, \frac{v_{15} - v_{14}}{v_{14} - v_{13}}\right)\right) $



Чтобы ее вычислить, необходимо предварительно сделать таблицу с тремя последними итерациями ray marching:

...
lastIters AS (
 SELECT
  it0.row,
  it0.col,
  it0.v AS v0,
  it1.v AS v1,
  it2.v AS v2
 FROM iters as it0
  JOIN iters AS it1
   ON it0.col = it1.col AND it0.row = it1.row
  JOIN iters AS it2
   ON it0.col = it2.col AND it0.row = it2.row
 WHERE it0.it = 15 AND it1.it = 14 AND it2.it = 13
),
...

И, собственно, сама формула (операции $\max$ и $\min$ будут применены в финальном запросе):

...
res AS (SELECT row, col, (v0 - v1) / (v1 - v2) as v FROM lastIters)
...

Генерируем «ascii-арт»


...
SELECT group_concat(
 substr(
  '$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/|()1{}[]?-_+~<>i!lI;:,"^. ',
  round(1 + max(0, min(66, v * 67))),
  1
 ) ||
 CASE
  WHEN col=88 THEN X'0A'
  ELSE ''
 END, 
'') FROM res;

Задача финального запроса состоит в том, чтобы конвертировать таблицу с интенсивностями пикселей в одну ascii-строку. На вход он получает только таблицу res, содержащую колонки col и v.

  • group_concat(s, delim) – агрегирующая функция, конкатенирующая выражение s для всех строк, используя строку delim в качестве разделителя.
  • CASE WHEN cond1 THEN val1 WHEN cond2 THEN val2 ... ELSE valN END – условная конструкция, аналог тернарного оператора.
  • X'0A' – символ переноса строки, который вставляется перед первым символом каждой строки.
  • || – оператор конкатенации строк.
  • substr(s, start, count) – функция, возвращающая count символов строки s, начиная с символа с номером start. Индексация символов идет с единицы.
  • '$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/|()1{}[]?-_+~<>i!lI;:,"^. ' – строка, содержащая «градиент» от «черного» ($) к «белому» (пробел) в ascii-символах. Взято с сайта http://paulbourke.net/dataformats/asciiart/.
  • round(1 + max(0, min(66, v * 67))) – преобразуем вещественные числа из интервала $[0, 1]$ в целое число в интервале $[1, 67]$ чтобы взять символ с соответствующим номером.
Взято отсюда