Рендеринг воды

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

Наиболее важными являются моделирование цвета воды, отражения поверхностью воды (а также в ряде случаев и преломления), анимация воды.

Цвет воды (считая ее достаточно глубокой, чтобы дна не было видно) зависит от угла между направлением на наблюдателя v и нормалью n в данной точке.

Рис 1.

Обычно цвет меняется от слегка зеленоватого оттенка (при углах близким к 90 градусам) до темно-синих (для углов близких к нулю).

Наиболее простым способом получения цвета воды, соответствующего заданному углу (удобнее вместо самого угла использовать его косинус, т.е. скалярное произведение (n,v)) является линейная интерполяция между двумя крайними значениями цвета:

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

При моделировании этих эффектов важную роль играет коэффициент Френеля, задающий какая доля падающей световой энергии отражается для заданного угла падения.

Существует несколько приближенных способов вычисления коэффициентов Френеля, но далее мы будем использовать следующую формулу:

Наиболее сложным элементом является анимация волн. В самом простейшем случае для этого можно использовать функцию turbulence, вводя в нее зависимость от трех переменных - двух пространственных (x,y) и времени (t).

Одной из наиболее правдоподобных моделей анимации воды является так называемая статистическая модель. Одним из первых ее применений был рендеринг воды в фильме "Титаник".

Ниже приводится краткое описание этой модели, полное ее описание можно найти в статьях Jerry Tessendorf. Simulating Ocean Water и "Lasse Jensen. Deep-Water Animation and Rendering".

Статистическая модель основана на анализе очень большого числа реальных наблюдений. При этом высота волны h(x,y,t) рассматривалась как случайная величина, закон распределения которой требуется установить.

При этом сама волна представлялась как сумма большого числа гармоник (т.е. фактически осуществлялось разложение в ряд Фурье по пространственным переменным x и y).

Через k обозначен двухмерный вектор следующего вида:

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

Существует несколько моделей для спектра Ph(k), одной из наиболее распространенных является спектр Филлипса:

Через обозначена наибольшая высота волны, возникающей от постоянного ветра со скоростью V, а через - направление ветра. Величина A является нормирующей константой.

Исходя из спектра Филлипса можно легко построить требуемые функции в виде:

Используемые комплекснозначные величины h0(k) можно построить по следующей формуле:

Через обозначены распределенные по Гауссу (со нулевым средним значением и единичным среднеквадратичным отклонением) случайные величины.

Связь между частотой (по времени) и вектором k обычно задается при помощи следующей формулы (через g обозначено ускорение свободного падения):

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

Основная задача заключается в быстром нахождении суммы гармоник с заданными коэффициентами, что на легко реализуется при помощи быстрого преобразование Фурье (FFT, Fast Fourier Transform).

В статье "Jason Mitchell'а Real-Time Synthesis and Rendering of Ocean Water." рассматривается реализация данного метода целиком на GPU (за исключением первоначального набора используемых далее коэффициентов). При этом быстрое преобразование Фурье также вычисляется на GPU.

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

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

Шумовая функция при этом может быть задана при помощи 3D-текстуры и реализация данной функции на GPU является простой и достаточно эффективной.

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

uniform samplerCube envMap;
uniform float       time;
uniform float       amplitude;

varying vec3  n;
varying vec3  e;
varying	vec3  p;

uniform sampler3D noiseMap;

vec3 turbulence ( const in vec3 p, const in float freqScale )
{
    float sum = 0.0;
    vec3  t   = vec3 ( 0.0 );
    float f   = 1.0;

    for ( int i = 0; i <= 3; i++ )
    {
        t   += abs ( 2.0 * texture3D ( noiseMap, f * p.xyz ).rgb - vec3 ( 1.0 ) ) / f;
        sum += 1.0 / f;
        f   *= freqScale;
    }
                                 // remap from [0,sum] to [-1,1]
    return 2.0 * t / sum - vec3 ( 1.0 );
}

void main (void)
{
    const vec4 c2 = 2.2 * vec4 ( 0.03, 0.15, 0.125, 1.0 );
    const vec4 c1 = 1.6 * vec4 ( 0.03, 0.2,  0.07, 1.0 );

                                 // Compute reflection vector
    vec3 arg = 0.25*p + vec3 ( 0.04101, -0.0612149, 0.071109 ) * time * 0.7;
    vec3 ns  = vec3 ( turbulence ( arg, 2.17 ).xy, 1.0 );
    vec3 en  = normalize ( e );
    vec3 nn  = normalize ( n + amplitude * ns );
    vec3 r   = reflect   ( en, nn );

                                 // Do a lookup into the environment map.
    vec4 envColor = textureCube ( envMap, r );

                                 // calc fresnels term.  This allows a view dependant
                                 // blend of reflection/refraction
    float fresnel = clamp ( abs ( dot ( en, nn ) ), 0.1, 0.9 );

    gl_FragColor = vec4 ( mix ( c1, c2, abs ( dot ( en, nn ) ) )*fresnel +  (1.0-fresnel) * envColor );
}

Ниже приводится скриншот воды, полученной таким образом.

Рис 2. Изображением воды, полученное при помощи EMBM с использованием функции turbulence для задания карты нормалей.

Одним из существенных минусов такого подхода является слишком "плоский" характер воды - вся анимация реализуется только на уровне bumpmap'а.

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

Ниже приводится фрагментный шейдер, осуществляющий построение массива вершин и нормалей.

#define   SIZE      16.0
#define   TEX_SIZE  128.0

uniform sampler3D noiseMap;
uniform float     time;
uniform float     amplitude;

vec3 turbulence ( const in vec3 p, const in float freqScale )
{
    float sum = 0.0;
    vec3  t   = vec3 ( 0.0 );
    float f   = 1.0;

    for ( int i = 0; i <= 3; i++ )
    {
        t   += abs ( 2.0 * texture3D ( noiseMap, f * p.xyz ).rgb - vec3 ( 1.0 ) ) / f;
        sum += 1.0 / f;
        f   *= freqScale;
    }
                                        // remap from [0,sum] to [0,1]
    return t / sum;
}

void	main ()
{                                       // deltas
    const vec3 dx  = vec3 ( 0.01, 0.0,  0.0 );
    const vec3 dy  = vec3 ( 0.0,  0.01, 0.0 );
    const vec3 vel = vec3 ( 0.02, 0.01, 0.0 );
	
                                        // compute height and delta height (in dx and dy)
    vec3 tex    = vec3       ( gl_TexCoord [0].xy, 0.03*time ) + vel * time;
    vec3 turb   = turbulence ( tex,      2.0 );
    vec3 turbX  = turbulence ( tex + dx, 2.0 );
    vec3 turbY  = turbulence ( tex + dy, 2.0 );
    vec3 n      = vec3       ( turbX.x - turb.x, turbY.x - turb.x, 0.8 );

    gl_FragData [0] = vec4 ( gl_TexCoord [0].xy * 2.0 * SIZE - vec2 ( SIZE ), -1.0 + amplitude * turb.x, 1.0 );
    gl_FragData [1] = vec4 ( normalize ( n ), 1.0 );
}

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

Рис 3. Изображение воды, полученное при помощи анимации вершин функцией turbulence.

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

Рис 4. Изображения воды, получаемые при Использовании функции turbulence как для анимации вершин, так и для задание карты нормалей.

Однако к сожалению, анимация, получаемая при помощи функции turbulence все-таки достаточно сильно отличается как от настоящей воды, так и от результатов статистической модели.

Интересным компромиссным вариантом является подход, описанный Юрием Крячко в книге GPU Gems 2 (гл.18 ). Описываемый им подход основывается на использовании набора карт высот для заранее рассчитанной анимации (на следующем рисунке приводится одна из таких карт высот).

Рис. 5. Одна из заранее рассчитанных карт высот для анимации воды.

При этом эти текстуры образуют периодическую (как в пространстве, так и по времени) анимацию достаточно высокого качества (и реализма, ведь их построение происходит заранее и поэтому может использовать достаточно сложные методы).

Самый простой способ использование подобных текстур заключается в том, что для каждого момента времени t берутся две соседние карты Hi и Hi+1 и между ними осуществляется линейная интерполяция.

Так если у нас задан период анимации T, то для момента времени t в качестве карт высот следует взять карты со следующими номерами (через N обозначено общее количество карт высоты):

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

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

Обозначим функцию высоты, получаемую описанным выше способом (т.е. при помощи интерполяции между двумя последовательными картами высот) через h(x,y,t).

Тогда довольно легко построить следующую функцию.:

Константы A и B обеспечивают "сдвиг по частоте и фазе", за счет чего период получаемой функции возрастает на порядки (хотя она по-прежнему является периодической).

Полученная таким образом функция H(x,y,t) не только практически непериодична (т.е. обладает очень большим периодом), но и содержит в себе сразу несколько уровней детализации (как и функция turbulence).

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

Достаточно простым способом вычисления функции H(x,y,t) является использование сразу нескольких пар соседних карт высот из анимационной последовательности, каждая такая пара соответствует одному слагаемому в итоговой сумме.

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

#define  SIZE      16.0
#define  TEX_SIZE  128.0

uniform sampler2D h1Map;                        // height maps
uniform sampler2D h2Map;
uniform sampler2D h3Map;
uniform sampler2D h1Map2;                       // next frame height maps
uniform sampler2D h2Map2;
uniform sampler2D h3Map2;
uniform float     time;
uniform float     amplitude;

vec3 clr2vec3 ( const vec3 clr )
{
    return 2.0 * clr - vec3 ( 1.0 );
}

void main ()
{                                               // deltas
    const vec3 dx  = vec3 ( 0.01, 0.0,  0.0 );
    const vec3 dy  = vec3 ( 0.0,  0.01, 0.0 );
    const vec2 vel = vec2 ( 0.02, 0.01 )*0.2;
    const vec2 b1  = vec2 ( 0.1074, 0.7083 );  // texCooord biases
    const vec2 b2  = vec2 ( 0.3571, 0.1144 );
//  const vec2 b3  = vec2 ( 0.5037, 0.1974 );

                                                // compute height and delta height (in dx and dy)
    vec2  tex    = vec2  ( gl_TexCoord [0].xy + vel * time );
    vec2  tex2   = vec2  ( gl_TexCoord [0].xy + vel * time );
    vec2  tex3   = vec2  ( gl_TexCoord [0].xy + vel * time );
    float t      = fract ( 2.0 * time );
    float t2     = fract ( 4.0 * time );
    float t3     = fract ( 8.0 * time );
    vec4  h1     = mix   ( texture2D ( h1Map, tex + b1 ),       texture2D ( h1Map2, tex + b1 ), t );
    vec4  h2     = mix   ( texture2D ( h2Map, 2.0 * tex2 + b2 ), texture2D ( h2Map2, 2.0 * tex2 + b2 ), t2 );
//  vec4  h3     = mix   ( texture2D ( h3Map, 4.0 * tex3 + b3 ), texture2D ( h1Map2, 4.0 * tex3 + b3 ), t3 );
    vec4  h      = 0.5*h1 + 0.25 * h2;
    vec3  n      = clr2vec3 ( h1.xyz ) + clr2vec3 ( h2.xyz );// + clr2vec3 ( h3 );

    gl_FragData [0] = vec4 ( gl_TexCoord [0].xy * 2.0 * SIZE - vec2 ( SIZE ), -2.0 + amplitude * h.w, 1.0 );
    gl_FragData [1] = vec4 ( normalize ( n ), 1.0 );
}

Ниже приводится фрагмент кода на С++, используемый для нахождения соответствующих пар индексов текстур.

float time   = 0.001f * glutGet ( GLUT_ELAPSED_TIME );  // time in seconds
int   index  = ((int) (SCALE1*curTime)) & 0x0F;
int   index2 = ((int) (SCALE2*curTime)) & 0x0F;
int   index3 = ((int) (SCALE3*curTime)) & 0x0F;

Точно также функция H(x,y,t) может использоваться и для построения анимированной карты нормалей, соответствующий фрагментный шейдер приводится ниже.

#define  SIZE       16.0
#define  TEX_SIZE   128.0
#define  SCALE1     2.073
#define  SCALE2     4.235
#define  SCALE3     8.0973

uniform sampler2D h1Map;                           // height maps
uniform sampler2D h2Map;
uniform sampler2D h3Map;
uniform sampler2D h1Map2;                          // next frame height maps
uniform sampler2D h2Map2;
uniform sampler2D h3Map2;
uniform float     time;
uniform float     amplitude;

vec3	clr2vec3 ( const vec3 clr )
{
    return 2.0 * clr - vec3 ( 1.0 );
}

void	main ()
{                                                  // deltas
    const vec3 dx  = vec3 ( 0.01, 0.0,  0.0 );
    const vec3 dy  = vec3 ( 0.0,  0.01, 0.0 );
    const vec2 vel = vec2 ( 0.02, 0.01 )*0.2;
    const vec2 b1   = vec2 ( 0.1074, 0.7083 );     // texCooord biases
    const vec2 b2   = vec2 ( 0.3571, 0.1144 );

                                                   // compute height and delta height (in dx and dy)
    vec2  tex    = vec2  ( gl_TexCoord [0].xy + vel * time );
    vec2  tex2   = vec2  ( gl_TexCoord [0].xy + vel * time );
    vec2  tex3   = vec2  ( gl_TexCoord [0].xy + vel * time );
    float t      = fract ( SCALE1 * time );
    float t2     = fract ( SCALE2 * time );
    float t3     = fract ( SCALE3 * time );
    vec4  h1     = mix   ( texture2D ( h1Map, SCALE1 * tex + b1 ),  texture2D ( h1Map2, SCALE1 * tex  + b1 ), t );
    vec4  h2     = mix   ( texture2D ( h2Map, SCALE2 * tex2 + b2 ), texture2D ( h2Map2, SCALE2 * tex2 + b2 ), t2 );
    vec4  h      = 0.5*h1 + 0.25 * h2;
    vec3  n      = clr2vec3 ( h1.xyz ) + clr2vec3 ( h2.xyz );

    gl_FragData [0] = vec4 ( gl_TexCoord [0].xy * 2.0 * SIZE - vec2 ( SIZE ), -2.0 + amplitude * h.w, 1.0 );
    gl_FragData [1] = vec4 ( normalize ( n ), 1.0 );
}

Ниже приводятся скриншоты воды, получаемые при таком подходе:

Рис 6. Изображения, получающиеся при использовании набора карт высоты как для анимации вершин, так и для задания анимированной карты нормалей.

Рис 7. HDR-версия.

По этой ссылке можно скачать весь исходный код к этой статье. Также доступны для скачивания откомпилированные версии для M$ Windows, Linux.

Используются технологии uCoz