Changeset 2663 for libpipi/trunk/pipi/dither/dbs.c

Timestamp:

08/03/08 20:36:17 (5 years ago)

Author:

sam

Message:

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

File:

• libpipi/trunk/pipi/dither/dbs.c (modified) (9 diffs)

libpipi/trunk/pipi/dither/dbs.c

@@ -25 +25 @@
 #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++)
@@ -46 +45 @@
             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;
@@ -72 +61 @@
     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);
@@ -83 +72 @@
     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;
@@ -90 +79 @@
     {
         int changes = 0;
-        int x, y, i, j, n;
+        int x, y, n;
 
         for(y = 0; y < h; y++)
@@ -114 +103 @@
                         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];
@@ -153 +142 @@
                         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];
@@ -183 +172 @@
             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);
                 }
@@ -193 +182 @@
                                 && 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);