WAVELET_NOISE.h

Go to the documentation of this file.

// This file is part of Wavelet Turbulence.
// 
// Wavelet Turbulence is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// Wavelet Turbulence is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
// 
// Copyright 2008 Theodore Kim and Nils Thuerey
// 
// Wavelet noise functions
//
// This code is based on the C code provided in the appendices of:
//
// @article{1073264,
//  author = {Robert L. Cook and Tony DeRose},
//  title = {Wavelet noise},
//  journal = {ACM Trans. Graph.},
//  volume = {24},
//  number = {3},
//  year = {2005},
//  issn = {0730-0301},
//  pages = {803--811},
//  doi = {http://doi.acm.org/10.1145/1073204.1073264},
//  publisher = {ACM},
//  address = {New York, NY, USA},
//  }
//

#ifndef WAVELET_NOISE_H
#define WAVELET_NOISE_H

#include <MERSENNETWISTER.h>

#ifdef WIN32
#include <float.h>
#define isnan _isnan
#endif

// Tile file header, update revision upon any change done to the noise generator
static const char tilefile_headerstring[] = "Noise Tile File rev. ";
static const char tilefile_revision[] =  "001";

#define NOISE_TILE_SIZE 128
static const int noiseTileSize = NOISE_TILE_SIZE;

// warning - noiseTileSize has to be 128^3!
#define modFast128(x) ((x) & 127)
#define modFast64(x)  ((x) & 63)
#define DOWNCOEFFS 0.000334f,-0.001528f, 0.000410f, 0.003545f,-0.000938f,-0.008233f, 0.002172f, 0.019120f, \
                  -0.005040f,-0.044412f, 0.011655f, 0.103311f,-0.025936f,-0.243780f, 0.033979f, 0.655340f, \
                   0.655340f, 0.033979f,-0.243780f,-0.025936f, 0.103311f, 0.011655f,-0.044412f,-0.005040f, \
                   0.019120f, 0.002172f,-0.008233f,-0.000938f, 0.003546f, 0.000410f,-0.001528f, 0.000334f

// Wavelet downsampling -- periodic boundary conditions
static void downsampleX(float *from, float *to, int n){
  // if these values are not local incorrect results are generated
  float downCoeffs[32] = { DOWNCOEFFS };
    const float *a = &downCoeffs[16];
    for (int i = 0; i < n / 2; i++) {
        to[i] = 0;
        for (int k = 2 * i - 16; k < 2 * i + 16; k++)
            to[i] += a[k - 2 * i] * from[modFast128(k)];
    }
}
static void downsampleY(float *from, float *to, int n){
  // if these values are not local incorrect results are generated
  float downCoeffs[32] = { DOWNCOEFFS };
    const float *a = &downCoeffs[16];
    for (int i = 0; i < n / 2; i++) {
        to[i * n] = 0;
        for (int k = 2 * i - 16; k < 2 * i + 16; k++)
            to[i * n] += a[k - 2 * i] * from[modFast128(k) * n];
    }
}
static void downsampleZ(float *from, float *to, int n){
  // if these values are not local incorrect results are generated
  float downCoeffs[32] = { DOWNCOEFFS };
    const float *a = &downCoeffs[16];
    for (int i = 0; i < n / 2; i++) {
        to[i * n * n] = 0;
        for (int k = 2 * i - 16; k < 2 * i + 16; k++)
            to[i * n * n] += a[k - 2 * i] * from[modFast128(k) * n * n];
    }
}

// Wavelet downsampling -- Neumann boundary conditions
static void downsampleNeumann(const float *from, float *to, int n, int stride)
{
  // if these values are not local incorrect results are generated
  float downCoeffs[32] = { DOWNCOEFFS };
  const float *const aCoCenter= &downCoeffs[16];
    for (int i = 0; i < n / 2; i++) {
        to[i * stride] = 0;
        for (int k = 2 * i - 16; k < 2 * i + 16; k++) { 
            // handle boundary
            float fromval; 
            if (k < 0) {
                fromval = from[0];
            } else if(k > n - 1) {
                fromval = from[(n - 1) * stride];
            } else {
                fromval = from[k * stride]; 
            } 
            to[i * stride] += aCoCenter[k - 2 * i] * fromval; 
        }
    }
}
static void downsampleXNeumann(float* to, const float* from, int sx,int sy, int sz) {
    for (int iy = 0; iy < sy; iy++) 
        for (int iz = 0; iz < sz; iz++) {
            const int i = iy * sx + iz*sx*sy;
            downsampleNeumann(&from[i], &to[i], sx, 1);
        }
}
static void downsampleYNeumann(float* to, const float* from, int sx,int sy, int sz) {
    for (int ix = 0; ix < sx; ix++) 
        for (int iz = 0; iz < sz; iz++) {
            const int i = ix + iz*sx*sy;
            downsampleNeumann(&from[i], &to[i], sy, sx);
    }
}
static void downsampleZNeumann(float* to, const float* from, int sx,int sy, int sz) {
    for (int ix = 0; ix < sx; ix++) 
        for (int iy = 0; iy < sy; iy++) {
            const int i = ix + iy*sx;
            downsampleNeumann(&from[i], &to[i], sz, sx*sy);
    }
}

// Wavelet upsampling - periodic boundary conditions
static float _upCoeffs[4] = {0.25f, 0.75f, 0.75f, 0.25f};
static void upsampleX(float *from, float *to, int n) {
    const float *p = &_upCoeffs[2];

    for (int i = 0; i < n; i++) {
        to[i] = 0;
        for (int k = i / 2; k <= i / 2 + 1; k++)
            to[i] += p[i - 2 * k] * from[modFast64(k)];
    }
}
static void upsampleY(float *from, float *to, int n) {
    const float *p = &_upCoeffs[2];

    for (int i = 0; i < n; i++) {
        to[i * n] = 0;
        for (int k = i / 2; k <= i / 2 + 1; k++)
            to[i * n] += p[i - 2 * k] * from[modFast64(k) * n];
    }
}
static void upsampleZ(float *from, float *to, int n) {
    const float *p = &_upCoeffs[2];

    for (int i = 0; i < n; i++) {
        to[i * n * n] = 0;
        for (int k = i / 2; k <= i / 2 + 1; k++)
            to[i * n * n] += p[i - 2 * k] * from[modFast64(k) * n * n];
    }
}

// Wavelet upsampling - Neumann boundary conditions
static void upsampleNeumann(const float *from, float *to, int n, int stride) {
  static const float *const pCoCenter = &_upCoeffs[2];
    for (int i = 0; i < n; i++) {
        to[i * stride] = 0;
        for (int k = i / 2; k <= i / 2 + 1; k++) {
            float fromval;
            if(k>n/2) {
                fromval = from[(n/2) * stride];
            } else {
                fromval = from[k * stride]; 
            }  
            to[i * stride] += pCoCenter[i - 2 * k] * fromval; 
        }
    }
}
static void upsampleXNeumann(float* to, const float* from, int sx, int sy, int sz) {
    for (int iy = 0; iy < sy; iy++) 
        for (int iz = 0; iz < sz; iz++) {
            const int i = iy * sx + iz*sx*sy;
            upsampleNeumann(&from[i], &to[i], sx, 1);
        }
}
static void upsampleYNeumann(float* to, const float* from, int sx, int sy, int sz) {
    for (int ix = 0; ix < sx; ix++) 
        for (int iz = 0; iz < sz; iz++) {
            const int i = ix + iz*sx*sy;
            upsampleNeumann(&from[i], &to[i], sy, sx);
        }
}
static void upsampleZNeumann(float* to, const float* from, int sx, int sy, int sz) {
    for (int ix = 0; ix < sx; ix++) 
        for (int iy = 0; iy < sy; iy++) {
            const int i = ix + iy*sx;
            upsampleNeumann(&from[i], &to[i], sz, sx*sy);
        }
}


// load in an existing noise tile
static bool loadTile(float* const noiseTileData, std::string filename)
{
    FILE* file;
    char headerbuffer[64];
    size_t headerlen;
    size_t bread;
    int endiantest = 1;
    char endianness;
    
    file = fopen(filename.c_str(), "rb");

    if (file == NULL) {
        printf("loadTile: No noise tile '%s' found.\n", filename.c_str());
        return false;
    }

    //Check header
    headerlen = strlen(tilefile_headerstring) + strlen(tilefile_revision) + 2;
    bread = fread((void*)headerbuffer, 1, headerlen, file);
    if (*((unsigned char*)&endiantest) == 1)
        endianness = 'L';
    else
        endianness = 'B';
    if ((bread != headerlen)
        || (strncmp(headerbuffer, tilefile_headerstring, strlen(tilefile_headerstring)))
        || (strncmp(headerbuffer+ strlen(tilefile_headerstring), tilefile_revision, strlen(tilefile_revision)))
        || (headerbuffer[headerlen-2] != endianness)
        || (headerbuffer[headerlen-1] != (char)((char)sizeof(long)+'0')))
    {
        printf("loadTile : Noise tile '%s' was generated on an incompatible platform.\n",filename.c_str());
        fclose(file);
        return false;
    }
    
    // dimensions
    size_t gridSize = noiseTileSize * noiseTileSize * noiseTileSize;

    // noiseTileData memory is managed by caller
    bread = fread((void*)noiseTileData, sizeof(float), gridSize, file);
    fclose(file);
    printf("Noise tile file '%s' loaded.\n", filename.c_str());

    if (bread != gridSize) {
        printf("loadTile: Noise tile '%s' is wrong size %d.\n", filename.c_str(), (int)bread);
        return false;
    } 

    // check for invalid nan tile data that could be generated. bug is now
    // fixed, but invalid files may still hang around
    if (isnan(noiseTileData[0])) {
        printf("loadTile: Noise tile '%s' contains nan values.\n", filename.c_str());
        return false;
    }

    return true;
}

// write out an existing noise tile
static void saveTile(float* const noiseTileData, std::string filename)
{
    FILE* file;
    file = fopen(filename.c_str(), "wb");
    int endiantest=1;
    char longsize;

    if (file == NULL) {
        printf("saveTile: Noise tile '%s' could not be saved.\n", filename.c_str());
        return;
    } 

    //Write file header
    fwrite(tilefile_headerstring, strlen(tilefile_headerstring), 1, file);
    fwrite(tilefile_revision, strlen(tilefile_revision), 1, file);
    //Endianness
    if (*((unsigned char*)&endiantest) == 1)
        fwrite("L", 1, 1, file); //Little endian
    else
        fwrite("B",1,1,file); //Big endian
    //32/64bit
    longsize = (char)sizeof(long)+'0';
    fwrite(&longsize, 1, 1, file);
    
    
    fwrite((void*)noiseTileData, sizeof(float), noiseTileSize * noiseTileSize * noiseTileSize, file);
    fclose(file);

    printf("saveTile: Noise tile file '%s' saved.\n", filename.c_str());
}

// create a new noise tile if necessary
static void generateTile_WAVELET(float* const noiseTileData, std::string filename) {
    // if a tile already exists, just use that
    if (loadTile(noiseTileData, filename)) return;

    const int n = noiseTileSize;
    const int n3 = n*n*n;
    std::cout <<"Generating new 3d noise tile size="<<n<<"^3 \n";
    MTRand twister;

    float *temp13 = new float[n3];
    float *temp23 = new float[n3];
    float *noise3 = new float[n3];

    // initialize
    for (int i = 0; i < n3; i++) {
        temp13[i] = temp23[i] = noise3[i] = 0.;
    }

    // Step 1. Fill the tile with random numbers in the range -1 to 1.
    for (int i = 0; i < n3; i++) 
        noise3[i] = twister.randNorm();

    // Steps 2 and 3. Downsample and upsample the tile
    for (int iy = 0; iy < n; iy++) 
        for (int iz = 0; iz < n; iz++) {
            const int i = iy * n + iz*n*n;
            downsampleX(&noise3[i], &temp13[i], n);
            upsampleX  (&temp13[i], &temp23[i], n);
        }
    for (int ix = 0; ix < n; ix++) 
        for (int iz = 0; iz < n; iz++) {
            const int i = ix + iz*n*n;
            downsampleY(&temp23[i], &temp13[i], n);
            upsampleY  (&temp13[i], &temp23[i], n);
        }
    for (int ix = 0; ix < n; ix++) 
        for (int iy = 0; iy < n; iy++) {
            const int i = ix + iy*n;
            downsampleZ(&temp23[i], &temp13[i], n);
            upsampleZ  (&temp13[i], &temp23[i], n);
        }

    // Step 4. Subtract out the coarse-scale contribution
    for (int i = 0; i < n3; i++) 
        noise3[i] -= temp23[i];

    // Avoid even/odd variance difference by adding odd-offset version of noise to itself.
    int offset = n / 2;
    if (offset % 2 == 0) offset++;

    int icnt=0;
    for (int ix = 0; ix < n; ix++)
        for (int iy = 0; iy < n; iy++)
            for (int iz = 0; iz < n; iz++) { 
                temp13[icnt] = noise3[modFast128(ix+offset) + modFast128(iy+offset)*n + modFast128(iz+offset)*n*n];
                icnt++;
            }

    for (int i = 0; i < n3; i++) 
        noise3[i] += temp13[i];

    for (int i = 0; i < n3; i++) 
        noiseTileData[i] = noise3[i];

    saveTile(noise3, filename); 
    delete[] temp13;
    delete[] temp23;
    delete[] noise3;
    std::cout <<"Generating new 3d noise done\n";
}

// x derivative of noise
static inline float WNoiseDx(Vec3 p, float* data) { 
  int c[3], mid[3], n = noiseTileSize;
  float w[3][3], t, result = 0;
  
  mid[0] = (int)ceil(p[0] - 0.5); 
  t = mid[0] - (p[0] - 0.5);
    w[0][0] = -t;
    w[0][2] = (1.f - t);
    w[0][1] = 2.0f * t - 1.0f;
  
  mid[1] = (int)ceil(p[1] - 0.5); 
  t = mid[1] - (p[1] - 0.5);
  w[1][0] = t * t / 2; 
  w[1][2] = (1 - t) * (1 - t) / 2;
  w[1][1] = 1 - w[1][0] - w[1][2];

  mid[2] = (int)ceil(p[2] - 0.5); 
  t = mid[2] - (p[2] - 0.5);
  w[2][0] = t * t / 2; 
  w[2][2] = (1 - t) * (1 - t)/2; 
  w[2][1] = 1 - w[2][0] - w[2][2];
 
  // to optimize, explicitly unroll this loop
  for (int z = -1; z <=1; z++)
    for (int y = -1; y <=1; y++)
      for (int x = -1; x <=1; x++)
      {
        float weight = 1.0f;
        c[0] = modFast128(mid[0] + x);
        weight *= w[0][x+1];
        c[1] = modFast128(mid[1] + y);
        weight *= w[1][y+1];
        c[2] = modFast128(mid[2] + z);
        weight *= w[2][z+1];
        result += weight * data[c[2]*n*n+c[1]*n+c[0]];
      }
 return result;
}

// y derivative of noise
static inline float WNoiseDy(Vec3 p, float* data) { 
  int c[3], mid[3], n=noiseTileSize; 
  float w[3][3], t, result =0;
  
  mid[0] = (int)ceil(p[0] - 0.5); 
  t = mid[0]-(p[0] - 0.5);
  w[0][0] = t * t / 2; 
  w[0][2] = (1 - t) * (1 - t) / 2;
  w[0][1] = 1 - w[0][0] - w[0][2];
  
  mid[1] = (int)ceil(p[1] - 0.5); 
  t = mid[1]-(p[1] - 0.5);
    w[1][0] = -t;
    w[1][2] = (1.f - t);
    w[1][1] = 2.0f * t - 1.0f;

  mid[2] = (int)ceil(p[2] - 0.5); 
  t = mid[2] - (p[2] - 0.5);
  w[2][0] = t * t / 2; 
  w[2][2] = (1 - t) * (1 - t)/2; 
  w[2][1] = 1 - w[2][0] - w[2][2];
  
  // to optimize, explicitly unroll this loop
  for (int z = -1; z <=1; z++)
    for (int y = -1; y <=1; y++)
      for (int x = -1; x <=1; x++)
      {
        float weight = 1.0f;
        c[0] = modFast128(mid[0] + x);
        weight *= w[0][x+1];
        c[1] = modFast128(mid[1] + y);
        weight *= w[1][y+1];
        c[2] = modFast128(mid[2] + z);
        weight *= w[2][z+1];
        result += weight * data[c[2]*n*n+c[1]*n+c[0]];
      }

  return result;
}

// z derivative of noise
static inline float WNoiseDz(Vec3 p, float* data) { 
  int c[3], mid[3], n=noiseTileSize; 
  float w[3][3], t, result =0;

  mid[0] = (int)ceil(p[0] - 0.5); 
  t = mid[0]-(p[0] - 0.5);
  w[0][0] = t * t / 2; 
  w[0][2] = (1 - t) * (1 - t) / 2;
  w[0][1] = 1 - w[0][0] - w[0][2];
  
  mid[1] = (int)ceil(p[1] - 0.5); 
  t = mid[1]-(p[1] - 0.5);
  w[1][0] = t * t / 2; 
  w[1][2] = (1 - t) * (1 - t) / 2;
  w[1][1] = 1 - w[1][0] - w[1][2];

  mid[2] = (int)ceil(p[2] - 0.5); 
  t = mid[2] - (p[2] - 0.5);
    w[2][0] = -t;
    w[2][2] = (1.f - t);
    w[2][1] = 2.0f * t - 1.0f;

  // to optimize, explicitly unroll this loop
  for (int z = -1; z <=1; z++)
    for (int y = -1; y <=1; y++)
      for (int x = -1; x <=1; x++)
      {
        float weight = 1.0f;
        c[0] = modFast128(mid[0] + x);
        weight *= w[0][x+1];
        c[1] = modFast128(mid[1] + y);
        weight *= w[1][y+1];
        c[2] = modFast128(mid[2] + z);
        weight *= w[2][z+1];
        result += weight * data[c[2]*n*n+c[1]*n+c[0]];
      }
  return result;
}

#endif