forked from bartvdbraak/blender
520 lines
16 KiB
C++
520 lines
16 KiB
C++
/** \file smoke/intern/WAVELET_NOISE.h
|
|
* \ingroup smoke
|
|
*/
|
|
//////////////////////////////////////////////////////////////////////
|
|
// 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>
|
|
|
|
#include <string.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.5f);
|
|
t = mid[0] - (p[0] - 0.5f);
|
|
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.5f);
|
|
t = mid[1] - (p[1] - 0.5f);
|
|
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.5f);
|
|
t = mid[2] - (p[2] - 0.5f);
|
|
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.5f);
|
|
t = mid[0]-(p[0] - 0.5f);
|
|
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.5f);
|
|
t = mid[1]-(p[1] - 0.5f);
|
|
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.5f);
|
|
t = mid[2] - (p[2] - 0.5f);
|
|
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.5f);
|
|
t = mid[0]-(p[0] - 0.5f);
|
|
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.5f);
|
|
t = mid[1]-(p[1] - 0.5f);
|
|
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.5f);
|
|
t = mid[2] - (p[2] - 0.5f);
|
|
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
|
|
|