+ free(medianbuf);
+}
+
+void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
+{
+ int x, y;
+ double fx, fy;
+ double nx, ny, nz;
+ double v;
+#ifndef C99
+ double save;
+#endif
+
+ fftw_complex *imgspace1 = (fftw_complex *) fftw_malloc(w*h * sizeof(fftw_complex));
+ fftw_complex *imgspace2 = (fftw_complex *) fftw_malloc(w*h * sizeof(fftw_complex));
+ fftw_complex *freqspace1 = (fftw_complex *) fftw_malloc(w*h * sizeof(fftw_complex));
+ fftw_complex *freqspace2 = (fftw_complex *) 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
+ if(v < 1)
+ v = 1; /* do not write alpha zero */
+ 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)
+ {
+ fx = x;
+ fy = y;
+ if(fx > w/2)
+ fx -= w;
+ if(fy > h/2)
+ fy -= h;
+#ifdef DISCONTINUOUS
+ fx = sin(fx * TWO_PI / w);
+ fy = sin(fy * TWO_PI / h);
+#else
+#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);
+#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);
+#endif
+#endif
+#ifdef C99
+ 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
+ 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, const double *filter, int filterw, int filterh)
+{
+ int x, y;
+ double nx, ny, nz;
+ double v;
+ int i, j;
+ double *img_reduced = (double *) malloc(w*h * sizeof(double));
+
+ 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;
+ if(v < 1)
+ v = 1; /* do not write alpha zero */
+ 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;
+ nx = ny = 0;
+
+ for(i = -filterh / 2; i <= filterh / 2; ++i)
+ for(j = -filterw / 2; j <= filterw / 2; ++j)
+ {
+ nx += img_reduced[w*((y+i+h)%h)+(x+j+w)%w] * filter[(i + filterh / 2) * filterw + j + filterw / 2];
+ ny += img_reduced[w*((y+j+h)%h)+(x+i+w)%w] * filter[(i + filterh / 2) * filterw + j + filterw / 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);