--- FFT_NOISE.h 2012-05-14 16:03:56.000000000 -0500 +++ FFT_NOISE.h 2012-05-15 19:41:14.093398574 -0500 @@ -41,40 +41,49 @@ ///////////////////////////////////////////////////////////////////////// static void shift3D(float*& field, int xRes, int yRes, int zRes) { - int xHalf = xRes / 2; - int yHalf = yRes / 2; - int zHalf = zRes / 2; - // int slabSize = xRes * yRes; - for (int z = 0; z < zHalf; z++) - for (int y = 0; y < yHalf; y++) - for (int x = 0; x < xHalf; x++) - { - int index = x + y * xRes + z * xRes * yRes; - float temp; - int xSwap = xHalf; - int ySwap = yHalf * xRes; - int zSwap = zHalf * xRes * yRes; - - // [0,0,0] to [1,1,1] - temp = field[index]; - field[index] = field[index + xSwap + ySwap + zSwap]; - field[index + xSwap + ySwap + zSwap] = temp; - - // [1,0,0] to [0,1,1] - temp = field[index + xSwap]; - field[index + xSwap] = field[index + ySwap + zSwap]; - field[index + ySwap + zSwap] = temp; - - // [0,1,0] to [1,0,1] - temp = field[index + ySwap]; - field[index + ySwap] = field[index + xSwap + zSwap]; - field[index + xSwap + zSwap] = temp; - - // [0,0,1] to [1,1,0] - temp = field[index + zSwap]; - field[index + zSwap] = field[index + xSwap + ySwap]; - field[index + xSwap + ySwap] = temp; - } + int x = 0, y = 0, z = 0; /* loop variables */ + + int index = 0; /* field index variable */ + + /* Divisions by 2 replaced with bitwise right shifts */ + int xHalf = (int)((unsigned int)xRes >> 1); + int yHalf = (int)((unsigned int)yRes >> 1); + int zHalf = (int)((unsigned int)zRes >> 1); + + /* Removal of redundant inline calculations (see notes) */ + int xSwap = xHalf; + int ySwap = yHalf * xRes; + int zSwap = zHalf * xRes * yRes; + + float temp; + + for(z = 0; z < zHalf; ++z) { + for(y = 0; y < yHalf; ++y) { + for(x = 0; x < xHalf; ++x) { + index = x + y * xRes + z * xRes * yRes; + + // [0,0,0] to [1,1,1] + temp = field[index]; + field[index] = field[index + xSwap + ySwap + zSwap]; + field[index + xSwap + ySwap + zSwap] = temp; + + // [1,0,0] to [0,1,1] + temp = field[index + xSwap]; + field[index + xSwap] = field[index + ySwap + zSwap]; + field[index + ySwap + zSwap] = temp; + + // [0,1,0] to [1,0,1] + temp = field[index + ySwap]; + field[index + ySwap] = field[index + xSwap + zSwap]; + field[index + xSwap + zSwap] = temp; + + // [0,0,1] to [1,1,0] + temp = field[index + zSwap]; + field[index + zSwap] = field[index + xSwap + ySwap]; + field[index + xSwap + ySwap] = temp; + } + } + } } static void generatTile_FFT(float* const noiseTileData, std::string filename) @@ -85,59 +94,78 @@ int xRes = res; int yRes = res; int zRes = res; + + int x = 0, y = 0, z = 0; + int index = 0; + int totalCells = xRes * yRes * zRes; - + float inv_totalCells = 1.0f / (float)totalCells; + + int xHalf = (int)((unsigned int)xRes >> 1); + int yHalf = (int)((unsigned int)yRes >> 1); + int zHalf = (int)((unsigned int)zRes >> 1); + + /* Time saving pre-calculations */ + float inv_log2 = 1.0f / log(2.0f); + float inv_pi = 1.0f / M_PI; + float pi_div2 = M_PI * 0.5f; + float pi_div4 = M_PI * 0.25f; + // create and shift the filter float* filter = new float[totalCells]; - for (int z = 0; z < zRes; z++) - for (int y = 0; y < yRes; y++) - for (int x = 0; x < xRes; x++) - { - int index = x + y * xRes + z * xRes * yRes; - float diff[] = {abs(x - xRes/2), - abs(y - yRes/2), - abs(z - zRes/2)}; - float radius = sqrtf(diff[0] * diff[0] + - diff[1] * diff[1] + - diff[2] * diff[2]) / (xRes / 2); + for (z = 0; z < zRes; ++z) { + for (y = 0; y < yRes; ++y) { + for (x = 0; x < xRes; ++x) { + + index = x + y * xRes + z * xRes * yRes; + + float diff[] = { + abs(x - xHalf), + abs(y - yHalf), + abs(z - zHalf) + }; + + float radius = sqrtf( diff[0] * diff[0] + \ + diff[1] * diff[1] + \ + diff[2] * diff[2]) / xHalf; + radius *= M_PI; - float H = cos((M_PI / 2.0f) * log(4.0f * radius / M_PI) / log(2.0f)); + + float H = cos(pi_div2 * log(4.0f * radius * inv_pi) * inv_log2); H = H * H; + float filtered = H; // clamp everything outside the wanted band - if (radius >= M_PI / 2.0f) + if (radius >= pi_div2) filtered = 0.0f; // make sure to capture all low frequencies - if (radius <= M_PI / 4.0f) + if (radius <= pi_div4) filtered = 1.0f; filter[index] = filtered; } + } + } + shift3D(filter, xRes, yRes, zRes); // create the noise float* noise = new float[totalCells]; - int index = 0; MTRand twister; - for (int z = 0; z < zRes; z++) - for (int y = 0; y < yRes; y++) - for (int x = 0; x < xRes; x++, index++) + for(index = 0; index < totalCells; ++index) noise[index] = twister.randNorm(); // create padded field fftw_complex* forward = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * totalCells); // init padded field - index = 0; - for (int z = 0; z < zRes; z++) - for (int y = 0; y < yRes; y++) - for (int x = 0; x < xRes; x++, index++) - { + for(index = 0; index < totalCells; ++index) { forward[index][0] = noise[index]; forward[index][1] = 0.0f; - } + } // forward FFT fftw_complex* backward = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * totalCells); @@ -146,14 +174,10 @@ fftw_destroy_plan(forwardPlan); // apply filter - index = 0; - for (int z = 0; z < zRes; z++) - for (int y = 0; y < yRes; y++) - for (int x = 0; x < xRes; x++, index++) - { - backward[index][0] *= filter[index]; - backward[index][1] *= filter[index]; - } + for(index = 0; index < totalCells; ++index) { + backward[index][0] *= filter[index]; + backward[index][1] *= filter[index]; + } // backward FFT fftw_plan backwardPlan = fftw_plan_dft_3d(xRes, yRes, zRes, backward, forward, FFTW_BACKWARD, FFTW_ESTIMATE); @@ -161,11 +185,8 @@ fftw_destroy_plan(backwardPlan); // subtract out the low frequency components - index = 0; - for (int z = 0; z < zRes; z++) - for (int y = 0; y < yRes; y++) - for (int x = 0; x < xRes; x++, index++) - noise[index] -= forward[index][0] / totalCells; + for(index = 0; index < totalCells; ++index) + noise[index] -= forward[index][0] * inv_totalCells; // save out the noise tile saveTile(noise, filename);