NoiseSampler.cpp

#include "NoiseSampler.h"
 
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <algorithm>
 
#include <glm/vec2.hpp>
#include <glm/vec4.hpp>
 
#include <random>
 
namespace
{
  char grad2D[4][2] =
  {
    { 1, 1 },
    { -1, 1 },
    { 1, -1 },
    { -1, -1 }
  };
}
 
namespace Cogs
{
  namespace Core
  {
    NoiseSampler::NoiseSampler(unsigned int seed)
    {
      //if (permutationTable.empty()) {
 
      // Create permuation table
      std::minstd_rand eng(seed);
      std::uniform_int_distribution<int> dist(0, RAND_MAX);
 
      permutationTable.resize(256);
      for (unsigned i = 0; i < 256; i++) {
        permutationTable[i] = static_cast<uint8_t>(i);
      }
      for (int i = 255; i > 0; i--) {
        int j = dist(eng) / ((RAND_MAX / i) + 1);  // use higher-order bits of rand
        std::swap(permutationTable[j], permutationTable[i]);
      }
 
      // create gradient table
      gradient2DTable.resize(256 * 2);
      for (int i = 0; i < 256; i++) {
        float theta = dist(eng) * (float((2.0f * glm::pi<float>()) / float(RAND_MAX))) - 1.f;
        gradient2DTable[2 * i + 0] = std::cos(theta);
        gradient2DTable[2 * i + 1] = std::sin(theta);
      }
      //}
    }
 
    float NoiseSampler::gradNoise2D(glm::vec2 t, glm::ivec2 period)
    {
      glm::vec2 ti = glm::floor(t); // integer part
 
      glm::ivec2 tl(int(ti.x) % period.x, int(ti.y) % period.y);
      glm::ivec2 tq((tl.x + 1) % period.x, (tl.y + 1) % period.y);
 
      glm::vec2 tf = t - ti;        // fractional part
      glm::vec2 sf = glm::vec2(1.f) - tf;
 
      // we do a reparameterization of the fractional part, to improve smoothness
      // over lattice lines. In Perlin's original 1984 paper, the
      // reparameterization was
      //
      //   phi   = 3t^2-2t^3,    phi(0)   = 0, phi(1)   =  1
      //   phi'  = 6t - 6t^2,    phi'(0)  = 0, phi'(1)  =  0,
      //   phi'' = 6 - 12t,      phi''(0) = 6, phi''(1) = -6
      //
      // which has nonzero second-order derivatives at 0 and 1. This becomes
      // visible in specular highlights, and in "Improving Noise", Perlin suggests
      // to use
      //
      //   phi   = 6t^5 - 15t^4 + 10t^3,   phi(0)   = 0, phi(1)   = 1
      //   phi'  = 30t^4 - 60t^3 + 30t^2,  phi'(0)  = 0, phi'(1)  = 0,
      //   phi'' = 120t^3 - 180t^2 + 60t,  phi'(0)  = 0, phi'(1)  = 0,
      //
      // which indeed has zero second-order derivatives.
      glm::vec2 phi_t = tf*tf*tf*(6.f*tf*tf - 15.f*tf + 10.f);
      glm::vec2 phi_s = glm::vec2(1.f) - phi_t;
 
      int l0x = permutationTable[tl.x & 0xff];
      int l00 = (permutationTable[(l0x + tl.y)&0xff]);
      int l01 = (permutationTable[(l0x + tq.y)&0xff]);
      int l1x = permutationTable[tq.x & 0xff];
      int l10 = (permutationTable[(l1x + tl.y)&0xff]);
      int l11 = (permutationTable[(l1x + tq.y)&0xff]);
 
#if 1
      float vy0 = (grad2D[l00&0x3][0] * tf.x + grad2D[l00&0x3][1] * tf.y)*phi_s.x
        + (grad2D[l10&0x3][0] * sf.x + grad2D[l10&0x3][1] * tf.y)*phi_t.x;
      float vy1 = (grad2D[l01&0x3][0] * tf.x + grad2D[l01&0x3][1] * sf.y)*phi_s.x
        + (grad2D[l11&0x3][0] * sf.x + grad2D[l11&0x3][1] * sf.y)*phi_t.x;
#else
      float vy0 = (gradient2DTable[2 * l00 + 0] * tf.x + gradient2DTable[2 * l00 + 1] * tf.y)*phi_s.x
        + (gradient2DTable[2 * l10 + 0] * sf.x + gradient2DTable[2 * l10 + 1] * tf.y)*phi_t.x;
 
      float vy1 = (gradient2DTable[2 * l01 + 0] * tf.x + gradient2DTable[2 * l01 + 1] * sf.y)*phi_s.x
        + (gradient2DTable[2 * l11 + 0] * sf.x + gradient2DTable[2 * l11 + 1] * sf.y)*phi_t.x;
#endif
 
      return vy0*phi_s.y + vy1*phi_t.y;
    }
 
    void NoiseSampler::gradNoise2D(float* dst, float frequency, int w, int h)
    {
      glm::ivec2 period(w/frequency, h/frequency);
 
      float a = 1.f / std::max(w, h);
      for (int j = 0; j < h; j++) {
        float v = a*static_cast<float>(j);
        for (int i = 0; i < w; i++) {
          float u = a* static_cast<float>(i);
          dst[j*w + i] = gradNoise2D(frequency*glm::vec2(u+0.5f, v+0.5f), period);
        }
      }
    }
 
    void NoiseSampler::gradNoise2DSum(float* dst, float frequency_min, float frequency_max, int w, int h)
    {
      std::vector<glm::vec4> param;   // shift u, shift v, frequency, weight
      std::vector<glm::ivec2> periods;
      for (float f = frequency_min; f < frequency_max; f = f * 2) {
        param.push_back(glm::vec4(1.f / f, 1.f / f, f, 2.f*frequency_min / f));
        periods.push_back(glm::ivec2(w / f, h / f));
      }
 
      float a = 1.f / std::max(w, h);
 
      for (int j = 0; j < h; j++) {
        float v = a*static_cast<float>(j);
        for (int i = 0; i < w; i++) {
          float u = a*static_cast<float>(i);
 
          float sum = 0.f;
          for (size_t k = 0; k < periods.size(); k++) {
            sum += param[k].w*gradNoise2D(param[k].z*glm::vec2(u + 0.5f, v + 0.5f), periods[k]);
          }
          dst[j*w + i] = sum;
        }
      }
    }
 
    void NoiseSampler::gradNoise2DTurbulence(float* dst, float frequency_min, float frequency_max, int w, int h)
    {
      std::vector<glm::vec4> param;   // shift u, shift v, frequency, weight
      std::vector<glm::ivec2> periods;
      for (float f = frequency_min; f < frequency_max; f = f * 2) {
        param.push_back(glm::vec4(1.f / f, 1.f / f, f, 2.f*frequency_min / f));
        periods.push_back(glm::ivec2(w / f, h / f));
      }
 
 
      float a = 1.f / std::max(w, h);
 
      for (int j = 0; j < h; j++) {
        float v = a*static_cast<float>(j);
        for (int i = 0; i < w; i++) {
          float u = a*static_cast<float>(i);
 
          float sum = 0.f;
          for (size_t k = 0; k < periods.size(); k++) {
            sum += std::abs(param[k].w*gradNoise2D(param[k].z*glm::vec2(u+0.5f, v+0.5f), periods[k]));
          }
          dst[j*w + i] = sum - 1.f;
        }
      }
    }
  } // namespace Core
} // namespace Cogs