* dbs.c: improve the DBS human visual system kernel by adding two Gaussian

kernels.

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

pipi/dither/dbs.c

@@ -24,46 +24,35 @@
 #include "pipi.h"
 #include "pipi_internals.h"
 
-#define N 5
+#define N 7
 #define NN ((N * 2 + 1))
 
-/* FIXME: use a better HVS than a simple gaussian; adding two gaussians
- * seemed to be a lot more efficient. Get rid of pipi_gaussian_blur calls
- * at the same time. */
 /* FIXME: though the algorithm is supposed to stop, we do not have a real,
  * guaranteed stop condition here. */
 
-static void makegauss(double mat[NN][NN], double sigma)
+pipi_image_t *pipi_dbs(pipi_image_t *src)
 {
-    double t = 0;
-    int i, j;
-
-    sigma = 2. * sigma * sigma;
+    double kernel[NN * NN];
+    double t = 0.;
+    pipi_image_t *dst, *tmp1, *tmp2;
+    pipi_pixels_t *srcp, *dstp, *tmp1p, *tmp2p;
+    float *srcdata, *dstdata, *tmp1data, *tmp2data;
+    int i, j, w, h;
 
     for(j = 0; j < NN; j++)
         for(i = 0; i < NN; i++)
         {
             double a = (double)(i - N);
             double b = (double)(j - N);
-            mat[i][j] = pow(M_E, - (a * a + b * b) / sigma);
-            t += mat[i][j];
+            kernel[j * NN + i] =
+                    1.0 * exp(-(a * a + b * b) / (2. * 1.5 * 1.5))
+                  + 0.1 * exp(-(a * a + b * b) / (2. * 0.6 * 0.6));
+            t += kernel[j * NN + i];
         }
 
     for(j = 0; j < NN; j++)
         for(i = 0; i < NN; i++)
-            mat[i][j] /= t;
-}
-
-pipi_image_t *pipi_dbs(pipi_image_t *src)
-{
-    double mat[NN][NN];
-    double sigma = 1.2;
-    pipi_image_t *dst, *tmp1, *tmp2;
-    pipi_pixels_t *srcp, *dstp, *tmp1p, *tmp2p;
-    float *srcdata, *dstdata, *tmp1data, *tmp2data;
-    int w, h;
-
-    makegauss(mat, sigma);
+            kernel[j * NN + i] /= t;
 
     w = src->w;
     h = src->h;
@@ -71,25 +60,25 @@ pipi_image_t *pipi_dbs(pipi_image_t *src)
     srcp = pipi_getpixels(src, PIPI_PIXELS_Y_F);
     srcdata = (float *)srcp->pixels;
 
-    tmp1 = pipi_gaussian_blur(src, sigma);
+    tmp1 = pipi_convolution(src, NN, NN, kernel);
     tmp1p = pipi_getpixels(tmp1, PIPI_PIXELS_Y_F);
     tmp1data = (float *)tmp1p->pixels;
 
     /* The initial dither is an empty image. So is its blurred version,
-     * but I leave the pipi_gaussian_blur() call here in case we choose
+     * but I leave the pipi_convolution() call here in case we choose
      * to change the way to create the initial dither. */
     dst = pipi_new(w, h);
     dstp = pipi_getpixels(dst, PIPI_PIXELS_Y_F);
     dstdata = (float *)dstp->pixels;
 
-    tmp2 = pipi_gaussian_blur(dst, sigma);
+    tmp2 = pipi_convolution(dst, NN, NN, kernel);
     tmp2p = pipi_getpixels(tmp2, PIPI_PIXELS_Y_F);
     tmp2data = (float *)tmp2p->pixels;
 
     for(;;)
     {
         int changes = 0;
-        int x, y, i, j, n;
+        int x, y, n;
 
         for(y = 0; y < h; y++)
         for(x = 0; x < w; x++)
@@ -113,7 +102,7 @@ pipi_image_t *pipi_dbs(pipi_image_t *src)
                     if(x + i < 0 || x + i >= w)
                         continue;
 
-                    m = mat[i + N][j + N];
+                    m = kernel[(j + N) * NN + i + N];
                     p = tmp1data[(y + j) * w + x + i];
                     q1 = tmp2data[(y + j) * w + x + i];
                     q2 = q1 - m * d + m * (1. - d);
@@ -152,8 +141,8 @@ pipi_image_t *pipi_dbs(pipi_image_t *src)
                             continue;
                         if(i - idx + N < 0 || i - idx + N >= NN)
                             continue;
-                        ma = mat[i + N][j + N];
-                        mb = mat[i - idx + N][j - idy + N];
+                        ma = kernel[(j + N) * NN + i + N];
+                        mb = kernel[(j - idy + N) * NN + i - idx + N];
                         p = tmp1data[(y + j) * w + x + i];
                         q1 = tmp2data[(y + j) * w + x + i];
                         q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
@@ -182,17 +171,17 @@ pipi_image_t *pipi_dbs(pipi_image_t *src)
             for(j = -N; j < N + 1; j++)
             for(i = -N; i < N + 1; i++)
             {
-                double m = mat[i + N][j + N];
+                double m = kernel[(j + N) * NN + i + N];
                 if(y + j >= 0 && y + j < h
                     && x + i >= 0 && x + i < w)
                 {
-                    double t = tmp2data[(y + j) * w + x + i];
+                    t = tmp2data[(y + j) * w + x + i];
                     tmp2data[(y + j) * w + x + i] = t + m * (d2 - d);
                 }
                 if((opx || opy) && y + opy + j >= 0 && y + opy + j < h
                                 && x + opx + i >= 0 && x + opx + i < w)
                 {
-                    double t = tmp2data[(y + opy + j) * w + x + opx + i];
+                    t = tmp2data[(y + opy + j) * w + x + opx + i];
                     tmp2data[(y + opy + j) * w + x + opx + i]
                             = t + m * (d - d2);
                 }