WaveSpectrum.cpp

#include "ClipmapTerrainTypes.h"
#include "RenderContext.h"
#include "WaveSpectrum.h"
 
#include "Rendering/ITextures.h"
#include "Rendering/IRenderTargets.h"
#include "Rendering/IGraphicsDevice.h"
#include "Rendering/IBuffers.h"
#include "Rendering/CommandGroupAnnotation.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cassert>
 
using std::complex;
using std::sqrt;
using std::min;
using std::max;
using std::pow;
using std::exp;
using std::tgamma;
using std::sin;
using std::cos;
using std::acos;
using std::log;
using std::polar;
using std::isfinite;
 
using glm::vec2;
using glm::vec3;
using glm::vec4;
using glm::ivec2;
using glm::length;
using glm::dot;
using glm::cross;
 
namespace {
 
  const float pi = float(M_PI);
  const float twoPi = float(2.0*M_PI);
  const float g = 9.80665f; // gravity of earth
 
  struct DispersionParameters{
    int N;
    float twoPiOverSideLength;
    float dt;
  };
 
  struct PackParamters{
    int N_minus_1;
    float L_over_twoN;
  };
 
  float waveSpectrumPhillips(const float omega,
                             const float alpha = 0.0081f)
  {
    const float omega2 = omega*omega;
    const float omega4 = omega2*omega2;
    const float omega5 = omega4*omega;
 
    return (alpha*g*g) / omega5;
  }
 
  float waveSpectrumPiersonMoskowitz(const float omega,
                                     const float omega_p, // frequency of the spectrum peak
                                     const float alpha = 0.0081f)
  {
    const float omega_p_over_omega = omega_p / omega;
    const float omega_p_over_omega4 = (omega_p_over_omega*omega_p_over_omega)*(omega_p_over_omega*omega_p_over_omega);
 
    return waveSpectrumPhillips(omega, alpha)*exp(float(-5.0 / 4.0)*omega_p_over_omega4);
  }
 
  float dispersionLonguetHiggins(const float omega,
                                 const float omega_p,
                                 const float U10,
                                 const float windWaveAngle) // in [-pi,pi]
  {
 
    // Sharpness of directional spreading (Mitsuyasu et al., 1975)
    float mu = omega <= omega_p ? 5.f : -2.5f;
    float s_omega = 11.5f*pow(g / (omega_p*U10), 2.5f)*pow(omega / omega_p, mu);
    assert(isfinite(s_omega));
 
    // Normalization factor of dispersion relation
 
    // Part of this is the ratio gamma(x+1)/gamma(x+0.5). Evaluating this
    // directly is begging for overflow since gamma rapidly becomes
    // ridiculously large, e.g. gamma(30)=8.8x10^30. However, since we are
    // interested in the ratio, the following asymptotic series,
    //
    // gamma(J+1/2)/gamma(J) = sqrt(gamma)*(1 - 1/8J + 1/128J^2
    //                                      5/1024J^3 - 21/32768J^4 + ... )
    //
    // is handy, as if we define J = s_omega + 0.5, we can use it directly
    // to evaluate the ratio.
    float J = s_omega + 0.5f;
    float gammaRatio = sqrt(J)*(1.f - 1.f / (8.f*J) + 1.f / (128.f*J*J) + 5.f / (1024.f*J*J*J) - 21.f / (32768.f*J*J*J*J));
    float N_s_omega = float(1.0 / (2.0*sqrt(M_PI)))*(gammaRatio);
    assert(isfinite(N_s_omega));
 
    //    x!(x / ((1 / 2 (2 x + 1))!) + 1 / (2 (1 / 2 (2 x + 1))!))
 
    // Wind-wave direction angle halved
    //    float theta_half = 0.5f*(waveDirection);
 
    // Dispersion relation (Longuet-Higgins, 1962)
    float D = N_s_omega*pow(max(0.0f, cos(0.5f*windWaveAngle)), 2.f*s_omega);
    assert(isfinite(D));
    return D;
  }
 
} // of anonymous namespace
 
float Cogs::WaveSpectrum::createDirectionalWaveSpectrum(std::vector<float>& E,
                                                        int N,
                                                        float L,
                                                        float waveNumberMute,
                                                        float waveNumberPass,
                                                        float windSpeed,
                                                        float windDirection,
                                                        float dominantWavePeriod,
                                                        float scale)
{
 
  if (windDirection >= pi) {
    windDirection -= twoPi;
  }
  if (windSpeed < 2.f) {
    windSpeed = 2.f;
  }
 
  float dominantAngularVelocity = float(2.0*M_PI) / dominantWavePeriod;
  float sumS = 0.f;
 
  glm::vec2 windDir(glm::cos(windDirection), glm::sin(windDirection));
 
  for (int j = 0; j < N; j++) {
    for (int i = 0; i < N; i++) {
      if ((i == 0) && (j == 0)) {
        continue;
      }
      // [0,N/2]     -> [0,N/2]
      // [N/2+1,N-1] -> [-N/2+1,-1]
      ivec2 ij(i <= N / 2 ? i : i - N,
               j <= N / 2 ? j : j - N);
 
      const vec2 K = (twoPi / L)*vec2(ij);  // K:     Wave vector
      float k = length(K);                  // k:     Wave number
 
      float angularVelocity = sqrt(g*k);    // omega: Angular velocity
      //float waveDirection = std::atan2(K.y, K.x); // theta_w: Wave direction
 
      float cosWindWaveAngle = std::max(-1.f, std::min(1.f, (1.f / k)*dot(K, windDir)));
      float windWaveAngle = std::acos(cosWindWaveAngle);
      assert(std::isfinite(windWaveAngle));
 
      float pass = 1.f;
      if (waveNumberMute != waveNumberPass) {
        pass = max(0.f, min(1.f, (k - waveNumberMute) / (waveNumberPass - waveNumberMute)));
      }
 
      const float S = waveSpectrumPiersonMoskowitz(angularVelocity,
                                                   dominantAngularVelocity,
                                                   0.0081f);
      const float D = dispersionLonguetHiggins(angularVelocity,
                                               dominantAngularVelocity,
                                               windSpeed,
                                               windWaveAngle);
 
      const float chainFactor = (1.f / (2.f*k))*sqrt(g / k);
 
      sumS += chainFactor*pass*S;
      E[N*j + i] = chainFactor*pass*S*D;
    }
  }
  E[0] = 0.f;
 
  // Adjust spectrum to match significant wave height
  float w = 1.f / sumS;// (significantWaveHeight*significantWaveHeight) / (sumS);
  if (std::numeric_limits<float>::epsilon() < std::abs(scale)) {
    w = scale;
  }
  for (int i = 0; i < N*N; i++) {
    E[i] = w*E[i];
  }
  return 1.f / sumS;
}
#define MINSTD_RAND_MAX ((1u<<31)-2u)
static uint32_t minstd_rand(uint32_t &seed)
{
  seed = ((uint64_t)seed * 48271u) % ((1u<<31)-1u);
  return seed;
}
void Cogs::WaveSpectrum::createRandomizedWaveSpectrumInstance(std::vector<glm::vec2>& H,
                                                              const std::vector<float>& E,
                                                              uint32_t seed,
                                                              const size_t N)
{
  for (size_t j = 0; j < N; j++) {
    for (size_t i = 0; i < N; i++) {
 
      // Check the sanity of this trying to comply with normal distribution
 
      float U0 = (float)((double)(minstd_rand(seed) + 1) / (double)(MINSTD_RAND_MAX + 2));
      float U1 = (float)((2.0*M_PI*(double)minstd_rand(seed)) / (double)(MINSTD_RAND_MAX + 1));
 
      // Box-Muller transform to create normal distributed numbers, mean 0, var 1.
      complex<float> eta = sqrt(-2.f*log(U0))*complex<float>(cos(U1), sin(U1));
 
      complex<float> res = sqrt(E[N*j + i] / 2.f) * eta;
 
      H[N*j + i] = vec2(res.real(), res.imag());
    }
  }
}
 
float Cogs::WaveSpectrum::setConditions(const float tileExtent,
                                        const float waveNumberMute,
                                        const float waveNumberPass,
                                        const float significantWavePeriod,
                                        const float windSpeed,
                                        const float /*windDirection*/,
                                        const float scale,
                                        const unsigned int /*seed*/)
{
  this->tileExtent = tileExtent;
 
  float significantWaveLength = (g*significantWavePeriod*significantWavePeriod) / (2.f*glm::pi<float>());
 
  int harmonic = std::max(1, static_cast<int>(std::round(std::log2(tileExtent / significantWaveLength))));
 
  tileExtentAdjust = ((1 << harmonic)*significantWaveLength) / tileExtent;
 
  float rv = createDirectionalWaveSpectrum(frequencyDomain.E,
                                           N,
                                           tileExtentAdjust*tileExtent,
                                           waveNumberMute,
                                           waveNumberPass,
                                           windSpeed,
                                           0.f,
                                           significantWavePeriod,
                                           scale);
  createRandomizedWaveSpectrumInstance(frequencyDomain.a, frequencyDomain.E, 42, N);
 
  ITextures* textures = device->getTextures();
  frequencyDomain.aTex = textures->loadTexture(reinterpret_cast<unsigned char*>(frequencyDomain.a.data()), N, N, TextureFormat::R32G32_FLOAT);
  textures->annotate(frequencyDomain.aTex, "Initial spectrum instance.");
 
  return rv;
}
 
void Cogs::WaveSpectrum::initialize(IGraphicsDevice* device)
{
  this->frame = 0;
  this->device = device;
 
  gpgpuQuadRenderer.initialize(device);
  fourierTransform.initialize(device, gpgpuQuadRenderer);
 
  auto ie = device->getEffects();
 
  {
    PreprocessorDefinitions definitions;
    //  definitions.push_back(PreprocessorDefinition("TWO_PI_OVER_L", float(2.0*M_PI/L)));
 
    disperse.effect = ie->loadEffect("Terrain/GPGPUPassThroughVS.hlsl",
                                     "Terrain/WaveSpectrumDispersionPS.hlsl",
                                     definitions);
    assert(ie->checkEffect(disperse.effect) == Cogs::ResourceStatus::Ready && "Expects synchronous effect loading");
    VertexFormatHandle handle = gpgpuQuadRenderer.vertexFormat();
    disperse.il = device->getBuffers()->loadInputLayout(&handle, 1, disperse.effect);
    disperse.constantBuffer = device->getBuffers()->loadBuffer(nullptr, sizeof(DispersionParameters), Usage::Dynamic, AccessMode::Write, BindFlags::ConstantBuffer);
  }
 
  {
    PreprocessorDefinitions definitions;
 
    texturePack.effect = ie->loadEffect("Terrain/GPGPUPassThroughVS.hlsl",
                                        "Terrain/OceanBuildTexPositionPS.hlsl",
                                        definitions);
    assert(ie->checkEffect(texturePack.effect) == Cogs::ResourceStatus::Ready && "Expects synchronous effect loading");
    VertexFormatHandle handle = gpgpuQuadRenderer.vertexFormat();
    texturePack.il = device->getBuffers()->loadInputLayout(&handle, 1, texturePack.effect);
    texturePack.constantBuffer = device->getBuffers()->loadBuffer(nullptr, sizeof(PackParamters), Usage::Dynamic, AccessMode::Write, BindFlags::ConstantBuffer);
  }
}
 
void Cogs::WaveSpectrum::setSize(const int NLog2)
{
  fourierTransform.setSize(NLog2);
 
  this->NLog2 = NLog2;
  N = 1u << NLog2;
 
  std::vector<float> zeros(N*N);
  frequencyDomain.E.resize(N*N);
  frequencyDomain.a.resize(N*N);
 
  ITextures* textures = device->getTextures();
  IRenderTargets* renderTargets = device->getRenderTargets();
 
  packed.xyzTex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32B32A32_FLOAT, TextureFlags::RenderTarget | TextureFlags::GenerateMipMaps);
  packed.dxdu_dydv_dzdu_dzdvTex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32B32A32_FLOAT, TextureFlags::RenderTarget | TextureFlags::GenerateMipMaps);
 
  for (int i = 0; i < 2; i++) {
    phaseTex[i] = textures->loadTexture(reinterpret_cast<uint8_t*>(zeros.data()), N, N, TextureFormat::R32_FLOAT, TextureFlags::RenderTarget);
  }
 
  frequencyDomain.xyTex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32B32A32_FLOAT, TextureFlags::RenderTarget);
  frequencyDomain.zTex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32_FLOAT, TextureFlags::RenderTarget);
  frequencyDomain.dzdu_dzdv_Tex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32B32A32_FLOAT, TextureFlags::RenderTarget);
 
  spatialDomain.xyTex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32B32A32_FLOAT, TextureFlags::RenderTarget);
  spatialDomain.zTex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32_FLOAT, TextureFlags::RenderTarget);
  spatialDomain.dzdu_dzdv_Tex = textures->loadTexture(nullptr, N, N, TextureFormat::R32G32B32A32_FLOAT, TextureFlags::RenderTarget);
 
  // Create render targets
  spatialDomain.xyTarget = renderTargets->createRenderTarget(spatialDomain.xyTex);
  spatialDomain.zTarget = renderTargets->createRenderTarget(spatialDomain.zTex);
  spatialDomain.dzduTarget = renderTargets->createRenderTarget(spatialDomain.dzdu_dzdv_Tex);
 
  for (int i = 0; i < 2; i++) {
    TextureHandle texs[4] = {
      phaseTex[i],
      frequencyDomain.xyTex,
      frequencyDomain.zTex,
      frequencyDomain.dzdu_dzdv_Tex
    };
    dispersionTarget[i] = renderTargets->createRenderTarget(texs, 4);
  }
 
  TextureHandle packedTexs[2] = {
    packed.xyzTex ,
    packed.dxdu_dydv_dzdu_dzdvTex
  };
  packed.packTarget = renderTargets->createRenderTarget(packedTexs, 2);
 
}
 
bool Cogs::WaveSpectrum::update(RenderContext& renderContext, const float dt)
{
  IContext* context = renderContext.context;
 
  {
    CommandGroupAnnotation preGroup(renderContext.context, "WaveSpectrum::Disperse");
 
 
    context->setRenderTarget(dispersionTarget[frame], DepthStencilHandle::InvalidHandle);
    context->setViewport(0, 0, float(N), float(N));
    context->setEffect(disperse.effect);
    context->setInputLayout(disperse.il);
 
    gpgpuQuadRenderer.bind(context);
 
    context->setTexture("aTex", 0, frequencyDomain.aTex);
    context->setTexture("phaseTex", 0, phaseTex[(frame + 1) & 1]);
 
    {
      MappedBuffer<DispersionParameters> constants(context, disperse.constantBuffer, MapMode::WriteDiscard);
      if (constants) {
        constants->N = N;
        constants->twoPiOverSideLength = float(2.0*M_PI / (tileExtentAdjust*tileExtent));
        constants->dt = dt;
      }
    }
    context->setConstantBuffer("DispersionParameters", disperse.constantBuffer);
 
    gpgpuQuadRenderer.draw(context);
  }
 
  {
    CommandGroupAnnotation preGroup(renderContext.context, "WaveSpectrum::iFFT");
 
    fourierTransform.inverseFourierTransform(renderContext, gpgpuQuadRenderer, spatialDomain.xyTarget, frequencyDomain.xyTex, true);
 
    fourierTransform.inverseFourierTransform(renderContext, gpgpuQuadRenderer, spatialDomain.zTarget, frequencyDomain.zTex, false);
 
    fourierTransform.inverseFourierTransform(renderContext, gpgpuQuadRenderer, spatialDomain.dzduTarget, frequencyDomain.dzdu_dzdv_Tex, true);
  }
 
  {
    CommandGroupAnnotation preGroup(renderContext.context, "WaveSpectrum::Pack");
 
    context->setRenderTarget(packed.packTarget, DepthStencilHandle::InvalidHandle);
    context->setViewport(0, 0, float(N), float(N));
 
    context->setEffect(texturePack.effect);
    {
      MappedBuffer<PackParamters> constants(context, texturePack.constantBuffer, MapMode::WriteDiscard);
      if (constants) {
        constants->N_minus_1 = N - 1;
        constants->L_over_twoN = tileExtent / (2.f*N);
      }
    }
    context->setConstantBuffer("PackParamters", texturePack.constantBuffer);
    context->setTexture("waveXY", 0, spatialDomain.xyTex);
    context->setTexture("waveZ", 1, spatialDomain.zTex);
    context->setTexture("wavedZdu_dZdv", 2, spatialDomain.dzdu_dzdv_Tex);
 
    gpgpuQuadRenderer.bind(context);
    context->setInputLayout(texturePack.il);
    gpgpuQuadRenderer.draw(context);
 
    auto textures = device->getTextures();
    textures->generateMipmaps(packed.xyzTex);
    textures->generateMipmaps(packed.dxdu_dydv_dzdu_dzdvTex);
  }
 
  frame = (frame + 1) & 1;
 
  return true;
}