make the algorithm a bit more customizable
[xonotic/xonotic.git] / misc / tools / fft-normalmap-to-heightmap.c
index dfbacde85c83e1ecbf5dc246cf9f82a452503432..909f591ff33667355fec0f52bebdcb5d77ab219d 100644 (file)
 
 #include <fftw3.h>
 
-void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
+#define TWO_PI (4*atan2(1,1) * 2)
+
+void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h, double scale, double offset)
 {
        int x, y;
+       int fx, fy;
+       int ffx, ffy;
        double nx, ny, nz;
        double v, vmin, vmax;
 #ifndef C99
@@ -46,9 +50,9 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
        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(w, h, imgspace1, freqspace1, FFTW_FORWARD, FFTW_ESTIMATE);
-       fftw_plan i22f2 = fftw_plan_dft_2d(w, h, imgspace2, freqspace2, FFTW_FORWARD, FFTW_ESTIMATE);
-       fftw_plan f12i1 = fftw_plan_dft_2d(w, h, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
+       fftw_plan i12f1 = fftw_plan_dft_2d(h, w, imgspace1, freqspace1, FFTW_FORWARD, FFTW_ESTIMATE);
+       fftw_plan i22f2 = fftw_plan_dft_2d(h, w, imgspace2, freqspace2, FFTW_FORWARD, FFTW_ESTIMATE);
+       fftw_plan f12i1 = fftw_plan_dft_2d(h, w, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
 
        for(y = 0; y < h; ++y)
        for(x = 0; x < w; ++x)
@@ -60,14 +64,14 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
                 * n_z = -dh/dh = -1
                 * BUT: darkplaces uses inverted normals, n_y actually is dh/dy by image pixel coordinates
                 */
-               nx = (int)map[(w*y+x)*4+2] - 127.5;
-               ny = (int)map[(w*y+x)*4+1] - 127.5;
-               nz = (int)map[(w*y+x)*4+0] - 127.5;
+               nx = ((int)map[(w*y+x)*4+2] - 127.5) / 128;
+               ny = ((int)map[(w*y+x)*4+1] - 127.5) / 128;
+               nz = ((int)map[(w*y+x)*4+0] - 127.5) / 128;
 
                /* reconstruct the derivatives from here */
 #ifdef C99
-               imgspace1[(w*y+x)] =  nx / nz; /* = dz/dx */
-               imgspace2[(w*y+x)] = -ny / nz; /* = dz/dy */
+               imgspace1[(w*y+x)] =  nx / nz * w; /* = dz/dx */
+               imgspace2[(w*y+x)] = -ny / nz * h; /* = dz/dy */
 #else
                imgspace1[(w*y+x)][0] =  nx / nz; /* = dz/dx */
                imgspace1[(w*y+x)][1] = 0;
@@ -84,23 +88,26 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
        for(y = 0; y < h; ++y)
        for(x = 0; x < w; ++x)
        {
-               int fx = x;
-               int fy = y;
+               fx = x;
+               fy = y;
                if(fx > w/2)
                        fx -= w;
                if(fy > h/2)
                        fy -= h;
+               /* these must have the same sign as fx and fy (so ffx*fx + ffy*fy is nonzero), otherwise do not matter */
+               ffx = fx;
+               ffy = fy;
 #ifdef C99
                if(fx||fy)
-                       freqspace1[(w*y+x)] = I * (fx * freqspace1[(w*y+x)] + fy * freqspace2[(w*y+x)]) / (fx*fx + fy*fy);
+                       freqspace1[(w*y+x)] = _Complex_I * (ffx * freqspace1[(w*y+x)] + ffy * freqspace2[(w*y+x)]) / (ffx*fx + ffy*fy) / TWO_PI;
                else
                        freqspace1[(w*y+x)] = 0;
 #else
                if(fx||fy)
                {
                        save = freqspace1[(w*y+x)][0];
-                       freqspace1[(w*y+x)][0] = -(fx * freqspace1[(w*y+x)][1] + fy * freqspace2[(w*y+x)][1]) / (fx*fx + fy*fy);
-                       freqspace1[(w*y+x)][1] = (fx * save + fy * freqspace2[(w*y+x)][0]) / (fx*fx + fy*fy);
+                       freqspace1[(w*y+x)][0] = -(ffx * freqspace1[(w*y+x)][1] + ffy * freqspace2[(w*y+x)][1]) / (ffx*fx + ffy*fy) / TWO_PI;
+                       freqspace1[(w*y+x)][1] =  (ffx * save + ffy * freqspace2[(w*y+x)][0]) / (ffx*fx + ffy*fy) / TWO_PI;
                }
                else
                {
@@ -112,41 +119,87 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
 
        fftw_execute(f12i1);
 
-       if(scale == 0)
+       /* renormalize, find min/max */
+       vmin = vmax = 0;
+       for(y = 0; y < h; ++y)
+       for(x = 0; x < w; ++x)
        {
 #ifdef C99
-               vmin = vmax = creal(imgspace1[0]);
+               v = creal(imgspace1[(w*y+x)] /= (w*h));
 #else
-               vmin = vmax = imgspace1[0][0];
+               v = (imgspace1[(w*y+x)][0] /= (w*h));
+               imgspace1[(w*y+x)][1] /= (w*h);
 #endif
+               if(v < vmin || (x == 0 && y == 0))
+                       vmin = v;
+               if(v > vmax || (x == 0 && y == 0))
+                       vmax = v;
+       }
+
+       if(refmap)
+       {
+               double f, a;
+               double o, s;
+               double sa, sfa, sffa, sfva, sva;
+               double mi, ma;
+               sa = sfa = sffa = sfva = sva = 0;
+               mi = 1;
+               ma = -1;
                for(y = 0; y < h; ++y)
                for(x = 0; x < w; ++x)
                {
+                       a = (int)refmap[(w*y+x)*4+3];
+                       v = (refmap[(w*y+x)*4+0]*0.114 + refmap[(w*y+x)*4+1]*0.587 + refmap[(w*y+x)*4+2]*0.299);
+                       v = (v - 128.0) / 127.0;
 #ifdef C99
-                       v = creal(imgspace1[(w*y+x)]);
+                       f = creal(imgspace1[(w*y+x)]);
 #else
-                       v = imgspace1[(w*y+x)][0];
+                       f = imgspace1[(w*y+x)][0];
 #endif
-                       if(v < vmin)
-                               vmin = v;
-                       if(v > vmax)
-                               vmax = v;
+                       if(a <= 0)
+                               continue;
+                       if(v < mi)
+                               mi = v;
+                       if(v > ma)
+                               ma = v;
+                       sa += a;
+                       sfa += f*a;
+                       sffa += f*f*a;
+                       sfva += f*v*a;
+                       sva += v*a;
+               }
+               if(mi < ma)
+               {
+                       /* linear regression ftw */
+                       o = (sfa*sfva - sffa*sva) / (sfa*sfa-sa*sffa);
+                       s = (sfa*sva - sa*sfva) / (sfa*sfa-sa*sffa);
+               }
+               else /* all values of v are equal, so we cannot get scale; we can still get offset */
+               {
+                       o = ((sva - sfa) / sa);
+                       s = 1;
                }
-
-               vmin /= (w*h);
-               vmax /= (w*h);
 
+               /*
+                * now apply user-given offset and scale to these values
+                * (x * s + o) * scale + offset
+                * x * s * scale + o * scale + offset
+                */
+               offset += o * scale;
+               scale *= s;
+       }
+       else if(scale == 0)
+       {
                /*
                 * map vmin to -1
                 * map vmax to +1
                 */
                scale = 2 / (vmax - vmin);
                offset = -(vmax + vmin) / (vmax - vmin);
-
-               printf("Autocomputed scale: %f\nAutocomputed offset: %f\n", scale, offset);
        }
 
-       scale /= (w*h);
+       printf("Min: %f\nAvg: %f\nMax: %f\nScale: %f\nOffset: %f\nScaled-Min: %f\nScaled-Avg: %f\nScaled-Max: %f\n", 
+               vmin, 0.0, vmax, scale, offset, vmin * scale + offset, offset, vmax * scale + offset);
 
        for(y = 0; y < h; ++y)
        for(x = 0; x < w; ++x)
@@ -178,7 +231,7 @@ 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, vmin, vmax;
+       double v;
 #ifndef C99
        double save;
 #endif
@@ -187,9 +240,9 @@ void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
        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(w, h, imgspace1, freqspace1, FFTW_FORWARD, FFTW_ESTIMATE);
-       fftw_plan f12i1 = fftw_plan_dft_2d(w, h, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
-       fftw_plan f22i2 = fftw_plan_dft_2d(w, h, freqspace2, imgspace2, FFTW_BACKWARD, FFTW_ESTIMATE);
+       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)
@@ -237,8 +290,8 @@ void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
                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)] = I * fy * freqspace1[(w*y+x)]; /* y derivative */
-               freqspace1[(w*y+x)] = I * fx * freqspace1[(w*y+x)]; /* x derivative */
+               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);
@@ -246,11 +299,11 @@ void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
                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] = -fy * freqspace1[(w*y+x)][1]; /* y derivative */
-               freqspace2[(w*y+x)][1] =  fy * freqspace1[(w*y+x)][0];
+               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] = -fx * freqspace1[(w*y+x)][1]; /* x derivative */
-               freqspace1[(w*y+x)][1] =  fx * save;
+               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
        }
 
@@ -269,15 +322,17 @@ void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
                nx = imgspace1[(w*y+x)][0];
                ny = imgspace2[(w*y+x)][0];
 #endif
-               nz = 1 / scale;
-               v = sqrt(nx*nx + ny*ny + nz*nz);
+               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(127.5 + 127.5 * nx);
-               map[(w*y+x)*4+1] = floor(127.5 + 127.5 * ny);
-               map[(w*y+x)*4+0] = floor(127.5 + 127.5 * nz);
+               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);
@@ -290,6 +345,69 @@ void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
        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);
+}
+
 unsigned char *FS_LoadFile(const char *fn, int *len)
 {
        unsigned char *buf = NULL;
@@ -773,22 +891,28 @@ int Image_WriteTGABGRA (const char *filename, int width, int height, const unsig
 
 int usage(const char *me)
 {
-       printf("Usage: %s <infile_norm.tga> <outfile_normandheight.tga> [<scale> [<offset>]] (get heightmap from normalmap)\n", me);
-       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -1 [<scale>] (read from R)\n", me);
-       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -2 [<scale>] (read from G)\n", me);
-       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -3 [<scale>] (read from R)\n", me);
-       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -4 [<scale>] (read from A)\n", me);
-       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -5 [<scale>] (read from (R+G+B)/3)\n", me);
-       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -6 [<scale>] (read from Y)\n", me);
+       printf("Usage: %s <infile_norm.tga> <outfile_normandheight.tga> [<scale> [<offset> [<infile_ref.tga>]]] (get heightmap from normalmap)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -1 [<scale>] (read from B, Diff)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -2 [<scale>] (read from G, Diff)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -3 [<scale>] (read from R, Diff)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -4 [<scale>] (read from A, Diff)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -5 [<scale>] (read from (R+G+B)/3, Diff)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -6 [<scale>] (read from Y, Diff)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -7 [<scale>] (read from B, FFT)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -8 [<scale>] (read from G, FFT)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -9 [<scale>] (read from R, FFT)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -10 [<scale>] (read from A, FFT)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -11 [<scale>] (read from (R+G+B)/3, FFT)\n", me);
+       printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> -12 [<scale>] (read from Y, FFT)\n", me);
        return 1;
 }
 
 int main(int argc, char **argv)
 {
-       const char *infile, *outfile;
+       const char *infile, *outfile, *reffile;
        double scale, offset;
-       int nmaplen;
-       unsigned char *nmapdata, *nmap;
+       int nmaplen, w, h;
+       unsigned char *nmapdata, *nmap, *refmap;
 
        if(argc > 1)
                infile = argv[1];
@@ -808,7 +932,12 @@ int main(int argc, char **argv)
        if(argc > 4)
                offset = atof(argv[4]);
        else
-               offset = 0;
+               offset = (scale<0) ? 1 : 0;
+
+       if(argc > 5)
+               reffile = argv[5];
+       else
+               reffile = NULL;
 
        nmapdata = FS_LoadFile(infile, &nmaplen);
        if(!nmapdata)
@@ -823,10 +952,39 @@ int main(int argc, char **argv)
                printf("LoadTGA_BGRA failed\n");
                return 2;
        }
-       if(scale < 0)
-               hmap_to_nmap(nmap, image_width, image_height, -scale-1, offset);
+       w = image_width;
+       h = image_height;
+
+       if(reffile)
+       {
+               nmapdata = FS_LoadFile(reffile, &nmaplen);
+               if(!nmapdata)
+               {
+                       printf("FS_LoadFile failed\n");
+                       return 2;
+               }
+               refmap = LoadTGA_BGRA(nmapdata, nmaplen);
+               free(nmapdata);
+               if(!refmap)
+               {
+                       printf("LoadTGA_BGRA failed\n");
+                       return 2;
+               }
+               if(image_width != w || image_height != h)
+               {
+                       printf("reference map must have same size as input normalmap\n");
+                       return 2;
+               }
+       }
+       else
+               refmap = NULL;
+
+       if(scale < -6)
+               hmap_to_nmap(nmap, image_width, image_height, -scale-7, offset);
+       else if(scale < 0)
+               hmap_to_nmap_local(nmap, image_width, image_height, -scale-1, offset);
        else
-               nmap_to_hmap(nmap, image_width, image_height, scale, offset);
+               nmap_to_hmap(nmap, refmap, image_width, image_height, scale, offset);
        if(!Image_WriteTGABGRA(outfile, image_width, image_height, nmap))
        {
                printf("Image_WriteTGABGRA failed\n");