* Various small optimisations that make us gain ~3% in speed. Not really

worth it, but since it's also better style, I'm keeping them.


git-svn-id: file:///srv/caca.zoy.org/var/lib/svn/research@2276 92316355-f0b4-4df1-b90c-862c8a59935f

2008-displacement/Makefile

@@ -4,7 +4,7 @@ xy2d: xy2d.c
 	$(CC) -W -Wall $^ -o $@
 
 main: main.c
-	$(CC) -Wall -O3 -funroll-loops -W -Wall $^ -o $@ $$(pkg-config --cflags --libs sdl) -lSDL_image -lm
+	$(CC) -Wall -O3 -funroll-loops -ffast-math -W -Wall $^ -o $@ $$(pkg-config --cflags --libs sdl) -lSDL_image -lm
 
 paper.pdf: paper.tex fs-histo.tex
 	latex paper.tex

2008-displacement/main.c

@@ -15,9 +15,9 @@
 #define true 1
 #define false 0
 
-int WIDTH, HEIGHT;
+static int WIDTH, HEIGHT;
 
-static inline float get(float *src, int x, int y)
+static inline float get(float const *src, int x, int y)
 {
     return src[y * WIDTH + x];
 }
@@ -32,7 +32,7 @@ static float *new(void)
     return malloc(WIDTH * HEIGHT * sizeof(float));
 }
 
-static float *copy(float *src)
+static float *copy(float const *src)
 {
     float *dest = malloc(WIDTH * HEIGHT * sizeof(float));
     memcpy(dest, src, WIDTH * HEIGHT * sizeof(float)); 
@@ -63,7 +63,7 @@ static void makegauss(float mat[NN][NN], float sigma, float dx, float dy)
             mat[i][j] /= t;
 }
 
-static float *gauss(float *src, float mat[NN][NN])
+static float *gauss(float const *src, float mat[NN][NN])
 {
     float *dest = new();
     int x, y, i, j;
@@ -83,7 +83,7 @@ static float *gauss(float *src, float mat[NN][NN])
     return dest;
 }
 
-static float *fullgauss(float *src, float mat[NN][NN])
+static float *fullgauss(float const *src, float mat[NN][NN])
 {
     float *dest = new();
     int x, y, i, j;
@@ -105,7 +105,7 @@ static float *fullgauss(float *src, float mat[NN][NN])
     return dest;
 }
 
-static float fulldist(float *p1, float *p2)
+static float fulldist(float const *p1, const float *p2)
 {
     float error = 0.;
     int x, y;
@@ -120,22 +120,32 @@ static float fulldist(float *p1, float *p2)
     return error / (WIDTH * HEIGHT);
 }
 
-static float dist(float *p1, float *p2)
+static float dist(float const *p1, float const *p2, float max)
 {
     float error = 0.;
     int x, y;
 
+    max *= (WIDTH - N) * (HEIGHT - N);
+
     for(y = N; y < HEIGHT - N; y++)
+    {
         for(x = N; x < WIDTH - N; x++)
         {
             float t = get(p1, x, y) - get(p2, x, y);
             error += t * t;
         }
 
+        /* Experience shows that this is almost useless, because the
+         * functions we manipulate are so small that the difference
+         * only happens in the last 1% of the image... */
+        if(error > max)
+            break;
+    }
+
     return error / ((WIDTH - N) * (HEIGHT - N));
 }
 
-static float *load(char *name)
+static float *load(char const *name)
 {
     SDL_Surface *tmp, *surface;
     uint32_t *pixels;
@@ -168,7 +178,7 @@ static float *load(char *name)
     return floats;
 }
 
-static void save(float *src, char *name)
+static void save(float const *src, char const *name)
 {
     SDL_Surface *surface;
     uint32_t *pixels;
@@ -189,7 +199,7 @@ static void save(float *src, char *name)
     SDL_SaveBMP(surface, name);
 }
 
-static float *ostromoukhov(float *src)
+static float *ostromoukhov(float const *src)
 {
     static int const table[][3] =
     {
@@ -275,7 +285,7 @@ static float *ostromoukhov(float *src)
      c d e f g
      h i j k l
 */
-static float *ed(float *src, int serpentine,
+static float *ed(float const *src, int serpentine,
                  int a, int b, int c, int d, int e, int f,
                  int g, int h, int i, int j, int k, int l)
 {
@@ -356,7 +366,8 @@ static float *ed(float *src, int serpentine,
 }
 
 /* XXX */
-static float *dbs(float *src, float *orig, float sigma, float dx, float dy)
+static float *dbs(float const *src, float const *orig,
+                  float sigma, float dx, float dy)
 {
     float mat[NN][NN];
     float *dest, *tmp, *tmp2;
@@ -369,7 +380,7 @@ static float *dbs(float *src, float *orig, float sigma, float dx, float dy)
     dest = copy(orig);
     tmp2 = fullgauss(dest, mat);
 
-    error = dist(tmp, tmp2);
+    error = dist(tmp, tmp2, 1.);
 
     for(;;)
     {
@@ -497,33 +508,30 @@ static float *dbs(float *src, float *orig, float sigma, float dx, float dy)
     return dest;
 }
 
-static void study(float *src, float *dest, float sigma, float precision)
+static void study(float const *src, float const *dest,
+                  float sigma, float precision)
 {
 #   define Z 3
     float mat[NN][NN];
     float *tmp, *tmp2;
-    float e0, best = 0., fx = -1., fy = -1., step = 2.;
+    float e, e0, best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.;
     int dx, dy;
 
     makegauss(mat, sigma, 0., 0.);
     tmp = gauss(src, mat);
     tmp2 = gauss(dest, mat);
 
-    e0 = dist(tmp, tmp2);
+    e0 = dist(tmp, tmp2, 1.);
     free(tmp2);
 
     while(step > precision)
     {
-        float bfx = 0., bfy = 0., e;
-
-        best = 9999.;
-
         for(dy = 0; dy <= Z; dy++)
             for(dx = 0; dx <= Z; dx++)
             {
                 makegauss(mat, sigma, fx + step * dx / Z, fy + step * dy / Z);
                 tmp2 = gauss(dest, mat);
-                e = dist(tmp, tmp2);
+                e = dist(tmp, tmp2, best);
                 free(tmp2);
                 if(e < best)
                 {
@@ -545,7 +553,7 @@ static void study(float *src, float *dest, float sigma, float precision)
     fflush(stdout);
 }
 
-static float *merge(float *im1, float *im2, float t)
+static float *merge(float const *im1, float const *im2, float t)
 {
     float *dest = new();
     int x, y;
@@ -569,11 +577,8 @@ static void usage(char *argv[])
 
 int main(int argc, char *argv[])
 {
-    float mat0[NN][NN];
-    float mat[NN][NN];
     float *src;
-    float sigma;
-    int mode, dx, dy, i;
+    int mode, i;
 
     if(argc < 2)
     {
@@ -630,6 +635,7 @@ int main(int argc, char *argv[])
         case 3:
         case 4:
         {
+            float mat0[NN][NN];
             float *dest, *tmp, *tmp2;
             int a, b, c, d, e;
 
@@ -670,12 +676,12 @@ int main(int argc, char *argv[])
                 {
                     tmp = gauss(src, mat0);
                     tmp2 = gauss(dest, mat0);
-                    printf("E = %.5g\n", 1000. * dist(tmp, tmp2));
+                    printf("E = %.5g\n", 1000. * dist(tmp, tmp2, 1.));
                     free(tmp);
                     free(tmp2);
                 }
                 else
-                    study(src, dest, 1.2, 0.01);
+                    study(src, dest, 1.2, 0.001);
                 fflush(stdout);
                 free(dest);
             }
@@ -728,7 +734,7 @@ int main(int argc, char *argv[])
 
             makegauss(mat, 1.2, fx, fy);
             tmp2 = gauss(dest, mat);
-            printf("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2));
+            printf("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.));
             fflush(stdout);
             free(tmp2);
         }
@@ -750,7 +756,7 @@ int main(int argc, char *argv[])
         dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0);
         tmp2 = gauss(dest, 0., 0.);
         printf("%.5g: (%02d %02d %02d %02d)\n",
-               1000. * dist(tmp, tmp2), a, b, c, d);
+               1000. * dist(tmp, tmp2, 1.), a, b, c, d);
         free(dest);
         free(tmp2);
     }
@@ -765,41 +771,41 @@ int main(int argc, char *argv[])
 
     dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
     tmp2 = gauss(dest, mat0);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     makegauss(mat, 1.2, 0.16, 0.28);
     tmp2 = gauss(dest, mat);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     free(dest);
 
     dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
     tmp2 = gauss(dest, mat0);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     makegauss(mat, 1.2, 0.26, 0.76);
     tmp2 = gauss(dest, mat);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     free(dest);
 
     dest = ostromoukhov(src);
     tmp2 = gauss(dest, mat0);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     makegauss(mat, 1.2, 0.0, 0.19);
     tmp2 = gauss(dest, mat);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     free(dest);
 
     dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0);
     tmp2 = gauss(dest, mat0);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     makegauss(mat, 1.2, 0.0, 0.34);
     tmp2 = gauss(dest, mat);
-    printf("%.5g ", 1000. * dist(tmp, tmp2));
+    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
     free(tmp2);
     free(dest);