|
|
Главная
Статьи
Ссылки
Скачать
Скриншоты
Юмор
Почитать
Tools
Проекты
Обо мне
Гостевая
Форум
|
В своей книге "Программирование компьютерной графики. Современный OpenGL" я рассматривал двухмерных полей расстояний со знаком (SDF, Signed Distance Field) для рендеринга текста. В трехмерном случае SDF также оказываются очень удобным и интересным инструментом. Они позволяют легко строить и визуализировать довольно сложную геометрию. В частности рендеринг ряда трехмерным фракталов основан именно на использовании SDF.
Итак 3D SDF - это фактически функция, задающая для каждой точки пространства расстояние со знаком до некоторой (иногда довольно сложной) поверхности. Если наша точка лежит снаружи, то значение функции больше нуля, если внутри - то меньше. На самой поверхности расстояние равно нулю (рис. 1).
Рис. 1. Функция расстояния со знаком
Проще всего задать SDF для плоскости и сферы:
float sphere ( in vec3 p )
{
return length ( p ) - 1.0;
}
float plane ( in vec3 p, vec4 n )
{
return dot ( p, n.xyz ) + n.w;
}
Также довольно легко можно задать SDF и для AABB (Axis Aligned Bounding Box):
float box ( in vec3 pos, in vec3 size )
{
vec3 pt = abs ( mv * pos ) - size;
return length ( max ( pt, 0.0 ) ) + min ( max ( pt.x, max ( pt.y, pt.z ) ), 0.0 );
}
Также можно слегка изменить эту функцию и получить SDF для AABB со скругленными краями и вершинами:
float boxRounded ( in vec3 pos, in vec3 size, float r )
{
vec3 q = abs ( mv * pos ) - size;
return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0) - r;
}
Можно довольно легко задать тор:
float torus ( in vec3 pos, in vec2 t )
{
vec3 pt = mv * pos;
vec2 q = vec2 ( length ( pt.xz) - t.x, pt.y );
return length ( q ) - t.y;
}
Ниже приводится задание SDF для цилиндра и конуса.
float cylinder ( in vec3 pos, in vec3 a, in vec3 b, in float r )
{
vec3 pt = abs ( mv * pos );
vec3 ba = b - a;
vec3 pa = pt - a;
float baba = dot(ba,ba);
float paba = dot(pa,ba);
float x = length(pa*baba-ba*paba) - r*baba;
float y = abs(paba-baba*0.5)-baba*0.5;
float x2 = x*x;
float y2 = y*y*baba;
float d = (max(x,y)<0.0)?-min(x2,y2):(((x>0.0)?x2:0.0)+((y>0.0)?y2:0.0));
return sign(d)*sqrt(abs(d))/baba;
}
float coneCapped ( in vec3 pos, in float h, in float r1, in float r2 )
{
vec3 pt = mv * pos;
vec2 q = vec2 ( length ( pt.xz ), pt.y );
vec2 k1 = vec2 ( r2, h );
vec2 k2 = vec2 ( r2 - r1, 2.0*h );
vec2 ca = vec2 ( q.x-min(q.x,(q.y<0.0)?r1:r2), abs(q.y)-h);
vec2 cb = q - k1 + k2*clamp( dot(k1-q,k2)/dot(k2, k2), 0.0, 1.0 );
float s = (cb.x<0.0 && ca.y<0.0) ? -1.0 : 1.0;
return s*sqrt( min(dot(ca, ca),dot(cb, cb)) );
}
Теперь давайте посмотрим на то, как мы можем визуализировать SDF. Оказывается самым простым способом этого будет трассировка лучей. Давайте посмотрим как это можно сделать.
Для того, чтобы строить изображения объекта методом трассировки лучей нам нужно уметь находить пересечения произвольного луча с этим объектом. Для SDF мы можем построить искомое пересечение итерационным способом, когда с каждым новым шагом мы будем все ближе подходить к точке пересечения (рис. 2).
Рис 2. Нахождение пересечения луча с SDF
Пусть у нас задана SDF f(p) и некоторый луч.
Тогда в качестве первого приближения мы выберем начало нашего луча p0.
Если 
, то мы считаем, что p0 и есть искомая точка пересечения.
Если нет, то заметим, что в шаре радиуса 
точки пересечения быть не может.
Это значит, что мы можем с качестве следующей точки пересечения взять пересечение луча с этой сферой.
Если наш луч задается уравнением 
, то следующее приближение p1 будет
задаваться формулой
Если 
, то p1 и есть искомая точка пересечения.
В противном случае мы знаем что внутри сферы с центром в p1 и радиусом 
пересечения быть не может.
Это значит, что мы можем сделать вдоль нашего луча шаг на :
Тем самым мы приходим к итерационному процессу.
Если для очередного приближения pi выполнено 
, то мы считаем, что
это и есть искомая точка пересечения.
Иначе мы строим следующее приближение pi+1 по формуле
Данный алгоритм можно реализовать в виде функции trace на языке GLSL:
#define MAX_ITERS 90
#define MAX_DIST 10.0
vec3 trace ( vec3 from, vec3 dir, out bool hit, out int steps )
{
vec3 p = from;
float totalDist = 0.0;
hit = false;
for ( steps = 0; steps < MAX_ITERS; steps++ )
{
float dist = dE ( p );
if ( dist < 0.001 )
{
hit = true;
break;
}
totalDist += dist;
if ( totalDist > MAX_DIST )
break;
p += dist * dir;
}
return p;
}
Эта функция возвращает точку пересечения с лучом, если она была найдена.
Через выходные параметры hit и steps
возвращается была ли найдена точка пересечения (hit=true) и за сколько шагов (steps).
Обратите внимание на ограничение по числу шагом - MAX_ITERS=90 и
максимальное расстояние от начала луча - MAX_DIST=10.
Они отвечают за остановку алгоритма в случае, когда пересечение не удается быстро найти.
Теперь давайте реализуем функцию main для нашего шейдера.
Она будет получать из вершинного шейдера во входной переменной p координаты точки на экране
и по ним и координатам начала eye будет строить луч для трассировки.
Далее построенный луч трассируется при помощи рассмотренной ранее функции trace и
по ее результатам определяется выходной цвет.
void main (void)
{
bool hit;
int steps;
vec3 org = eye;
vec3 dir = normalize ( p - eye );
vec3 light = vec3 ( 0.0, 0.0, 5.0 );
vec3 p = trace ( org, dir, hit, steps );
if ( hit )
color = vec4 ( 1.0 );
else
color = vec4 ( 0.0, 0.0, 0.0, 1.0 );
}
На следующем скриншоте приводится изображение для SDF AABB.
Рис 3. Построенное изображение для AABB.
Теперь давайте рассмотрим, как можно добавить освещение к нашему рендерингу SDF. Поскольку мы знаем точку пересечения, то мы можем в ней легко найти вектор к источнику света l. В качестве вектора v будет выступать направление, обратное направлению распространения луча.
Для расчета освещения нам не хватает только вектора нормали к поверхности в точке пересечения. Однако поскольку поверхность фактически является изоповерхностью (т.е. поверхностью уровня), то нормаль к ней будет нормированным градиентов SDF в точке пересечения. Для нахождения градиента мы воспользуемся разностным приближением к первой производной:
vec3 generateNormal ( vec3 z, float d )
{
float e = max(d * 0.5, EPS);
float dx1 = dE(z + vec3(e, 0, 0));
float dx2 = dE(z - vec3(e, 0, 0));
float dy1 = dE(z + vec3(0, e, 0));
float dy2 = dE(z - vec3(0, e, 0));
float dz1 = dE(z + vec3(0, 0, e));
float dz2 = dE(z - vec3(0, 0, e));
return normalize ( vec3 ( dx1 - dx2, dy1 - dy2, dz1 - dz2 ) );
}
Теперь мы моем добавить в наш код для трассировки полноценный расчет освещения.
void main (void)
{
bool hit;
int steps;
vec3 org = eye;
vec3 dir = normalize ( p - eye );
vec3 light = vec3 ( 0.0, 0.0, 5.0 );
vec3 p = trace ( org, dir, hit, steps );
if ( hit )
{
vec3 l = normalize ( light - p );
vec3 n = generateNormal ( p, 0.001 );
float nl = max ( 0.0, dot ( n, l ) );
color = vec4 ( nl );
}
else
color = vec4 ( 0.0, 0.0, 0.0, 1.0 );
}
Рис 4. Изображение SDF с рассчитанным освещением
В ряде SDF (например для тора) мы используем 2D и 3D метрику (расстояния)ю
На самом деле можно почить очень интересные фигуры, если заменить стандартную евклидову метрику 
на другую метрику.
В качестве возможных вариантов можно использовать метрики вида 
.
Ниже приводится пример для тора с нестандартной метрикой (код на GLSL и получившиеся изображения)
float lengthInf ( in vec2 p )
{
return max ( abs ( p.x ), abs ( p.y ) );
}
float length1 ( in vec2 p )
{
return abs ( p.x ) + abs ( p.y );
}
float length8 ( in vec2 p )
{
return pow ( pow ( abs(p.x), 8.0 ) + pow ( abs(p.y), 8.0 ), 1.0 / 8.0 );
}
float torus ( in vec3 pos, in vec2 t )
{
vec3 pt = mv * pos;
vec2 q = vec2 ( length ( pt.xz ) - t.x, pt.y );
return length8 ( q ) - t.y;
}
Еще одним "визуальным украшательством" обычно является наложение текстуры. Поскольку нас нет заранее заданных текстурных координат, то мы не можем воспользоваться традиционным подходом. Вместо этого мы применим так называемый triplanar texture mapping.
Сама идея triplanar texture mapping крайне проста. Давайте для каждой точки прочтем сразу три значения из текстуры, используя в качестве текстурных координат компоненты трехмерных координат точки пересечения - xy, xz и yz. После этого мы объединим их все вместе используя в качестве весов абсолютные значения компонент нормали в этой точке.
vec3 getTriplanarWeights ( in vec3 n )
{
vec3 w = abs ( n );
w *= w;
return w / ( w.x + w.y + w.z );
}
void main (void)
{
bool hit;
int steps;
vec3 org = eye;
vec3 dir = normalize ( p - eye );
vec3 light = vec3 ( 0.0, 0.0, 5.0 );
vec3 p = trace ( org, dir, hit, steps );
if ( hit )
{
vec3 l = normalize ( light - p );
vec3 n = generateNormal ( p, 0.001 );
float nl = max ( 0.0, dot ( n, l ) );
vec3 tx = getTriplanarWeights ( n );
vec3 q = mv * p;
vec4 cx = texture ( image, q.yz );
vec4 cy = texture ( image, q.zx );
vec4 cz = texture ( image, q.xy );
color = nl * (tx.x * cx + tx.y + cy + tx.z * cz);
}
else
color = vec4 ( 0.0, 0.0, 0.0, 1.0 );
}
Рис. 5. Текстурированное изображение тора.
Одним из больших плюсов SDF является то, что многие операции над множествами очень легко выражаются в терминах SDF. Одними из самых базовых операций над множествами являются операции объединения (union), пересечения (intersect) и вычитания (subtract). SDF для объединения двух (или большего числа) объектов, заданных при помощи SDF, определяется как минимум из соответствующих SDF:
float box ( in vec3 pos, in vec3 size )
{
vec3 pt = abs ( pos ) - size;
return length ( max ( pt, 0.0 ) ) + min ( max ( pt.x, max ( pt.y, pt.z ) ), 0.0 );
}
float sphere ( vec3 p, in vec3 center, in float radius )
{
return length ( p - center ) - radius;
}
float dE ( in vec3 p )
{
vec3 pt = mv * p;
float d1 = sphere ( pt, vec3 ( 0.0, 0.5, 0.0 ), 0.55 );
float d2 = box ( pt, vec3 ( 0.6, 0.2, 0.7 ) );
return max ( d1, d2 ); // intersection
//return min ( -d1, d2 ); // union
//return max ( -d1, d2 ); // subtraction
}
Рис 6. Вычитание двух объектов
С пересечением и вычитанием дело обстоит несколько сложнее - соответствующие SDF можно очень легко определить во внешней области (где значения больше 0). Во внутренней области эти формулы дают неверные результаты. Поэтому если вам нужно просто отрендерить результат этих операций снаружи, то эти формулы отлично работают. Однако если вы хотите использовать их так, что ван нужно и значение во внутренней области, т о вам нужно в этой области данные формулы подкорректировать.
Для пересечения такой (частично верной) функцией будет просто минимум из двух значений. Для вычитания - .
Рис 7. Результат пересечения двух объектов, заданный SDF
При помощи этих операций можно построить и довольно сложные объекты, как показано на следующем рисунке.
Ниже приводится код, задающий SDF для этого объекта.
float box ( in vec3 pos, in vec3 size )
{
vec3 pt = abs ( pos ) - size;
return length ( max ( pt, 0.0 ) ) + min ( max ( pt.x, max ( pt.y, pt.z ) ), 0.0 );
}
float sphere ( in vec3 pos, in vec3 center, in float radius )
{
return length ( pos - center ) - radius;
}
float cylinder ( in vec3 pos, in vec3 a, in vec3 b, in float r )
{
vec3 pt = abs ( pos );
vec3 ba = b - a;
vec3 pa = pt - a;
float baba = dot ( ba, ba );
float paba = dot ( pa, ba );
float x = length ( pa*baba - ba*paba ) - r*baba;
float y = abs ( paba - baba*0.5 ) - baba*0.5;
float x2 = x*x;
float y2 = y*y*baba;
float d = (max(x,y)<0.0)?-min(x2,y2):(((x>0.0)?x2:0.0)+((y>0.0)?y2:0.0));
return sign(d)*sqrt(abs(d))/baba;
}
float dE ( in vec3 pos )
{
vec3 pt = mv * pos;
float v1 = cylinder ( pt, vec3 ( -1, 0, 0 ), vec3 ( 1, 0, 0 ), 0.11 );
float v2 = cylinder ( pt, vec3 ( 0, -1, 0 ), vec3 ( 0, 1, 0 ), 0.11 );
float v3 = cylinder ( pt, vec3 ( 0, 0, -1 ), vec3 ( 0, 0, 1 ), 0.11 );
float vb = box ( pt, vec3 ( 0.28 ) );
float vs = sphere ( pt, vec3 ( 0.0 ), 0.35 );
float v = max ( vb, vs );
return max ( max ( max ( v, -v1 ), -v2 ), -v3 );
}
У операций объединения, пересечения и вычитания есть также их "сглаженные" варианты - в них возникающие ребра будут сглаженными.
Рис 8. Сглаженные варианты операций над объектами.
Ниже приводится код для сглаженных операций над SDF-объектами.
float smoothUnion ( float d1, float d2, float k )
{
float h = clamp( 0.5 + 0.5*(d2-d1)/k, 0.0, 1.0 );
return mix( d2, d1, h ) - k*h*(1.0-h);
}
float smoothSubtraction ( float d1, float d2, float k )
{
float h = clamp( 0.5 - 0.5*(d2+d1)/k, 0.0, 1.0 );
return mix( d2, -d1, h ) + k*h*(1.0-h);
}
float smoothIntersection ( float d1, float d2, float k )
{
float h = clamp( 0.5 - 0.5*(d2-d1)/k, 0.0, 1.0 );
return mix( d2, d1, h ) + k*h*(1.0-h);
}
Давайте сейчас рассмотрим как преобразование поверхности влияет на SDF. Пусть у нас есть SDF f(p) и мы хотим получить SDF соответствующую переносу f на заданный вектор t. Как несложно убедиться поверхность, полученная путем переноса, будет задаваться SDF f(p+t).
Чуть сложнее обстоит дело в поворотом. Пусть у нас есть преобразование поворота, задаваемое матрицей 3x3 R. Тогда соответствующая функция будет равна f(R-1 p).
Однородное масштабирование на коэффициент s дает нам следующую SDF: f(p/s)*s.
Теперь представьте, что у нас есть некоторый небольшой объект, заданный при помощи SDF f(p) и мы хотим построить бесконечную матрицу из этого объекта с шагом по каждой из осей hx, hy и hz. Ниже приводится соответствующий код.
float box ( in vec3 pos, in vec3 size )
{
vec3 pt = abs ( pos ) - size;
return length ( max ( pt, 0.0 ) ) + min ( max ( pt.x, max ( pt.y, pt.z ) ), 0.0 );
}
float dE ( in vec3 p )
{
const vec3 c = vec3 ( 0.5, 0.5, 0.2 );
vec3 pt = mv * p;
vec3 q = mod ( pt + 0.5 * c, c ) - 0.5 * c;
return box ( q, vec3 ( 0.3, 0.3, 0.1 ) );
}
При помощи SDF мы можем легко "округлить" примитив - для этого достаточно просто из SDF примитива вычесть радиус скругления.
Рис 9. Исходный и скругленный примитивы
Кроме традиционный операций объединения, пересечения и вычитания при помощи SDF можно очень легко реализовать следующие "сглаженные" версии этих операций.
Рис 10. Сглаженные операции
float box ( in vec3 pos, in vec3 size )
{
vec3 pt = abs ( pos ) - size;
return length ( max ( pt, 0.0 ) ) + min ( max ( pt.x, max ( pt.y, pt.z ) ), 0.0 );
}
float sphere ( vec3 p, in vec3 center, in float radius )
{
return length ( p - center ) - radius;
}
float smoothUnion ( float d1, float d2, float k )
{
float h = clamp( 0.5 + 0.5*(d2-d1)/k, 0.0, 1.0 );
return mix( d2, d1, h ) - k*h*(1.0-h);
}
float smoothSubtraction ( float d1, float d2, float k )
{
float h = clamp( 0.5 - 0.5*(d2+d1)/k, 0.0, 1.0 );
return mix( d2, -d1, h ) + k*h*(1.0-h);
}
float smoothIntersection ( float d1, float d2, float k )
{
float h = clamp( 0.5 - 0.5*(d2-d1)/k, 0.0, 1.0 );
return mix( d2, d1, h ) + k*h*(1.0-h);
}
float dE ( in vec3 p )
{
vec3 pt = mv * p;
float d1 = sphere ( pt, vec3 ( 0.0, 0.5, 0.0 ), 0.55 );
float d2 = box ( pt, vec3 ( 0.6, 0.2, 0.7 ) );
//return smoothUnion ( d1, d2, 0.3 ); // intersection
return smoothSubtraction ( d1, d2, 0.2 ); // union
//return smoothIntersection ( d1, d2, 0.2 ); // subtraction
}
Очень красивый эффект получается, когда мы добавляем какой-либо шум к SDF. То, что мы при этом получаем, строго говоря уже не является SDF. Это значит, что если мы будем трассировать эту поверхность рассмотренным ранее способом, то мы можем легко проткнуть получившуюся поверхность насквозь. К счастью почти всегда можно решить эту проблему очень просто - нужно просто делать шаг вдоль луча не на абсолютное значение функции, а на 70-80% от нее. Ниже приводится пример поверхности, к которой добавлена турбулентность.
Рис 11. Поверхность с турбулентностью.
float smin( float a, float b, float k )
{
float h = max(k-abs(a-b),0.0);
return min(a, b) - h*h*0.25/k;
}
float smax( float a, float b, float k )
{
float h = max(k-abs(a-b),0.0);
return max(a, b) + h*h*0.25/k;
}
float sph( vec3 i, vec3 f, vec3 c )
{
// random radius at grid vertex i+c (please replace this hash by
// something better if you plan to use this for a real application)
vec3 p = 17.0*fract( (i+c)*0.3183099+vec3(0.11,0.17,0.13) );
float w = fract( p.x*p.y*p.z*(p.x+p.y+p.z) );
float r = 0.7*w*w;
// distance to sphere at grid vertex i+c
return length(f-c) - r;
}
float sdBase( vec3 p )
{
ivec3 i = ivec3(floor(p));
vec3 f = fract(p);
// distance to the 8 corners spheres
return min(min(min(sph(i,f,ivec3(0,0,0)),
sph(i,f,ivec3(0,0,1))),
min(sph(i,f,ivec3(0,1,0)),
sph(i,f,ivec3(0,1,1)))),
min(min(sph(i,f,ivec3(1,0,0)),
sph(i,f,ivec3(1,0,1))),
min(sph(i,f,ivec3(1,1,0)),
sph(i,f,ivec3(1,1,1)))));
}
float plane ( vec3 p, vec4 n )
{
return dot ( p, n.xyz ) + n.w;
}
float dE ( vec3 pos )
{
vec3 p = mv * pos;
float s = 1.0;
float d = plane ( p, vec4 ( 0, 1, 0, 1 ) );
for ( int i = 0; i < 8; i++ )
{
// evaluate new octave
float n = s*sdBase(p);
// add
n = smax(n, d-0.1*s,0.3*s);
d = smin(n, d ,0.3*s);
// prepare next octave
p = mat3( 0.00, 1.60, 1.20,
-1.60, 0.72,-0.96,
-1.20,-0.96, 1.28 )*p;
s = 0.5*s;
}
return d;
}
Можно, используя шум (точнее функцию, построенную на основе него - fBm) построить и такой объект:
Ниже приводится соответствующий код:
float fbm( vec3 p, float d )
{
float s = 0.5;
for( int i = 0; i < 7; i++ )
{
// evaluate new octave
float n = s*base ( p );
// subtract
d = smax( d, -n, 0.2*s );
// prepare next octave
p = m * p;
s = 0.5 * s;
}
return d;
}
float dE ( in vec3 pos )
{
vec3 p = mv * pos;
float d = box ( p, vec3 ( 0.5 ) );
float dt = fbm ( p + 0.5, d );
return dt;
}
При помощи SDF можно довольно легко строить различные фрактальные структуры. Одной из таких структур является губка Менгера (Menger Spone). Процедура ее построения довольно проста. Мы начинаем с единичного куба, который делится плоскостями, параллельными координатным, на 27 одинаковых подкубов.
После этого мы удаляем центральный подкуб и все подкубы, имеющие общую с ним грань.
Поучаем множество, состоящее из 20 кубов.
Далее к каждому из этих кубов мы применяем ту же самую процедуру и т.д.
Этот алгоритм довольно легко выразить в терминах SDF и реализовать на GLSL
(параметр NUM_ITER отвечает за число итераций данного алгоритма).
float box ( in vec2 pos, in vec2 size )
{
vec2 pt = abs ( pos ) - size;
return length ( max ( pt, 0.0 ) ) + min ( max ( pt.x, pt.y ), 0.0 );
}
float sdCross ( in vec3 p )
{
float da = box ( p.xy, vec2 ( 1.0 ) );
float db = box ( p.yz, vec2 ( 1.0 ) );
float dc = box ( p.zx, vec2 ( 1.0 ) );
return min ( da, min ( db, dc ) );
}
float dE ( in vec3 pos )
{
vec3 p = mv * pos;
float d = box ( p, vec3(1.0) );
vec3 res = vec3 ( d, 1.0, 0.0 );
float s = 1.0;
for( int m = 0; m < 5; m++ )
{
vec3 a = mod ( p*s, 2.0 ) - 1.0;
vec3 r = abs ( 1.0 - 3.0*abs(a) );
s *= 3.0;
float da = max ( r.x, r.y );
float db = max ( r.y, r.z );
float dc = max ( r.z, r.x );
float c = (min(da,min(db,dc)) - 1.0) / s;
if ( c > d )
{
d = c;
res = vec3 ( d, 0.2*da*db*dc, (1.0+float(m))/4.0 );
}
}
return res.x;
}
Рис 12. Губка Менгера
Есть очень простой и красивый метод, добавляющий реализма изображению - расчет фонового затенения (ambient occlusion). И для объектов, заданных при помощи SDF это можно очень легко реализовать. Для этого мы будем двигаться от текущей точки наружу вдоль вектора нормали и при помощи SDF проверять есть и поблизости закрывающая геометрия. Если мы отошли от исходной точки p на расстояние d, то мы смотрим насколько сильно отличается d от абсолютного значения SDF в этой точке. Чем больше разница, тем больше затеняющей геометрии поблизости. Обычно делается несколько шагов и суммируется вклад каждого из них в затенение (рис 13).
Рис 13. Расчет фонового затенения
float ambientOcclusion ( in vec3 pos, in vec3 normal )
{
float occ = 0.0;
float sca = 1.0;
for ( int i = ZERO; i < 5; i++ )
{
float h = 0.01 + 0.12*float(i)/4.0;
float d = dE ( pos + h*normal );
occ += (h-d)*sca;
sca *= 0.95;
if ( occ > 0.35 )
break;
}
return clamp ( 1.0 - 3.0*occ, 0.0, 1.0 ) * (0.5+0.5*normal.y);
}
Рис 14. Губка Менгера с фоновым затенением
Весь исходный код к этой статье можно найти в каталоге python-code в репозитории.
Copyright © Alexey V. Boreskov 2003-2022