/*
 *  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;
}