/*
 *  edd           error diffusion displacement
 *  Copyright (c) 2008 Sam Hocevar <sam@zoy.org>
 *                All Rights Reserved
 *
 *  $Id$
 *
 *  This program 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.
 */

/* This program computes the Floyd-Steinberg error diffusion algorithm's
 * displacement on the given input image. Error diffusion displacement is
 * introduced in the paper "Reinstating Floyd-Steinberg: Improved Metrics
 * for Quality Assessment of Error Diffusion Algorithms.", ICISP 2008
 * Proceedings. Lecture Notes in Computer Science 5099 Springer 2008, ISBN
 * 978-3-540-69904-0.
 *
 * The resulting dx/dy values are usually around 0.16/0.26 for images that
 * are not entirely black and white. */

#include "config.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <pipi.h>

#define Z 3

int main(int argc, char *argv[])
{
    double sigma = 1.2, precision = 0.001, step = 2.;
    double best = 1., fx = -1., fy = -1., bfx = 0., bfy = 0.;
    double e, e0, e1;
    pipi_image_t *img, *kernel, *gauss, *dither, *tmp;
    int dx, dy;

    if(argc < 2)
    {
        fprintf(stderr, "%s: too few arguments\n", argv[0]);
        fprintf(stderr, "Usage: %s <image>\n", argv[0]);
        return EXIT_FAILURE;
    }

    /* Load image, convert it to grayscale, dither it with Floyd-Steinberg */
    img = pipi_load(argv[1]);
    pipi_set_colorspace(img, PIPI_PIXELS_Y_F32);
    gauss = pipi_gaussian_blur(img, sigma);
    kernel = pipi_load("ediff:fs");
    dither = pipi_dither_ediff(img, kernel, PIPI_SCAN_RASTER);
    pipi_free(kernel);
    pipi_free(img);

    /* Compute the standard error */
    tmp = pipi_gaussian_blur(dither, sigma);
    e0 = pipi_measure_msd(gauss, tmp);
    pipi_free(tmp);

    /* Compute the fast error */
    tmp = pipi_gaussian_blur_ext(dither, sigma, sigma, 0.0, 0.16, 0.26);
    e1 = pipi_measure_msd(gauss, tmp);
    pipi_free(tmp);

    /* Compute displacement */
    while(step > precision)
    {
        for(dy = 0; dy <= Z; dy++)
            for(dx = 0; dx <= Z; dx++)
            {
                tmp = pipi_gaussian_blur_ext(dither, sigma, sigma, 0.0,
                                             fx + step * dx / Z,
                                             fy + step * dy / Z);
                e = pipi_measure_msd(gauss, tmp);
                pipi_free(tmp);
                if(e < best)
                {
                    best = e;
                    bfx = fx + step * dx / Z;
                    bfy = fy + step * dy / Z;
                }
            }

        fx = bfx - step / Z;
        fy = bfy - step / Z;
        step = step * 2 / Z;
    }

    printf("E: %g E_fast: %g E_min: %g dx: %g dy: %g\n", e0, e1, best, fx, fy);

    pipi_free(dither);
    pipi_free(gauss);

    return 0;
}