* dbs.c: new dithering function: Direct Binary Search. One of the best

dithering methods as of today, though extremely slow.

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

pipi/Makefile.am

@@ -30,6 +30,7 @@ libpipi_la_SOURCES = \
 	$(codec_sources) \
 	filter/blur.c \
 	dither/floydsteinberg.c \
+	dither/dbs.c \
 	$(NULL)
 libpipi_la_CFLAGS = $(codec_cflags)
 libpipi_la_LDFLAGS = $(codec_libs) \

pipi/dither/dbs.c

@@ -0,0 +1,213 @@
+/*
+ *  libpipi       Proper image processing implementation library
+ *  Copyright (c) 2004-2008 Sam Hocevar <sam@zoy.org>
+ *                All Rights Reserved
+ *
+ *  $Id$
+ *
+ *  This library is free software. It comes without any warranty, to
+ *  the extent permitted by applicable law. You can redistribute it
+ *  and/or modify it under the terms of the Do What The Fuck You Want
+ *  To Public License, Version 2, as published by Sam Hocevar. See
+ *  http://sam.zoy.org/wtfpl/COPYING for more details.
+ */
+
+/*
+ * dbs.c: Direct Binary Search dithering functions
+ */
+
+#include "config.h"
+#include "common.h"
+
+#include <math.h>
+
+#include "pipi.h"
+#include "pipi_internals.h"
+
+#define N 5
+#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)
+{
+    double t = 0;
+    int i, j;
+
+    sigma = 2. * sigma * sigma;
+
+    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];
+        }
+
+    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);
+
+    w = src->w;
+    h = src->h;
+
+    srcp = pipi_getpixels(src, PIPI_PIXELS_Y_F);
+    srcdata = (float *)srcp->pixels;
+
+    tmp1 = pipi_gaussian_blur(src, sigma);
+    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
+     * 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);
+    tmp2p = pipi_getpixels(tmp2, PIPI_PIXELS_Y_F);
+    tmp2data = (float *)tmp2p->pixels;
+
+    for(;;)
+    {
+        int changes = 0;
+        int x, y, i, j, n;
+
+        for(y = 0; y < h; y++)
+        for(x = 0; x < w; x++)
+        {
+            double d, d2, e, best = 0.;
+            int opx = -1, opy = -1;
+
+            d = dstdata[y * w + x];
+
+            /* Compute the effect of a toggle */
+            e = 0.;
+            for(j = -N; j < N + 1; j++)
+            {
+                if(y + j < 0 || y + j >= h)
+                    continue;
+
+                for(i = -N; i < N + 1; i++)
+                {
+                    double m, p, q1, q2;
+
+                    if(x + i < 0 || x + i >= w)
+                        continue;
+
+                    m = mat[i + N][j + N];
+                    p = tmp1data[(y + j) * w + x + i];
+                    q1 = tmp2data[(y + j) * w + x + i];
+                    q2 = q1 - m * d + m * (1. - d);
+                    e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
+                }
+            }
+            if(e > best)
+            {
+                best = e;
+                opx = opy = 0;
+            }
+
+            /* Compute the effect of swaps */
+            for(n = 0; n < 8; n++)
+            {
+                static int const step[] =
+                    { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
+                int idx = step[n * 2], idy = step[n * 2 + 1];
+                if(y + idy < 0 || y + idy >= h
+                     || x + idx < 0 || x + idx >= w)
+                    continue;
+                d2 = dstdata[(y + idy) * w + x + idx];
+                if(d2 == d)
+                    continue;
+                e = 0.;
+                for(j = -N; j < N + 1; j++)
+                {
+                    if(y + j < 0 || y + j >= h)
+                        continue;
+                    if(j - idy + N < 0 || j - idy + N >= NN)
+                        continue;
+                    for(i = -N; i < N + 1; i++)
+                    {
+                        double ma, mb, p, q1, q2;
+                        if(x + i < 0 || x + i >= w)
+                            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];
+                        p = tmp1data[(y + j) * w + x + i];
+                        q1 = tmp2data[(y + j) * w + x + i];
+                        q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
+                        e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
+                    }
+                }
+                if(e > best)
+                {
+                    best = e;
+                    opx = idx;
+                    opy = idy;
+                }
+            }
+
+            /* Apply the change if interesting */
+            if(best <= 0.)
+                continue;
+            if(opx || opy)
+            {
+                d2 = dstdata[(y + opy) * w + x + opx];
+                dstdata[(y + opy) * w + x + opx] = d;
+            }
+            else
+                d2 = 1. - d;
+            dstdata[y * w + x] = d2;
+            for(j = -N; j < N + 1; j++)
+            for(i = -N; i < N + 1; i++)
+            {
+                double m = mat[i + N][j + N];
+                if(y + j >= 0 && y + j < h
+                    && x + i >= 0 && x + i < w)
+                {
+                    double 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];
+                    tmp2data[(y + opy + j) * w + x + opx + i]
+                            = t + m * (d - d2);
+                }
+            }
+
+            changes++;
+        }
+
+        if(changes == 0)
+            break;
+    }
+
+    pipi_free(tmp1);
+    pipi_free(tmp2);
+
+    return dst;
+}
+

pipi/pipi.h

@@ -69,6 +69,7 @@ extern pipi_image_t *pipi_gaussian_blur_ext(pipi_image_t *,
                                             float, float, float, float);
 
 extern pipi_image_t *pipi_floydsteinberg(pipi_image_t *);
+extern pipi_image_t *pipi_dbs(pipi_image_t *);
 extern void pipi_dither_24to16(pipi_image_t *);
 
 extern void pipi_test(pipi_image_t *);