make the algorithm a bit more customizable
[xonotic/xonotic.git] / misc / tools / fft-normalmap-to-heightmap.c
index 95782ec83cf1ea8d6175d47496800eb1e7fcafb3..909f591ff33667355fec0f52bebdcb5d77ab219d 100644 (file)
@@ -38,6 +38,8 @@
 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
@@ -48,9 +50,9 @@ void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h,
        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)
@@ -86,23 +88,26 @@ void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h,
        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) / TWO_PI;
+                       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) / TWO_PI;
-                       freqspace1[(w*y+x)][1] =  (fx * save + fy * freqspace2[(w*y+x)][0]) / (fx*fx + fy*fy) / TWO_PI;
+                       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
                {
@@ -114,24 +119,25 @@ void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h,
 
        fftw_execute(f12i1);
 
-       if(refmap)
+       /* renormalize, find min/max */
+       vmin = vmax = 0;
+       for(y = 0; y < h; ++y)
+       for(x = 0; x < w; ++x)
        {
-               // refmap: a reference map to define the heights
-               // alpha = weight, color = value
-               // if more than one color value is used, colors are also matched
-
-               // we do linear regression, basically
-               // f'(x, y) = f(x, y) * scale + offset
-               // sum((f(x, y) * scale + offset - ref_y(x, y))^2 * ref_a(x, y)) minimize
-
-               // diff by offset:
-               // sum(-2*ref_y(x,y)*ref_a(x,y) + 2*scale*f(x,y)*ref_a(x,y) + 2*offset*ref_a(x,y)) = 0
-               // diff by scale:
-               // sum(-2*f(x,y)*ref_a(x,y) + 2*scale*f(x,y)^2*ref_a(x,y) + 2*offset*f(x,y)*ref_a(x,y)) = 0
-               // ->
-               // offset = (sfa*sfya - sffa*sya) / (sfa*sfa-sa*sffa)
-               // scale  = (sfa*sya - sa*sfya) / (sfa*sfa-sa*sffa)
+#ifdef C99
+               v = creal(imgspace1[(w*y+x)] /= (w*h));
+#else
+               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;
@@ -142,9 +148,9 @@ void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h,
                for(y = 0; y < h; ++y)
                for(x = 0; x < w; ++x)
                {
-                       a = (int)refmap[(w*y+x)*4+0];
-                       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);
-                       v = (v - 128.0) / 127.0; // value 0 is forbidden, 1 -> -1, 255 -> 1
+                       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
                        f = creal(imgspace1[(w*y+x)]);
 #else
@@ -152,71 +158,48 @@ void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h,
 #endif
                        if(a <= 0)
                                continue;
-                       if(y < mi)
-                               mi = y;
-                       if(y > ma)
-                               ma = y;
+                       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;
                }
-               sfa /= (w*h);
-               sffa /= (w*h);
-               sffa /= (w*h);
-               sfva /= (w*h);
                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
+               else /* all values of v are equal, so we cannot get scale; we can still get offset */
                {
                        o = ((sva - sfa) / sa);
                        s = 1;
                }
-               // now apply user-given offset and scale to these values
-               // (x * s + o) * scale + offset
-               // x * s * scale + o * scale + offset
+
+               /*
+                * 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)
        {
-#ifdef C99
-               vmin = vmax = creal(imgspace1[0]);
-#else
-               vmin = vmax = imgspace1[0][0];
-#endif
-               for(y = 0; y < h; ++y)
-               for(x = 0; x < w; ++x)
-               {
-#ifdef C99
-                       v = creal(imgspace1[(w*y+x)]);
-#else
-                       v = imgspace1[(w*y+x)][0];
-#endif
-                       if(v < vmin)
-                               vmin = v;
-                       if(v > vmax)
-                               vmax = v;
-               }
-
-               vmin /= (w*h);
-               vmax /= (w*h);
-
                /*
                 * 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)
@@ -257,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)
@@ -307,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)] = TWO_PI*I * fy * freqspace1[(w*y+x)]; /* y derivative */
-               freqspace1[(w*y+x)] = TWO_PI*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);
@@ -974,7 +957,7 @@ int main(int argc, char **argv)
 
        if(reffile)
        {
-               nmapdata = FS_LoadFile(infile, &nmaplen);
+               nmapdata = FS_LoadFile(reffile, &nmaplen);
                if(!nmapdata)
                {
                        printf("FS_LoadFile failed\n");
@@ -993,6 +976,8 @@ int main(int argc, char **argv)
                        return 2;
                }
        }
+       else
+               refmap = NULL;
 
        if(scale < -6)
                hmap_to_nmap(nmap, image_width, image_height, -scale-7, offset);