From 2a3c7ca9903b7131e552f792942e278d7df20b43 Mon Sep 17 00:00:00 2001 From: sam Date: Sun, 3 Aug 2008 18:36:17 +0000 Subject: [PATCH] * 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 | 57 +++++++++++++++++++---------------------------- 1 file changed, 23 insertions(+), 34 deletions(-) diff --git a/pipi/dither/dbs.c b/pipi/dither/dbs.c index 7f9ba85..a046514 100644 --- a/pipi/dither/dbs.c +++ b/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); }