+void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
+{
+ int x, y;
+ double nx, ny, nz;
+ double v;
+#ifndef C99
+ double save;
+#endif
+
+ fftw_complex *imgspace1 = fftw_malloc(w*h * sizeof(fftw_complex));
+ fftw_complex *imgspace2 = fftw_malloc(w*h * sizeof(fftw_complex));
+ fftw_complex *freqspace1 = fftw_malloc(w*h * sizeof(fftw_complex));
+ fftw_complex *freqspace2 = fftw_malloc(w*h * sizeof(fftw_complex));
+ fftw_plan i12f1 = fftw_plan_dft_2d(h, w, imgspace1, freqspace1, FFTW_FORWARD, FFTW_ESTIMATE);
+ fftw_plan f12i1 = fftw_plan_dft_2d(h, w, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
+ fftw_plan f22i2 = fftw_plan_dft_2d(h, w, freqspace2, imgspace2, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+ for(y = 0; y < h; ++y)
+ for(x = 0; x < w; ++x)
+ {
+ switch(src_chan)
+ {
+ case 0:
+ case 1:
+ case 2:
+ case 3:
+ v = map[(w*y+x)*4+src_chan];
+ break;
+ case 4:
+ v = (map[(w*y+x)*4+0] + map[(w*y+x)*4+1] + map[(w*y+x)*4+2]) / 3;
+ break;
+ default:
+ case 5:
+ v = (map[(w*y+x)*4+0]*0.114 + map[(w*y+x)*4+1]*0.587 + map[(w*y+x)*4+2]*0.299);
+ break;
+ }
+#ifdef C99
+ imgspace1[(w*y+x)] = (v - 128.0) / 127.0;
+#else
+ imgspace1[(w*y+x)][0] = (v - 128.0) / 127.0;
+ imgspace1[(w*y+x)][1] = 0;
+#endif
+ map[(w*y+x)*4+3] = floor(v + 0.5);
+ }
+
+ /* see http://www.gamedev.net/community/forums/topic.asp?topic_id=561430 */
+
+ fftw_execute(i12f1);
+
+ for(y = 0; y < h; ++y)
+ for(x = 0; x < w; ++x)
+ {
+ int fx = x;
+ int fy = y;
+ if(fx > w/2)
+ fx -= w;
+ if(fy > h/2)
+ fy -= h;
+#ifdef C99
+ /* a lowpass to prevent the worst */
+ freqspace1[(w*y+x)] *= 1 - pow(abs(fx) / (double)(w/2), 1);
+ freqspace1[(w*y+x)] *= 1 - pow(abs(fy) / (double)(h/2), 1);
+
+ freqspace2[(w*y+x)] = TWO_PI*_Complex_I * fy * freqspace1[(w*y+x)]; /* y derivative */
+ freqspace1[(w*y+x)] = TWO_PI*_Complex_I * fx * freqspace1[(w*y+x)]; /* x derivative */
+#else
+ /* a lowpass to prevent the worst */
+ freqspace1[(w*y+x)][0] *= 1 - pow(abs(fx) / (double)(w/2), 1);
+ freqspace1[(w*y+x)][1] *= 1 - pow(abs(fx) / (double)(w/2), 1);
+ freqspace1[(w*y+x)][0] *= 1 - pow(abs(fy) / (double)(h/2), 1);
+ freqspace1[(w*y+x)][1] *= 1 - pow(abs(fy) / (double)(h/2), 1);
+
+ freqspace2[(w*y+x)][0] = -TWO_PI * fy * freqspace1[(w*y+x)][1]; /* y derivative */
+ freqspace2[(w*y+x)][1] = TWO_PI * fy * freqspace1[(w*y+x)][0];
+ save = freqspace1[(w*y+x)][0];
+ freqspace1[(w*y+x)][0] = -TWO_PI * fx * freqspace1[(w*y+x)][1]; /* x derivative */
+ freqspace1[(w*y+x)][1] = TWO_PI * fx * save;
+#endif
+ }
+
+ fftw_execute(f12i1);
+ fftw_execute(f22i2);
+
+ scale /= (w*h);
+
+ for(y = 0; y < h; ++y)
+ for(x = 0; x < w; ++x)
+ {
+#ifdef C99
+ nx = creal(imgspace1[(w*y+x)]);
+ ny = creal(imgspace2[(w*y+x)]);
+#else
+ nx = imgspace1[(w*y+x)][0];
+ ny = imgspace2[(w*y+x)][0];
+#endif
+ nx /= w;
+ ny /= h;
+ nz = -1 / scale;
+ v = -sqrt(nx*nx + ny*ny + nz*nz);
+ nx /= v;
+ ny /= v;
+ nz /= v;
+ ny = -ny; /* DP inverted normals */
+ map[(w*y+x)*4+2] = floor(128 + 127.5 * nx);
+ map[(w*y+x)*4+1] = floor(128 + 127.5 * ny);
+ map[(w*y+x)*4+0] = floor(128 + 127.5 * nz);
+ }
+
+ fftw_destroy_plan(i12f1);
+ fftw_destroy_plan(f12i1);
+ fftw_destroy_plan(f22i2);
+
+ fftw_free(freqspace2);
+ fftw_free(freqspace1);
+ fftw_free(imgspace2);
+ fftw_free(imgspace1);
+}
+
+void hmap_to_nmap_local(unsigned char *map, int w, int h, int src_chan, double scale)
+{
+ int x, y;
+ double nx, ny, nz;
+ double v;
+ int i, j;
+ double *img_reduced = malloc(w*h * sizeof(double));
+ static const double filter[3][3] = { /* filter to derive one component */
+ { -1, 0, 1 },
+ { -2, 0, 2 },
+ { -1, 0, 1 }
+ };
+ static const double filter_mult = 0.125;
+
+ for(y = 0; y < h; ++y)
+ for(x = 0; x < w; ++x)
+ {
+ switch(src_chan)
+ {
+ case 0:
+ case 1:
+ case 2:
+ case 3:
+ v = map[(w*y+x)*4+src_chan];
+ break;
+ case 4:
+ v = (map[(w*y+x)*4+0] + map[(w*y+x)*4+1] + map[(w*y+x)*4+2]) / 3;
+ break;
+ default:
+ case 5:
+ v = (map[(w*y+x)*4+0]*0.114 + map[(w*y+x)*4+1]*0.587 + map[(w*y+x)*4+2]*0.299);
+ break;
+ }
+ img_reduced[(w*y+x)] = (v - 128.0) / 127.0;
+ map[(w*y+x)*4+3] = floor(v + 0.5);
+ }
+
+ for(y = 0; y < h; ++y)
+ for(x = 0; x < w; ++x)
+ {
+ nz = -1 / (scale * filter_mult);
+ nx = ny = 0;
+
+ for(i = -(int)(sizeof(filter) / sizeof(*filter)) / 2; i <= (int)(sizeof(filter) / sizeof(*filter)) / 2; ++i)
+ for(j = -(int)(sizeof(*filter) / sizeof(**filter)) / 2; j <= (int)(sizeof(*filter) / sizeof(**filter)) / 2; ++j)
+ {
+ nx += img_reduced[w*((y+i+h)%h)+(x+j+w)%w] * filter[i+(sizeof(filter) / sizeof(*filter)) / 2][j+(sizeof(*filter) / sizeof(**filter)) / 2];
+ ny += img_reduced[w*((y+j+h)%h)+(x+i+w)%w] * filter[i+(sizeof(filter) / sizeof(*filter)) / 2][j+(sizeof(*filter) / sizeof(**filter)) / 2];
+ }
+
+ v = -sqrt(nx*nx + ny*ny + nz*nz);
+ nx /= v;
+ ny /= v;
+ nz /= v;
+ ny = -ny; /* DP inverted normals */
+ map[(w*y+x)*4+2] = floor(128 + 127.5 * nx);
+ map[(w*y+x)*4+1] = floor(128 + 127.5 * ny);
+ map[(w*y+x)*4+0] = floor(128 + 127.5 * nz);
+ }
+
+ free(img_reduced);
+}
+