diff --git a/pipi/Makefile.am b/pipi/Makefile.am index b525425..83e69a3 100644 --- a/pipi/Makefile.am +++ b/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) \ diff --git a/pipi/dither/dbs.c b/pipi/dither/dbs.c new file mode 100644 index 0000000..7f9ba85 --- /dev/null +++ b/pipi/dither/dbs.c @@ -0,0 +1,213 @@ +/* + * libpipi Proper image processing implementation library + * Copyright (c) 2004-2008 Sam Hocevar + * 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 + +#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; +} + diff --git a/pipi/pipi.h b/pipi/pipi.h index f06e982..d4b3fbf 100644 --- a/pipi/pipi.h +++ b/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 *);