За последние два месяца я написал несколько маленьких GLSL-демо. О первом из них, Red Alp, я написал статью. В ней я подробно расписал весь процесс, поэтому рекомендую прочитать её, если вам незнакома эта сфера.

Мы рассмотрим четыре демо: Moonlight, Entrance 3, Archipelago и Cutie. Но на этот раз я расскажу лишь о паре уроков, которые извлёк из каждого. Мы не будем углубляться во все аспекты, потому что это было бы излишне.

Moonlight

Демо Moonlight состоит из 460 символов

// Moonlight [460] by bµg
// License: CC BY-NC-SA 4.0
void main(){vec3 o,p,u=vec3((P+P-R)/R.y,1),Q;Q++;for(float d,a,m,i,t;i++<1e2;p=t<7.2?Q:vec3(2,1,0),d=abs(d)*.15+.1,o+=p/m+(t>9.?d=9.,Q:p/d),t+=min(m,d))for(p=normalize(u)*t,p.z-=5e1,m=max(length(p)-1e1,.01),p.z+=T,d=5.-length(p.xy*=mat2(cos(t*.2+vec4(0,33,11,0)))),a=.01;a<1.;a+=a)p.xz*=mat2(8,6,-6,8)*.1,d-=abs(dot(sin(p/a*.6-T*.3),p-p+a)),m+=abs(dot(sin(p/a/5.),p-p+a/5.));o/=4e2;O=vec4(tanh(mix(vec3(-35,-15,8),vec3(118,95,60),o-o*length(u.xy*.5))*.01),1);}

Его можно посмотреть на официальной странице или поэкспериментировать с кодом в порте в Shadertoy.

Для движения через облака и туман в Red Alp я использовал volumetric raymarching, и достаточно большая часть кода потребовалась для того, чтобы поглощение и излучение выглядели убедительно. Но есть и альтернативная методика, которая на удивление проста.

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

Объяснение 1/d

Давайте посмотрим, как это работает на практике, начав с raymarch простого куба, в котором мы применим этот принцип:

void main() {
    float d, t;
    vec3 o, p,
         u = normalize(vec3(P+P-R,R.y)); // координаты screen to world

    for (int i = 0; i < 30; i++) {
        p = u * t; // позиция луча

        p.z -= 3.; // делаем шаг назад

        // Поворот Родрига с произвольно выбранным углом π/2
        // и невыровненной осью
        vec3 a = normalize(cos(T+vec3(0,2,4)));
        p = a*dot(a,p)-cross(a,p);

        // Функция расстояний со знаком куба размером 1
        p = abs(p)-1.;
        d = length(max(p,0.)) + min(max(p.x,max(p.y,p.z)),0.);

        // Берём максимум, чтобы не войти в сам куб
        d = max(d,.001);

        t += d; // делаем шаги вперёд на это расстояние

        // Наш загадочный вклад цвета в вывод
        o += 1./d;
    }

    // Произвольное масштабирование в видимом диапазоне
    O = vec4(o/200., 1);
}

Функция расстояний со знаком для куба взята с классической страницы Иниго Килеза. Информацию о повороте можно изучить в статье Xor или Blackle. Чтобы получить общее представление о коде, прочитайте мою статью про Red Alp.

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

Давайте упростим задачу при помощи следующего рисунка:

Луч, проходящий рядом с излучающим объектом

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

Допустим, мы испустили луч и хотим знать, сколько фотонов встречается на протяжении всего его пути. Мы должны «суммировать», а точнее интегрировать величины плотностей всех этих фотонов вдоль луча. Так как мы выполняем дискретное сэмплирование (точки на рисунке), нам необходимо интерполировать плотность фотонов и между точками сэмплирования.

Пусть у нас есть две произвольные точки сэмплирования с расстояниями и $d_{n+1}$ ; тогда любое промежуточное расстояние можно линейно интерполировать, как $r=\mathrm{mix}(d_n,d_{n+1},t)$ , где находится в интервале . Применив закон обратных квадратов, получим, что интегрированная плотность фотонов между этими двумя точками может быть выражена следующей формулой (в интервале ):

$v = \Delta t \int \frac{1}{\mathrm{mix}(d_n,d_{n+1},t)^2} dt$

нормализуется, поэтому iобозначает настоящее расстояние между отрезками. При помощи Sympy мы можем выполнить интегрирование:

>>> a, b, D, t = symbols('a b D t', real=True)
>>> mix = a*(1-t) + b*t
>>> D * integrate(1/mix**2, (t,0,1)).simplify()
 D
───
a⋅b

То есть результатом этого интегрирования будет

$v = \frac{\Delta t}{d_{n}d_{n+1}}$

В цикле шаг — это $d_{n+1}$ , то есть мы получаем:

$v = \frac{\Delta t}{d_{n}\Delta t} = \frac{1}{d_n}$

И так мы получаем наше загадочное . Оно «физически корректно», если считать, что свет распространяется в вакууме. Разумеется, в реальности всё сложнее, и мы даже не обязаны придерживаться этой формулы, но было здорово понять, что эта простая дробь оказывается достаточно хорошей моделью реальности.

Прохождение сквозь объект

В примере с кубом свет не проходит через объект, потому что мы используем max(d, .001). Но если бы мы хотели добавить прозрачности, то можно было бы использовать d = A*abs(d)+B, где A можно интерпретировать, как поглощение, а B — как прохождение, или прозрачность.

Впервые упоминание этой формулы я встретил в статье Xor о volumetric. Вот, как понимаю её я: +B приводит к потенциальному проникновению в твёрдое тело на следующей итерации, которое в противном случае бы не произошло (или произошло бы на малую глубину). При нахождении внутри тела abs(d) приводит к тому, что луч движется дальше (на величину расстояния до ближайшей грани). Умножение на A гарантирует, что мы не проникнем в тело слишком быстро; это поглощение, или «затухание».

По сути, именно эту методику я и использовал в Moonlight, чтобы не писать сложный код поглощения/испускания.

Entrance 3

Демо Entrance 3 состоит из 465 символов

// Entrance 3 [465] by bµg
// License: CC BY-NC-SA 4.0
#define V for(s++;d<l&&s>.001;q=abs(p+=v*s)-45.,b=abs(p+vec3(mod(T*5.,80.)-7.,45.+sin(T*10.)*.2,12))-vec3(1,7,1),d+=s=min(max(p.y,-min(max(abs(p.y+28.)-17.,abs(p.z+12.)-4.),max(q.x,max(q.y,q.z)))),max(b.x,max(b.y,b.z))))
void main(){float d,s,r=1.7,l=2e2;vec3 b,v=b-.58,q,p=mat3(r,0,-r,-1,2,-1,b+1.4)*vec3((P+P-R)/R.y*20.4,30);V;r=exp(-d*d/1e4)*.2;l=length(v=-vec3(90,30,10)-p);v/=l;d=1.;V;r+=50.*d/l/l;O=vec4(pow(mix(vec3(0,4,9),vec3(80,7,2),r*r)*.01,p-p+.45),1);}

Посмотреть его можно на его официальной странице, а поэкспериментировать с кодом — в порте на Shadertoy .

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

Изначально я решил попробовать воксели, но не смог уместить работу со светом в менее чем 512 символов (код инициализации был слишком большим, это не пошаговый DDA без ветвления). Также возникли неприятные ограничения, поэтому я вернулся к классическому raymarching.

Первым, что отличается в моём коде, стало использование нормы L-∞ вместо евклидовой нормы для функции расстояния: все твёрдые тела — это кубы, поэтому допустимо применять более простые формулы.

Свет — это не иллюзия, а настоящий свет: после первого raymarch в твёрдое тело направление луча переориентируется в сторону источника света, и raymarching выполняется снова (это макрос V). Попадание в твёрдое тело определяет, должен ли фрагмент освещаться.

Баги на мобильных платформах

Неприятной неожиданностью для меня стало обнаружение двух багов драйверов на мобильных платформах:

Один со сложными составными циклами for на Snapdragon/Adreno, потому что я всеми силами пытался избежать макросов и функций.
Второй с цепочными присваиваниями на Imagination/PowerVR (обычно касается Google Pixel Pro 10).

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

Изометрия

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

# Поворот против часовой стрелки
a, ax0, ax1, ax2 = symbols('a ax0:3')
c, s = cos(a), sin(a)
k = 1-c
rot = Matrix(3,3, [
    # столбец 1        столбец 2          # столбец 3
    ax0*ax0*k + c,     ax0*ax1*k + ax2*s, ax0*ax2*k - ax1*s, # строка 1
    ax1*ax0*k - ax2*s, ax1*ax1*k + c,     ax1*ax2*k + ax0*s, # строка 2
    ax2*ax0*k + ax1*s, ax2*ax1*k - ax0*s, ax2*ax2*k + c      # строка 3
])

# Поворот на 45° по оси Y
m45 = rot.subs({a:rad(-45), ax0:0, ax1:1, ax2:0})

# Применяем второй поворот по оси X, чтобы получить матрицы преобразований
# для двух классических проекций
# Примечание: asin(tan(rad(30))) эквивалентно atan(sin(rad(45)))
isometric = m45 * rot.subs({a:asin(tan(rad(30))), ax0:1, ax1:0, ax2:0})
dimetric  = m45 * rot.subs({a:         rad(30),   ax0:1, ax1:0, ax2:0})

Изучив матрицы и вынеся общие члены, мы получим следующие матрицы преобразований:

$M_{iso} = \sqrt{2}\sqrt{3}\begin{bmatrix} \sqrt{3} & -1 & \sqrt{2} \\ 0 & 2 & \sqrt{2} \\ -\sqrt{3} & -1 & \sqrt{2} \end{bmatrix} \text{ и } M_{dim} = \frac{4}{\sqrt{2}}\begin{bmatrix} 2 & -1 & \sqrt{3} \\ 0 & \sqrt{6} & \sqrt{2} \\ -2 & -1 & \sqrt{3} \end{bmatrix}$

Направление луча общее для всех фрагментов, поэтому мы используем в качестве опорной точки центральную координату UV (0,0). Для удобства перемещаем её вперёд: (0,0,1) и преобразуем при помощи матрицы. Так мы получаем центральную экранную координату в пространстве мира. Поскольку полученная координата точки относительна к точке начала координат мира, чтобы перейти от этой точки к точке начала координат, нужно поменять её знак. Формула направления луча принимает следующий вид:

$d_{iso} = -M_{iso} \begin{bmatrix}0 \\ 0 \\ 1\end{bmatrix} = -\frac{\sqrt{3}}{3}\begin{bmatrix}1 \\ 1 \\ 1\end{bmatrix} \text{ и } d_{dim} = -M_{dim} \begin{bmatrix}0 \\ 0 \\ 1\end{bmatrix} = -\frac{1}{4} \begin{bmatrix}\sqrt{6} \\ 2 \\ \sqrt{6}\end{bmatrix}$

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

Это требование можно смоделировать при помощи выражения

$M \begin{bmatrix}x \\ -1 \\ z\end{bmatrix} > 0$

Где -1 — самая нижняя экранная координата Y (которая не должна оказаться в «полу»). Я лентяй, поэтому просто попросил Sympy решить его за меня:

x, z = symbols("x z", real=True)
u = m * Matrix([x, -1, z])
uz = solve(u[1] > 0, z)

Мы получаем z>2 для изометрии и для диметрии.

При произвольном масштабе S координаты мы получим следующее:

const float S = 50.;
vec2 u = (P+P-R)/R.y * S; // масштабированные экранные координаты

float A=sqrt(2.), B=sqrt(3.);

// Изометрия
vec3 rd = -vec3(1)/B,
     ro = mat3(B,0,-B,-1,2,-1,A,A,A)/A/B * vec3(u, A*S + eps);

// Диметрия
vec3 rd = -vec3(B,A,B)/A/2.,
     ro = mat3(2,0,-2,-1,A*B,-1,B,A,B)/A/2. * vec3(u, B*S + eps);

eps — это произвольное малое значение, гарантирующее, что координата Y окажется выше 0.

В Entrance 3 я использовал грубую аппроксимацию изометрической системы.

Archipelago

Демо Archipelago состоит из 472 символов

// Archipelago [472] by bµg
// License: CC BY-NC-SA 4.0
#define r(a)*=mat2(cos(a+vec4(0,11,33,0))),
void main(){vec3 p,q,k;for(float w,x,a,b,i,t,h,e=.1,d=e,z=.001;i++<50.&&d>z;h+=k.y,w=h-d,t+=d=min(d,h)*.8,O=vec4((w>z?k.zxx*e:k.zyz/20.)+i/1e2+max(1.-abs(w/e),z),1))for(p=normalize(vec3(P+P-R,R.y))*t,p.zy r(1.)p.z+=T+T,p.x+=sin(w=T*.4)*2.,p.xy r(cos(w)*e)d=p.y+=4.,h=d-2.3+abs(p.x*.2),q=p,k-=k,a=e,b=.8;a>z;a*=.8,b*=.5)q.xz r(.6)p.xz r(.6)k.y+=abs(dot(sin(q.xz*.4/b),R-R+b)),k.x+=w=a*exp(sin(x=p.x/a*e+T+T)),p.x-=w*cos(x),d-=w;}

Посмотреть его можно на его официальной странице, а поэкспериментировать с кодом — в порте на Shadertoy.

Этой бесконечной процедурно генерируемой Японией я хотел ознаменовать свой разрыв с одержимостью красным и оранжевым. Строго говоря, на самом деле она довольно проста. Я использовал тот же шум из Red Alp для гор/островов, но для воды применён другой шум.

Кривая октавного шума имеет вид w=exp(sin(x)), где сдвиг координаты x выполняется при помощи её производной : x-=w*cos(x). Это некая разновидность domain warping, создающая здесь красивый эффект. Под x я на самом деле подразумеваю позицию по оси X. Для работы с компонентой Z это не нужно (xz образует плоскую поверхность), потому что каждая октава фрактального броуновского движения имеет поворот, «смешивающий» обе оси, поэтому z на самом деле копируется в x.

Примечание: я не придумал эту формулу, а нашёл её в видео Acerola. Не знаю, она ли её автор, но видел, что эту формулу повторяли во многих местах.

Cutie

Демо Cutie состоит из 602 символов

// Cutie [602] by bµg
// License: CC BY-NC-SA 4.0
#define V vec3
#define L length(p
#define C(A,B,X,Y)d=min(d,-.2*log2(exp2(X-L-A)/.2)+exp2(Y-L-B)/.2)))
#define H(Z)S,k=fract(T*1.5+s),a=V(1.3,.2,Z),b=V(1,.3*max(1.-abs(3.*k-1.),z),Z*.75+3.*max(-k*S,k-1.)),q=b*S,q+=a+sqrt(1.-dot(q,q))*normalize(V(-b.y,b.x,0)),C(a,q,3.5,2.5),C(q,a-b,2.5,2.)
void main(){float i,t,k,z,s,S=.5,d=S;for(V p,q,a,b;i++<5e1&&d>.001;t+=d=min(d,s=L+V(S-2.*p.x,-1,S))-S))p=normalize(V(P+P-R,R.y))*t,p.z-=5.,p.zy*=mat2(cos(vec4(1,12,34,1))),p.xz*=mat2(cos(sin(T)+vec4(0,11,33,0))),d=1.+p.y,C(z,V(z,z,1.2),7.5,6.),s=p.x<z?p.x=-p.x,z:H(z),s+=H(1.);O=vec4(V(exp(-i/(s>d?1e2:9.))),1);}

Посмотреть его можно на его официальной странице, а поэкспериментировать с кодом — в порте на Shadertoy.

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

Затем я использовал простую инверсную кинематику для анимации. Благодаря тому, что части ног имеют размер 1, формула упростилась и стала короче.

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

Что дальше

Я продолжу писать подобные демо, придушив свои художественные амбиции из-за выбранного ограничения в 512 символов.

Возможно, вы не понимаете, почему я одержим 512 символами, и многие люди задавали мне этот вопрос. На самом деле, на то есть множество причин:

В крошечном демо можно сделать упор только на один-два узких аспекта компьютерной графики, благодаря чему оно отлично подходит для обучения.
Это часть художественного исполнения: дело не только в техниках и визуальной составляющей, впечатляет ещё и волшебство, творящееся в коде. Мы живём в эпоху визуала, люди перекормлены самыми безумными VFX. Но видели ли они, как их можно реализовать всего в нескольких сотнях байт кода?
Такое ограничение помогает мне закончить работу: в творчестве всегда возникает вопрос, на чём следует остановиться. В данном случае это непреодолимая граница, дальше которой я ничего не могу сделать и вынужден переходить к чему-то новому.
Кроме того, это препятствует моим амбициям утянуть меня в какой-нибудь колоссальный проект, который я никогда не закончу, а может, даже не начну. Этот формат обладает кучей ограничений, и в этом его плюс.
Работа над таким крошечным фрагментом кода в течение дней/недель просто приносит мне удовольствие. Я ощущаю себя ремесленником, тратящим неразумное количество времени на совершенствование своего изделия ради красоты процесса.
Я пытаюсь создать портфолио, и для меня важно обеспечить его целостность. Если бы ограничение на размер было другим, я бы делал всё по-другому, поэтому сейчас уже ничего не могу поменять. Если бы у меня были ещё сотни символов в запасе, то в Red Alp могли появиться птицы, лучи, пробивающиеся сквозь облака и падающие на горы, и так далее.

Почему именно 512? Это размер toot в моём инстансе Mastodon, поэтому я могу уместить в него код; мне показалось, что он обеспечивает хороший баланс.

Трюки, которым я научился при создании маленьких GLSL-демо