Бесплатная реклама

Автор: Дмитрий Трифонов «DikobrAz»

Опубликовано: 11.04.2010

Изменено: 11.04.2010

Постоянная ссылка

Комментарии [24]

Визуализация водной поверхности. Быстрое преобразование Фурье на GPU





Страницы: 1 2 3 4

3. Быстрое преобразование Фурье на GPU

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

    3.1. Дискретное преобразование Фурье

Рассмотрим дискретное одномерное комплексное преобразование Фурье, далее ДПФ:



Формула 3. Прямое и обратное дискретное комплексное преобразование Фурье.
Весовые коэффициенты.

Если вычислять преобразование Фурье напрямую, потребуется O(N2) операций. При реализации статистической модели синтеза поверхности океана будут использоваться комплексные сетки размера N=128x128(N2 = 268 435 456) и больше, - такая сложность не приемлема для вычислений в реальном времени.

    3.2 Быстрое преобразование Фурье

Алгоритм быстрого преобразование Фурье(Cooley and Tukey "Decimation In Time") – далее БПФ, основывается на стратегии “разделяй и властвуй”. Он не самый эффективный по количеству операций, но легко распараллеливается и не требует сложных вспомогательных вычислений, реализация которых на GPU может оказаться затруднительной или даже невозможной.

Рассмотрим первый шаг алгоритма. Просуммируем отдельно слагаемые с четными и нечетными степенями весовых коэффициентов, полагая что N делится на два:


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



Рассмотрим ещё несколько вспомогательных свойств весовых коэффициентов:




Используя последние равенствa несложно свести задачу ДПФ от N коэф-фициентов к ДПФ от N/2 коэффициентов:


Формула 4. Разбиение ДПФ от N элементов на два.

Таким образом, вычисляя F(k) только для k=0..N/2 мы можем заполнить весь массив коэффициентов, уменьшив тем самым сложность вдвое.

Оба получившихся ДПФ можно разбить ещё на два, если N кратно четырем. Если N является степенью двойки, то можно свести искомое преобразование Фурье к ДПФ от одного элемента. Сложность БПФ в таком случае будет Θ(N log N).

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

    3.3 Реализация с использованием шейдеров

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

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

На нулевом шаге алгоритма в исходном массиве N коэффициентов ДПФ, фактически N преобразований Фурье размера 1. На первом необходимо посчитать N/2 преобразований Фурье размера 2, на втором N/4 размера 4 и так далее. Итого, на i’м шаге алгоритм работает N/2i разбиениями исходной последовательности длины 2i.

Для слияния двух ДПФ в одно, более высокого порядка, можно воспользо-ваться формулой 4, но необходимо ещё учитывать, что на каждом следующем шаге алгоритма будет меняться расположение коэффициентов преобразования. Классический метод индексации для БПФ называется bit-reversing. Суть его состоит в предварительной перестановке исходных элементов массива так, чтобы на каждой итерации коэффициенты для любого ДПФ были расположены последовательно.

Для реализации на GPU лучше избавиться от предварительной перестановки, используемой в этом методе. Есть способ определения индексов для итеративного вычисления БПФ без каких либо модификаций исходного массива.


Рисунок 12. Разбиение ДПФ. Сначала поделили элементы на две группы: с четными индексами и с нечетными.
Затем каждую группу разбили ещё на две, уже смотря на четность или нечетность индекса внутри группы, и т. д.

Выше на рисунке изображено последовательное разбиение исходного массива: синяя группа разбивается на синюю и желтую, полученная синяя на синюю и красную, а желтая на желтую и зеленую и т. д. Если же идти снизу вверх, объединяя в группы ячейки простым суммированием, обратно тому как производилось их разбиение, то на любом уровне в каждой ячейке будет хранится сумма тех элементов, соответствующие ячейки которых имеют такой же цвет. Таким образом, сумма значений 2 и 3 ячейки на втором уровне – есть сумма всех элементов массива. На самом верхнем уровне, в каждой ячейке хранится сумма всех элементов массива.

Теперь необходимо, чтобы объединение ячеек производило объединение ДПФ. Это возможно сделать, если при суммировании на каждом этапе домножать слагаемые на необходимые весовые коэффициенты. Формула 5, как и формула 4, объединяет два преобразования Фурье, полученных на предыдущем этапе, но уже с учетом индексации в массиве.


Формула 5. Определение индексов элементов для БПФ, i – номер итерации, N – размер ДПФ. Индексы берем по модулю N.

Строгое обоснование приведенной выше формулы и более подробное описание метода описано в статье "The FFT on a GPU" Kenneth Moreland & Edward Angel. Ниже приведен формальный пример вычисления второго элемента ДПФ, с использованием этой формулы, для N=8:


Алгоритм для вычисления двумерного БПФ средствами графического API выглядит следующим образом:

renderTarget->SetDrawBuffer(0, pingPongTexture[0]);
renderTarget->SetDrawBuffer(1, pingPongTexture[1]);
BindRenderTarget(renderTarget);
   
int numIterations = (int)(log( (double)N ) / log(2.0) + 0.5); // log2(N)
fftProgram->Bind();  // Шейдер, выполняющий формулу 5
sizeUniform->Set(N); // Размер преобразования

// БПФ по строкам
int curPing = 0;
directionUniform->Set(1.0f, 0.0f); // Направление БПФ
for(int i = 1; i <= numIterations; ++i)
{
    inputUniform->Set(0, pingPongTexture[curPing]); // Входной массив
    numIterationUniform->Set(i); // Номер итерации
    curPing = !curPing;
    renderTarget->SetDrawBuffer(curPing); // Меняем выходную текстуру
    gridQuad->Draw(QUADS);
}

// БПФ по столбцам
directionUniform->Set(0.0f, 1.0f); // Установить направление по столбцам
for(int i = 1; i <= numIterations; ++i)
{
    inputUniform->Set(0, pingPongTexture[curPing]);
    numIterationUniform ->Set(i);
    curPing = !curPing;
    renderTarget->SetDrawBuffer(curPing);
    gridQuad->Draw(QUADS);
}

Шейдеры, выполняющие преобразование Фурье(выполняется два двумерных преобразования Фурье, записанных в одной текстуре формата RGBA32F или RGBA16F - каждый элемент занимает два числа, т.к. элементы комплексные):

/=============== Vertex Shader ===============/
#extension GL_EXT_gpu_shader4 : enable // используем, если есть

uniform vec2  viewport;
uniform vec2  direction; // направление: столбцы или строки
uniform float N;         // размер преобразования(одномерного)

#if GL_EXT_gpu_shader4
noperspective varying vec2  index;
noperspective varying float linearIndex;
#else
varying vec2    index;       // двумерный индекс
varying float   linearIndex; // одномерный индекс
#endif

void main()
{
    vec2 gridVertex = N * gl_Vertex.xy;

    linearIndex = gridVertex.x;
    index       = gridVertex.x * direction.xy + gridVertex.y * direction.yx;
    gl_Position = vec4( 2.0 * ( index.xy / N - vec2(0.5) ), 0.0, 1.0 );
}

/=============== Fragment Shader ===============/
#extension GL_EXT_gpu_shader4 : enable

uniform sampler2D fftInput;

uniform vec2  direction;
uniform float N;
uniform float numPartitions; // колличество выполняемых ДПФ. Равно N / 2^i, i - номер прохода

#if GL_EXT_gpu_shader4
noperspective varying vec2  index;
noperspective varying float linearIndex;
#else
varying vec2    index;
varying float   linearIndex;
#endif

#define PI     3.141592653
#define PI_2   (2.0 * PI)

// Просто умножение комплексных чисел
vec2 complex_mult(vec2 a, vec2 b)
{
    return vec2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x);
}

// Рассчет весового коэффициента
vec2 W(float partitionSize, float X)
{
    float angle = PI_2 * X / partitionSize;
    return vec2( cos(angle), sin(angle) );
}

void main()
{
    float k = floor(linearIndex / numPartitions);

    // fftIndex = ( iIndex + k * n ) % N
    vec2 fftIndex = index + direction * k * numPartitions;
    fftIndex     -= floor(fftIndex / N) * N;

    // Выбираем коэффициенты из текстуры
#if GL_EXT_gpu_shader4
    vec4 He = texelFetch(fftInput, ivec2(fftIndex), 0);  // Выборка без фильтрации
    vec4 Ho = texelFetch(fftInput, ivec2(fftIndex + direction * numPartitions), 0);
#else
    vec4 He = texture2DLod(fftInput, fftIndex / N, 0.0); // Выборка без фильтрации
    vec4 Ho = texture2DLod(fftInput, (fftIndex + direction * numPartitions) / N, 0.0);
#endif

    // Преобразуем коэффициенты по формуле 5
    vec2 Wn     = W(N / numPartitions, k);
    vec4 WHo    = vec4( complex_mult(Wn, Ho.xy), complex_mult(Wn, Ho.zw) ); // vec4(H*W, H*W)
    vec4 result = He + WHo;

    /* Сдвигаем все преобразование на некоторую величину.
     * Данный фрагмент понадобится при реализации синтеза водной поверхности.
     * Просто для вычисления ДПФ он не нужен.
    if ( numPartitions == 1.0 )
    {
        Wn = W(N, -linearIndex * N / 2.0);
        result.xy = complex_mult(Wn, result.xy);
        result.zw = complex_mult(Wn, result.zw);
    }*/

    gl_FragColor = result; // H = He + W * Ho
}

    3.4 Производительность

В следующей таблице приведено сравнение производительности различных реализаций БПФ на GeForce 8800 GTX(Cory Quammen April 11, 2007). Рассмотренная реализация была протестирована на GeForce 8800 GT, при этом она производила два БПФ параллельно, а затраченное время было разделено на два. Видеокарты примерно одного уровня, поэтому результаты включены в ту же таблицу.

1D 16,384
milliseconds
1D 8,388,608
milliseconds
2D 10242
milliseconds
2D 20482
milliseconds
GPUFFTW 169 198 200 288
CUDA 0.701 ? 6.781 32
Mitchell2003 ? ? 12.205 ?
Moreland2003 Cg error Cg error Cg error Cg error
Sumanaweera2005 pbuffer error pbuffer error pbuffer error pbuffer error
Рассмотренная реализация ? ? 12 47
Таблица 2. Сравнение производительности различных реализаций БПФ на GeForce 8800 GTX(GT)

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

Стоит обратить внимание, что узким местом является пропускная способность памяти. При изменении формата текстуры для вычислений с GL_RGBA32F на GL_RGBA16F скорость работы увеличивается практически ровно в два раза. А использование же возможностей GL_EXT_gpu_shader4 не дает никаких преимуществ по производительности.

Исходя из этого, напрашивается оптимизация: переупорядочить выполнение операций таким образом, чтобы повторяющиеся выборки из текстуры находились в кэше. Из [5] видно, что например, при n=0 и n=N/2i выборки будут повторяться. На самом деле ровно половина выборок из текстуры актуальна.

Вероятно, такая оптимизация была сделана в реализации на CUDA, что и привело к двукратному увеличению производительности, но без CUDA такую оптимизацию провести не так просто. Кроме того, при небольших размерах сетки(128x128, 256x256), которые обычно используются при визуализации водной поверхности в реальном времени, вычисление БПФ занимает совсем немного времени.

4. Статистическая модель синтеза поверхности океана

    4.1 Основные положения. Примеры использования

В основе статистической модели синтеза поверхности океана лежит дискретное преобразование Фурье. Волна представляется как сумма большого числа гармоник (т.е. фактически осуществлялось разложение в ряд Фурье по пространственным переменным x и y - [6]).

Проводя ДПФ над комплексной сеткой можно получить значения амплитуд волн в её узлах, то есть, карту высот. Полученная с помощью преобразования Фурье карта высот периодична и её можно "натягивать" на сколь угодно большую поверхность. Соответственно, чем больше размер сетки для преобразования Фурье, тем больше полученная карта высот, тем более изощренная поверхность и менее заметна периодичность волн, что сильно увеличивает реалистичность.




Формула 6. Статистическая модель синтеза океанской поверхности.
L(Lx,Lz) – размер фрагмента водной поверхности. (N,M) – размеры сетки для преобразования Фурье.

Впервые модель была описана Джерри Тессендорфом (Jerry Tessendorf. "Simulating Ocean Water") и использовалась впоследствии почти во всех проектах, требовавших визуализации реалистичной поверхности воды. Один из наиболее полных обзоров реализации этой модели и разнообразных эффектов, таких как каустики и пена, был произведен в работе Lasse Staff Jensen and Robert Golias. "Deep-Water Animation and Rendering".

Выше был описан алгоритм вычисления ДПФ полностью на GPU, поэтому реализацию описанной модели можно так же переложить на графический ускоритель. Первая реализация, полностью производимая на GPU, была предложена в 2005 году: Jason L. Mitchell. "Real-Time Synthesis and Rendering of Ocean Water". Но в тот момент вычислительная мощность 3D ускорителей была ещё недостаточно велика, чтобы модель начала широко применяться в приложениях реального времени. На самом деле, таких приложений и сейчас немного.

  


Рисунок 13. Фильм "Титаник" и игра Crysis.
Пример использования статистической модели для визуализации океана.

    4.2 Спектр Филлипса

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

Обычно коэффициенты ДПФ выбираются исходя из спектра Филлипса - [7]. В работе Тессендорфа есть и другие примеры спектров.





Формула 7. Спектр Филлипса. ξri - выборки из нормального распределения. w – направление и сила ветра. L – размер фрагмента водной поверхности. (N,M) – размер ДПФ. A - нормирующая константа.

Замечание: Жирным выделены векторы, обычным шрифтом – скалярные величины либо длинны соответствующих векторов, символы с крышечкой – нормализованные векторы. Такая мнемоника будет использована и в дальнейшем.

Синтезировать спектр Филлипса каждый кадр не нужно. Можно просто ”анимировать” спектр с фиксированным по времени периодом.



Формула 8. Анимация спектра Филлипса во времени. g – гравитационная постоянная.

Замечание: - сопряжение . – исходные коэффициенты для [6].

Можно заметить, что коэффициенты преобразования Фурье комплексные, хотя высота волны представляется вещественным числом. На самом деле, комплексная компонента не лишняя, она означает фазу волны. Формула [6] как раз и описывает сдвиг фазы волны. Таким образом, сдвиг фаз всех гармоник дает плавную, и в то же время, непредсказуемую, анимацию поверхности.

    4.3 Нормали и "острые" волны

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


Формула 9. Вычисление градиента для водной поверхности. Обозначения такие же, как и в
[4].

Замечание: Для восстановление нормали по градиенту в данном случае можно воспользоваться формулой:
normal = normalize( vec3(-slope.x, 1.0, -slope.y) ).
Выводится из утверждения о том, что нормаль в точке (x,y,z) для функции
F(x,y,z) = 0 совпадает с градиентом в этой точке.

Из-за того, что преобразование Фурье представляет собой сумму гармоник, то есть достаточно гладких функций sin, cos, получить остроконечные штормовые волны трудно. Для этого придется использовать очень большие размеры сеток коэффициентов для преобразования Фурье, то есть очень большое число гармоник. К счастью есть "искусственный" способ придать волнам остроконечный профиль.

Суть метода состоит в сдвиге точек поверхности воды так, чтобы "сжать" волны. Величину сдвига можно определить из следующего преобразования Фурье, которое очень напоминает ДПФ для вычисления градиента.


Формула 10. Вычисление величины сдвига точек поверхности воды, для придания волнам остроконечного профиля.

Замечание: – нормализованный вектор k.

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