Blender V2.61 - r43446

WAVELET_NOISE.h

Go to the documentation of this file.
00001 
00004 
00005 // This file is part of Wavelet Turbulence.
00006 // 
00007 // Wavelet Turbulence is free software: you can redistribute it and/or modify
00008 // it under the terms of the GNU General Public License as published by
00009 // the Free Software Foundation, either version 3 of the License, or
00010 // (at your option) any later version.
00011 // 
00012 // Wavelet Turbulence is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // GNU General Public License for more details.
00016 // 
00017 // You should have received a copy of the GNU General Public License
00018 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
00019 // 
00020 // Copyright 2008 Theodore Kim and Nils Thuerey
00021 // 
00023 // Wavelet noise functions
00024 //
00025 // This code is based on the C code provided in the appendices of:
00026 //
00027 // @article{1073264,
00028 //  author = {Robert L. Cook and Tony DeRose},
00029 //  title = {Wavelet noise},
00030 //  journal = {ACM Trans. Graph.},
00031 //  volume = {24},
00032 //  number = {3},
00033 //  year = {2005},
00034 //  issn = {0730-0301},
00035 //  pages = {803--811},
00036 //  doi = {http://doi.acm.org/10.1145/1073204.1073264},
00037 //  publisher = {ACM},
00038 //  address = {New York, NY, USA},
00039 //  }
00040 //
00042 
00043 #ifndef WAVELET_NOISE_H
00044 #define WAVELET_NOISE_H
00045 
00046 #include <MERSENNETWISTER.h>
00047 
00048 #ifdef WIN32
00049 #include <float.h>
00050 #define isnan _isnan
00051 #endif
00052 
00053 // Tile file header, update revision upon any change done to the noise generator
00054 static const char tilefile_headerstring[] = "Noise Tile File rev. ";
00055 static const char tilefile_revision[] =  "001";
00056 
00057 #define NOISE_TILE_SIZE 128
00058 static const int noiseTileSize = NOISE_TILE_SIZE;
00059 
00060 // warning - noiseTileSize has to be 128^3!
00061 #define modFast128(x) ((x) & 127)
00062 #define modFast64(x)  ((x) & 63)
00063 #define DOWNCOEFFS 0.000334f,-0.001528f, 0.000410f, 0.003545f,-0.000938f,-0.008233f, 0.002172f, 0.019120f, \
00064                   -0.005040f,-0.044412f, 0.011655f, 0.103311f,-0.025936f,-0.243780f, 0.033979f, 0.655340f, \
00065                    0.655340f, 0.033979f,-0.243780f,-0.025936f, 0.103311f, 0.011655f,-0.044412f,-0.005040f, \
00066                    0.019120f, 0.002172f,-0.008233f,-0.000938f, 0.003546f, 0.000410f,-0.001528f, 0.000334f
00067 
00069 // Wavelet downsampling -- periodic boundary conditions
00071 static void downsampleX(float *from, float *to, int n){
00072   // if these values are not local incorrect results are generated
00073   float downCoeffs[32] = { DOWNCOEFFS };
00074     const float *a = &downCoeffs[16];
00075     for (int i = 0; i < n / 2; i++) {
00076         to[i] = 0;
00077         for (int k = 2 * i - 16; k < 2 * i + 16; k++)
00078             to[i] += a[k - 2 * i] * from[modFast128(k)];
00079     }
00080 }
00081 static void downsampleY(float *from, float *to, int n){
00082   // if these values are not local incorrect results are generated
00083   float downCoeffs[32] = { DOWNCOEFFS };
00084     const float *a = &downCoeffs[16];
00085     for (int i = 0; i < n / 2; i++) {
00086         to[i * n] = 0;
00087         for (int k = 2 * i - 16; k < 2 * i + 16; k++)
00088             to[i * n] += a[k - 2 * i] * from[modFast128(k) * n];
00089     }
00090 }
00091 static void downsampleZ(float *from, float *to, int n){
00092   // if these values are not local incorrect results are generated
00093   float downCoeffs[32] = { DOWNCOEFFS };
00094     const float *a = &downCoeffs[16];
00095     for (int i = 0; i < n / 2; i++) {
00096         to[i * n * n] = 0;
00097         for (int k = 2 * i - 16; k < 2 * i + 16; k++)
00098             to[i * n * n] += a[k - 2 * i] * from[modFast128(k) * n * n];
00099     }
00100 }
00101 
00103 // Wavelet downsampling -- Neumann boundary conditions
00105 static void downsampleNeumann(const float *from, float *to, int n, int stride)
00106 {
00107   // if these values are not local incorrect results are generated
00108   float downCoeffs[32] = { DOWNCOEFFS };
00109   const float *const aCoCenter= &downCoeffs[16];
00110     for (int i = 0; i < n / 2; i++) {
00111         to[i * stride] = 0;
00112         for (int k = 2 * i - 16; k < 2 * i + 16; k++) { 
00113             // handle boundary
00114             float fromval; 
00115             if (k < 0) {
00116                 fromval = from[0];
00117             } else if(k > n - 1) {
00118                 fromval = from[(n - 1) * stride];
00119             } else {
00120                 fromval = from[k * stride]; 
00121             } 
00122             to[i * stride] += aCoCenter[k - 2 * i] * fromval; 
00123         }
00124     }
00125 }
00126 static void downsampleXNeumann(float* to, const float* from, int sx,int sy, int sz) {
00127     for (int iy = 0; iy < sy; iy++) 
00128         for (int iz = 0; iz < sz; iz++) {
00129             const int i = iy * sx + iz*sx*sy;
00130             downsampleNeumann(&from[i], &to[i], sx, 1);
00131         }
00132 }
00133 static void downsampleYNeumann(float* to, const float* from, int sx,int sy, int sz) {
00134     for (int ix = 0; ix < sx; ix++) 
00135         for (int iz = 0; iz < sz; iz++) {
00136             const int i = ix + iz*sx*sy;
00137             downsampleNeumann(&from[i], &to[i], sy, sx);
00138     }
00139 }
00140 static void downsampleZNeumann(float* to, const float* from, int sx,int sy, int sz) {
00141     for (int ix = 0; ix < sx; ix++) 
00142         for (int iy = 0; iy < sy; iy++) {
00143             const int i = ix + iy*sx;
00144             downsampleNeumann(&from[i], &to[i], sz, sx*sy);
00145     }
00146 }
00147 
00149 // Wavelet upsampling - periodic boundary conditions
00151 static float _upCoeffs[4] = {0.25f, 0.75f, 0.75f, 0.25f};
00152 static void upsampleX(float *from, float *to, int n) {
00153     const float *p = &_upCoeffs[2];
00154 
00155     for (int i = 0; i < n; i++) {
00156         to[i] = 0;
00157         for (int k = i / 2; k <= i / 2 + 1; k++)
00158             to[i] += p[i - 2 * k] * from[modFast64(k)];
00159     }
00160 }
00161 static void upsampleY(float *from, float *to, int n) {
00162     const float *p = &_upCoeffs[2];
00163 
00164     for (int i = 0; i < n; i++) {
00165         to[i * n] = 0;
00166         for (int k = i / 2; k <= i / 2 + 1; k++)
00167             to[i * n] += p[i - 2 * k] * from[modFast64(k) * n];
00168     }
00169 }
00170 static void upsampleZ(float *from, float *to, int n) {
00171     const float *p = &_upCoeffs[2];
00172 
00173     for (int i = 0; i < n; i++) {
00174         to[i * n * n] = 0;
00175         for (int k = i / 2; k <= i / 2 + 1; k++)
00176             to[i * n * n] += p[i - 2 * k] * from[modFast64(k) * n * n];
00177     }
00178 }
00179 
00181 // Wavelet upsampling - Neumann boundary conditions
00183 static void upsampleNeumann(const float *from, float *to, int n, int stride) {
00184   static const float *const pCoCenter = &_upCoeffs[2];
00185     for (int i = 0; i < n; i++) {
00186         to[i * stride] = 0;
00187         for (int k = i / 2; k <= i / 2 + 1; k++) {
00188             float fromval;
00189             if(k>n/2) {
00190                 fromval = from[(n/2) * stride];
00191             } else {
00192                 fromval = from[k * stride]; 
00193             }  
00194             to[i * stride] += pCoCenter[i - 2 * k] * fromval; 
00195         }
00196     }
00197 }
00198 static void upsampleXNeumann(float* to, const float* from, int sx, int sy, int sz) {
00199     for (int iy = 0; iy < sy; iy++) 
00200         for (int iz = 0; iz < sz; iz++) {
00201             const int i = iy * sx + iz*sx*sy;
00202             upsampleNeumann(&from[i], &to[i], sx, 1);
00203         }
00204 }
00205 static void upsampleYNeumann(float* to, const float* from, int sx, int sy, int sz) {
00206     for (int ix = 0; ix < sx; ix++) 
00207         for (int iz = 0; iz < sz; iz++) {
00208             const int i = ix + iz*sx*sy;
00209             upsampleNeumann(&from[i], &to[i], sy, sx);
00210         }
00211 }
00212 static void upsampleZNeumann(float* to, const float* from, int sx, int sy, int sz) {
00213     for (int ix = 0; ix < sx; ix++) 
00214         for (int iy = 0; iy < sy; iy++) {
00215             const int i = ix + iy*sx;
00216             upsampleNeumann(&from[i], &to[i], sz, sx*sy);
00217         }
00218 }
00219 
00220 
00222 // load in an existing noise tile
00224 static bool loadTile(float* const noiseTileData, std::string filename)
00225 {
00226     FILE* file;
00227     char headerbuffer[64];
00228     size_t headerlen;
00229     size_t bread;
00230     int endiantest = 1;
00231     char endianness;
00232     
00233     file = fopen(filename.c_str(), "rb");
00234 
00235     if (file == NULL) {
00236         printf("loadTile: No noise tile '%s' found.\n", filename.c_str());
00237         return false;
00238     }
00239 
00240     //Check header
00241     headerlen = strlen(tilefile_headerstring) + strlen(tilefile_revision) + 2;
00242     bread = fread((void*)headerbuffer, 1, headerlen, file);
00243     if (*((unsigned char*)&endiantest) == 1)
00244         endianness = 'L';
00245     else
00246         endianness = 'B';
00247     if ((bread != headerlen)
00248         || (strncmp(headerbuffer, tilefile_headerstring, strlen(tilefile_headerstring)))
00249         || (strncmp(headerbuffer+ strlen(tilefile_headerstring), tilefile_revision, strlen(tilefile_revision)))
00250         || (headerbuffer[headerlen-2] != endianness)
00251         || (headerbuffer[headerlen-1] != (char)((char)sizeof(long)+'0')))
00252     {
00253         printf("loadTile : Noise tile '%s' was generated on an incompatible platform.\n",filename.c_str());
00254         fclose(file);
00255         return false;
00256     }
00257     
00258     // dimensions
00259     size_t gridSize = noiseTileSize * noiseTileSize * noiseTileSize;
00260 
00261     // noiseTileData memory is managed by caller
00262     bread = fread((void*)noiseTileData, sizeof(float), gridSize, file);
00263     fclose(file);
00264     printf("Noise tile file '%s' loaded.\n", filename.c_str());
00265 
00266     if (bread != gridSize) {
00267         printf("loadTile: Noise tile '%s' is wrong size %d.\n", filename.c_str(), (int)bread);
00268         return false;
00269     } 
00270 
00271     // check for invalid nan tile data that could be generated. bug is now
00272     // fixed, but invalid files may still hang around
00273     if (isnan(noiseTileData[0])) {
00274         printf("loadTile: Noise tile '%s' contains nan values.\n", filename.c_str());
00275         return false;
00276     }
00277 
00278     return true;
00279 }
00280 
00282 // write out an existing noise tile
00284 static void saveTile(float* const noiseTileData, std::string filename)
00285 {
00286     FILE* file;
00287     file = fopen(filename.c_str(), "wb");
00288     int endiantest=1;
00289     char longsize;
00290 
00291     if (file == NULL) {
00292         printf("saveTile: Noise tile '%s' could not be saved.\n", filename.c_str());
00293         return;
00294     } 
00295 
00296     //Write file header
00297     fwrite(tilefile_headerstring, strlen(tilefile_headerstring), 1, file);
00298     fwrite(tilefile_revision, strlen(tilefile_revision), 1, file);
00299     //Endianness
00300     if (*((unsigned char*)&endiantest) == 1)
00301         fwrite("L", 1, 1, file); //Little endian
00302     else
00303         fwrite("B",1,1,file); //Big endian
00304     //32/64bit
00305     longsize = (char)sizeof(long)+'0';
00306     fwrite(&longsize, 1, 1, file);
00307     
00308     
00309     fwrite((void*)noiseTileData, sizeof(float), noiseTileSize * noiseTileSize * noiseTileSize, file);
00310     fclose(file);
00311 
00312     printf("saveTile: Noise tile file '%s' saved.\n", filename.c_str());
00313 }
00314 
00316 // create a new noise tile if necessary
00318 static void generateTile_WAVELET(float* const noiseTileData, std::string filename) {
00319     // if a tile already exists, just use that
00320     if (loadTile(noiseTileData, filename)) return;
00321 
00322     const int n = noiseTileSize;
00323     const int n3 = n*n*n;
00324     std::cout <<"Generating new 3d noise tile size="<<n<<"^3 \n";
00325     MTRand twister;
00326 
00327     float *temp13 = new float[n3];
00328     float *temp23 = new float[n3];
00329     float *noise3 = new float[n3];
00330 
00331     // initialize
00332     for (int i = 0; i < n3; i++) {
00333         temp13[i] = temp23[i] = noise3[i] = 0.;
00334     }
00335 
00336     // Step 1. Fill the tile with random numbers in the range -1 to 1.
00337     for (int i = 0; i < n3; i++) 
00338         noise3[i] = twister.randNorm();
00339 
00340     // Steps 2 and 3. Downsample and upsample the tile
00341     for (int iy = 0; iy < n; iy++) 
00342         for (int iz = 0; iz < n; iz++) {
00343             const int i = iy * n + iz*n*n;
00344             downsampleX(&noise3[i], &temp13[i], n);
00345             upsampleX  (&temp13[i], &temp23[i], n);
00346         }
00347     for (int ix = 0; ix < n; ix++) 
00348         for (int iz = 0; iz < n; iz++) {
00349             const int i = ix + iz*n*n;
00350             downsampleY(&temp23[i], &temp13[i], n);
00351             upsampleY  (&temp13[i], &temp23[i], n);
00352         }
00353     for (int ix = 0; ix < n; ix++) 
00354         for (int iy = 0; iy < n; iy++) {
00355             const int i = ix + iy*n;
00356             downsampleZ(&temp23[i], &temp13[i], n);
00357             upsampleZ  (&temp13[i], &temp23[i], n);
00358         }
00359 
00360     // Step 4. Subtract out the coarse-scale contribution
00361     for (int i = 0; i < n3; i++) 
00362         noise3[i] -= temp23[i];
00363 
00364     // Avoid even/odd variance difference by adding odd-offset version of noise to itself.
00365     int offset = n / 2;
00366     if (offset % 2 == 0) offset++;
00367 
00368     int icnt=0;
00369     for (int ix = 0; ix < n; ix++)
00370         for (int iy = 0; iy < n; iy++)
00371             for (int iz = 0; iz < n; iz++) { 
00372                 temp13[icnt] = noise3[modFast128(ix+offset) + modFast128(iy+offset)*n + modFast128(iz+offset)*n*n];
00373                 icnt++;
00374             }
00375 
00376     for (int i = 0; i < n3; i++) 
00377         noise3[i] += temp13[i];
00378 
00379     for (int i = 0; i < n3; i++) 
00380         noiseTileData[i] = noise3[i];
00381 
00382     saveTile(noise3, filename); 
00383     delete[] temp13;
00384     delete[] temp23;
00385     delete[] noise3;
00386     std::cout <<"Generating new 3d noise done\n";
00387 }
00388 
00390 // x derivative of noise
00392 static inline float WNoiseDx(Vec3 p, float* data) { 
00393   int c[3], mid[3], n = noiseTileSize;
00394   float w[3][3], t, result = 0;
00395   
00396   mid[0] = (int)ceil(p[0] - 0.5); 
00397   t = mid[0] - (p[0] - 0.5);
00398     w[0][0] = -t;
00399     w[0][2] = (1.f - t);
00400     w[0][1] = 2.0f * t - 1.0f;
00401   
00402   mid[1] = (int)ceil(p[1] - 0.5); 
00403   t = mid[1] - (p[1] - 0.5);
00404   w[1][0] = t * t / 2; 
00405   w[1][2] = (1 - t) * (1 - t) / 2;
00406   w[1][1] = 1 - w[1][0] - w[1][2];
00407 
00408   mid[2] = (int)ceil(p[2] - 0.5); 
00409   t = mid[2] - (p[2] - 0.5);
00410   w[2][0] = t * t / 2; 
00411   w[2][2] = (1 - t) * (1 - t)/2; 
00412   w[2][1] = 1 - w[2][0] - w[2][2];
00413  
00414   // to optimize, explicitly unroll this loop
00415   for (int z = -1; z <=1; z++)
00416     for (int y = -1; y <=1; y++)
00417       for (int x = -1; x <=1; x++)
00418       {
00419         float weight = 1.0f;
00420         c[0] = modFast128(mid[0] + x);
00421         weight *= w[0][x+1];
00422         c[1] = modFast128(mid[1] + y);
00423         weight *= w[1][y+1];
00424         c[2] = modFast128(mid[2] + z);
00425         weight *= w[2][z+1];
00426         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
00427       }
00428  return result;
00429 }
00430 
00432 // y derivative of noise
00434 static inline float WNoiseDy(Vec3 p, float* data) { 
00435   int c[3], mid[3], n=noiseTileSize; 
00436   float w[3][3], t, result =0;
00437   
00438   mid[0] = (int)ceil(p[0] - 0.5); 
00439   t = mid[0]-(p[0] - 0.5);
00440   w[0][0] = t * t / 2; 
00441   w[0][2] = (1 - t) * (1 - t) / 2;
00442   w[0][1] = 1 - w[0][0] - w[0][2];
00443   
00444   mid[1] = (int)ceil(p[1] - 0.5); 
00445   t = mid[1]-(p[1] - 0.5);
00446     w[1][0] = -t;
00447     w[1][2] = (1.f - t);
00448     w[1][1] = 2.0f * t - 1.0f;
00449 
00450   mid[2] = (int)ceil(p[2] - 0.5); 
00451   t = mid[2] - (p[2] - 0.5);
00452   w[2][0] = t * t / 2; 
00453   w[2][2] = (1 - t) * (1 - t)/2; 
00454   w[2][1] = 1 - w[2][0] - w[2][2];
00455   
00456   // to optimize, explicitly unroll this loop
00457   for (int z = -1; z <=1; z++)
00458     for (int y = -1; y <=1; y++)
00459       for (int x = -1; x <=1; x++)
00460       {
00461         float weight = 1.0f;
00462         c[0] = modFast128(mid[0] + x);
00463         weight *= w[0][x+1];
00464         c[1] = modFast128(mid[1] + y);
00465         weight *= w[1][y+1];
00466         c[2] = modFast128(mid[2] + z);
00467         weight *= w[2][z+1];
00468         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
00469       }
00470 
00471   return result;
00472 }
00473 
00475 // z derivative of noise
00477 static inline float WNoiseDz(Vec3 p, float* data) { 
00478   int c[3], mid[3], n=noiseTileSize; 
00479   float w[3][3], t, result =0;
00480 
00481   mid[0] = (int)ceil(p[0] - 0.5); 
00482   t = mid[0]-(p[0] - 0.5);
00483   w[0][0] = t * t / 2; 
00484   w[0][2] = (1 - t) * (1 - t) / 2;
00485   w[0][1] = 1 - w[0][0] - w[0][2];
00486   
00487   mid[1] = (int)ceil(p[1] - 0.5); 
00488   t = mid[1]-(p[1] - 0.5);
00489   w[1][0] = t * t / 2; 
00490   w[1][2] = (1 - t) * (1 - t) / 2;
00491   w[1][1] = 1 - w[1][0] - w[1][2];
00492 
00493   mid[2] = (int)ceil(p[2] - 0.5); 
00494   t = mid[2] - (p[2] - 0.5);
00495     w[2][0] = -t;
00496     w[2][2] = (1.f - t);
00497     w[2][1] = 2.0f * t - 1.0f;
00498 
00499   // to optimize, explicitly unroll this loop
00500   for (int z = -1; z <=1; z++)
00501     for (int y = -1; y <=1; y++)
00502       for (int x = -1; x <=1; x++)
00503       {
00504         float weight = 1.0f;
00505         c[0] = modFast128(mid[0] + x);
00506         weight *= w[0][x+1];
00507         c[1] = modFast128(mid[1] + y);
00508         weight *= w[1][y+1];
00509         c[2] = modFast128(mid[2] + z);
00510         weight *= w[2][z+1];
00511         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
00512       }
00513   return result;
00514 }
00515 
00516 #endif
00517