Войти
ПрограммированиеСтатьиФизика

Имитация раскачивания дерева под воздействием ветра

Автор:

Статья демонстрирует простой метод имитации воздействия ветра на дерево (куст), реализованный на вертексном шейдере.

tree | Имитация раскачивания дерева под воздействием ветра

К статье прилагаются три примера, содержащие скомпилированные EXE файлы и их исходники на языке Visual Basic из состава Microsoft Visual Studio 2008 Express. Для запуска EXE файлов необходимо иметь установленными NET Framework 2.0, DirectX9.0c и Managed DirectX. Видеоадаптер должен поддерживать вертексные шейдеры vs_1_1.
Примеры следуют один за другим, шаг за шагом описывая построение программы.

Step1

Этот пример подготовительный. В нем мы готовим скелет дерева, с которым будем работать в дальнейшем. То есть мы будем использовать не готовую модель, а только ее имитацию, достаточную для оценки работы физики. Дерево можно создать на основе фрактала по такой схеме:
  1.  Имеем 2 вектора, начальную позицию и направление.
  2.  Позиция задает начало, а сумма позиции и направления конец отрезка, который будет добавлен в нашу модель.
  3.  Считаем новой позицией конец отрезка.
  4.  Трижды вычисляем новое направление, отклоняя старое направление в разные стороны на некоторый угол, длина вектора нового направления должна быть несколько короче старого (вводим коэффициент). Рекурсивно переходим к п.2 для каждой новой пары позиция-направление.
Таким образом, первый отрезок становится стволом дерева, от ствола отходит три ветви, несколько короче самого ствола, от каждой ветви еще три ветви и т. д. Глубину рекурсии необходимо ограничить. На языке Visual Basic это выглядит так:

  Private Sub AddStick(ByVal Level As Integer, ByVal pV As Integer, ByVal Dir As Vector3)
    Dim d, vS, vC As Vector3
    Dim p As Integer, s As Single
    pVert = pVert + 1
    Vert(pVert).Position = Vert(pV).Position + Dir
    Ind(pInd) = pV
    pInd = pInd + 1
    Ind(pInd) = pVert
    pInd = pInd + 1

    If Level = 1 Then Exit Sub

    StickLen = StickLen * 0.7
    d = Vector3.Normalize(Dir) * StickLen
    s = StickLen * 0.6
    vS = Vector3.Normalize(Vector3.Cross(Dir, New Vector3(Rnd, Rnd, Rnd))) * s
    vC = Vector3.Normalize(Vector3.Cross(Dir, vS)) * s

    p = pVert
    AddStick(Level - 1, p, vS * 0.5 + vC * 0.866 + d * (Rnd() + 1.5) * 0.5)
    AddStick(Level - 1, p, vS * -1 + d * (Rnd() + 1.5) * 0.5)
    AddStick(Level - 1, p, vS * 0.5 + vC * -0.866 + d * (Rnd() + 1.5) * 0.5)
    StickLen = StickLen / 0.7
  End Sub

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

vs_1_1

dcl_position v0
def c8, 1, 1, 1, 0

m4x4 oPos, v0, c0
mov oD0, c8

Просто перерасчитываем координаты и выводим в oD0 белый цвет, заданный в константе c8. В константы c0-c3 загружаем матрицу, обобщающую трансформации проекции, видовую и мировую:

    Dim M1 As Matrix
    Dim M2 As Matrix

    M1 = Dev.Transform.World
    M2 = Dev.Transform.View
    M1 = Matrix.Multiply(M1, M2)
    M2 = Dev.Transform.Projection
    M1 = Matrix.Multiply(M1, M2)
    M2 = Matrix.TransposeMatrix(M1)
    Dev.SetVertexShaderConstant(0, M2)

Вот и все, посмотрите пример в папке Step1, в нем двойным щелчком мыши можно генерировать новое дерево, левой кнопкой вращаем камеру вокруг дерева, которое пока статично. Уровень рекурсии 7, что соответствует 1094 вертексов (1093 ветви).

Step2

Теперь немного теории. Зададимся целью вносить в координаты вертексов некоторую регулируемую поправку, имитирующую сгибание ветки от ветра. Можно вычислять матрицу поворота вокруг оси, проходящей через начало ветки и перпендикулярной вектору силы (ветра), это было бы правильно, но очень громоздко. При уровне рекурсии 7 крайние ветви должны испытывать отклонение по семи матрицам! Это тяжелые расчеты, и это сильно «потяжелевший» вертекс. Попробуем найти более простой, пусть и не столь точный метод. Будем просто смещать вершину ветки в нужном направлении, наша задача найти это направление. Ветер не может быть вертикальным, то есть наша задача найти вектор смещения для ветра единичной силы для каждой ветки для двух направлений ветра, – вдоль оси X и вдоль оси Z, – и записать в вертексы эти вектора. Далее в шейдере останется умножить первый вектор на компоненту X ветра, второй на компоненту Z и просуммировать полученные вектора – получим результирующее смещение.
Смещение может быть только перпендикулярно направлению ветки (ведь ветка не может изменять длину), оно должно лежать в одной плоскости с векторами направлений ветки и ветра. Чтобы найти такой вектор, достаточно найти Cross продукт от вектора, перпендикулярного ветке и ветру, и вектора самой ветки. А перпендикулярный вектор так же находится Cross продуктом векторов ветра и ветки:

v = Vector3.Normalize(Vector3.Cross(Dir, Vector3.Cross(New Vector3(1, 0, 0), Dir))) * Dir.Length

Длина вектора смещения должна быть пропорциональна длине ветки, поэтому мы нормализовали результат и умножили на длину ветки.
Это для ветра вдоль оси X. То же самое для ветра вдоль оси Z:
v = Vector3.Normalize(Vector3.Cross(Dir, Vector3.Cross(New Vector3(0, 0, 1), Dir))) * Dir.Length

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

s = Dir.Length * (0.3 + meLevels - Level)

Здесь meLevels –  уровень рекурсии ствола, а Level – уровень рекурсии текущей ветки.
Теперь в вертексе кроме координат содержатся еще два вектора, задающие смещения при ветре вдоль осей X и Z. В вертексный шейдер их можно передать, например, как нормаль и бинормаль (а можно как-нибудь по-другому, вектор – он и есть вектор). Направление ветра передаем в шейдер через константу c4:

Dim w As Vector4
w.X = Wind.X
w.Y = 0
w.Z = Wind.Z
w.W = 0
Dev.SetVertexShaderConstant(4, w)

А сам шейдер выглядит так:

vs_1_1

dcl_position v0
dcl_normal v1
dcl_binormal v2

def c5, 1, 1, 1, 0

mov r0.w, c5.w
mul r0.xyz, v1.xyz, c4.x
add r1, r0, v0

mul r0.xyz, v2.xyz, c4.z
add r1, r1, r0

m4x4 oPos, r1, c0
mov oD0, c5

Тут вычисляется, как уже было сказано чуть выше, сумма произведений векторов-коэффициентов на соответствующие  компоненты вектора ветра, компонента w приравнивается нулю, и само смещение суммируется с координатами вершины. В примере Step2, как и в предыдущем, левой кнопкой вращаем камеру, правой направление ветра, а клавишами  1 – 9 задаем силу ветра. Мы видим сгибающееся в разных направлениях дерево, но на ветер это пока не похоже – не хватает хаотичности. При смене силы ветра дерево меняет вид мгновенно – то есть не хватает еще и инерции.

Step3

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

  Private Function GetNoise(ByVal s As Single) As Single
    Dim n1, n2 As Integer
    Dim t As Single
    n1 = Int(s)
    t = s - n1
    n1 = n1 And 255
    n2 = (n1 + 1) And 255
    GetNoise = Noise(n1) * (1 - t) + Noise(n2) * t
  End Function

  Private Sub GenNoise()
    Dim n As Integer, dn As Integer
    dn = 32
    Do
      For n = dn To 255 Step dn * 2
        Noise(n) = (Noise(n - dn) + Noise((n + dn) And 255)) * 0.5 + (Rnd() - 0.5) * dn / 16
      Next n
      dn = dn >> 1
      If dn = 0 Then Exit Do
    Loop
  End Sub

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

sv4 = sv4 * (1 - dt) + dv4 * dt

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

dv4.X = v.X * (1 - 0.5 * GetNoise(t * 63))
dv4.Y = 0
dv4.Z = v.Z * (1 - 0.5 * GetNoise(t * 50))
dv4.W = 0

Константы 63 и 50 подобраны имперически, это скорость движения по функции шума компонент X и Z, они заданы разными, таким образом, ветер будет менять не только силу, но и направление.
Уже стало реальнее, но добавим последний штрих – резонансные случайные движения тонких веток, зависящие от силы ветра. Добавим в вертекс три коэффициента, вместе они прекрасно уместятся в незанятый еще Tangent. Для трех верхних ступеней рекурсии заполним их случайными числами так, чтобы для тонких веток значения были больше:

    If Level = 1 Then
      s = Dir.Length * 0.3
    ElseIf Level = 2 Then
      s = Dir.Length * 0.15
    ElseIf Level = 3 Then
      s = Dir.Length * 0.075
    Else
      s = 0
    End If
    v = New Vector3(s * (Rnd() - 0.5), s * (Rnd() - 0.5), s * (Rnd() - 0.5))

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

    v4.X = GetNoise(QTime() * 33) * 10 * s
    v4.Y = GetNoise(QTime() * 52) * 10 * s
    v4.Z = GetNoise(QTime() * 94) * 10 * s
    Dev.SetVertexShaderConstant(5, v4)

Константы 33, 52, 94 определяют скорость дрожания веток, можно подобрать другие величины на свой вкус. Константа 10 определяет амплитуду колебаний, а s – сила ветра (длина вектора).
Мы не зря записали в вертексы случайные коэффициенты, на колебания разных веток выборки шума будут сказываться в разной степени, ветки будут двигаться не синхронно. Наша задача – в шейдере перемножить соответствующие коэффициенты на выборки шума и просуммировать эти величины. И с этой задачей прекрасно справится единственная инструкция dp3! Итак, наш шейдер:

vs_1_1

dcl_position v0
dcl_normal v1
dcl_binormal v2
dcl_tangent v3

def c6, 1, 1, 1, 0

dp3 r2, v3, c5

add r3, r2, c4.xxxw
mul r0, v1, r3
add r1, r0, v0

add r3, r2, c4.zzzw
mul r0, v2, r3
add r1, r1, r0

m4x4 oPos, r1, c0
mov oD0, c6

В примере Step3 глубина рекурсии увеличена до 8 (3281 ветка), управление такое же, как и в Step2. Реализованы инерция, неравномерность и случайные колебания тонких веток.
Статья окончена, в качестве творческого задания предлагаю разработать алгоритм записи таких же коэффициентов в уже готовую, созданную моделлерами, модель дерева.

Примеры к статье: Дерево на ветру

#шейдеры, #Visual Basic

9 октября 2020

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