/*
 * Copyright 2012, Red Hat, Inc.
 * Copyright 2012, Soren Sandmann
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice (including the next
 * paragraph) shall be included in all copies or substantial portions of the
 * Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *
 * Author: Soren Sandmann <soren.sandmann@gmail.com>
 */
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#ifdef HAVE_CONFIG_H
#include <pixman-config.h>
#endif
#include "pixman-private.h"

typedef double (* kernel_func_t) (double x);

typedef struct
{
    pixman_kernel_t	kernel;
    kernel_func_t	func;
    double		width;
} filter_info_t;

static double
impulse_kernel (double x)
{
    return (x == 0.0)? 1.0 : 0.0;
}

static double
box_kernel (double x)
{
    return 1;
}

static double
linear_kernel (double x)
{
    return 1 - fabs (x);
}

static double
gaussian_kernel (double x)
{
#define SQRT2 (1.4142135623730950488016887242096980785696718753769480)
#define SIGMA (SQRT2 / 2.0)
    
    return exp (- x * x / (2 * SIGMA * SIGMA)) / (SIGMA * sqrt (2.0 * M_PI));
}

static double
sinc (double x)
{
    if (x == 0.0)
	return 1.0;
    else
	return sin (M_PI * x) / (M_PI * x);
}

static double
lanczos (double x, int n)
{
    return sinc (x) * sinc (x * (1.0 / n));
}

static double
lanczos2_kernel (double x)
{
    return lanczos (x, 2);
}

static double
lanczos3_kernel (double x)
{
    return lanczos (x, 3);
}

static double
nice_kernel (double x)
{
    return lanczos3_kernel (x * 0.75);
}

static double
general_cubic (double x, double B, double C)
{
    double ax = fabs(x);

    if (ax < 1)
    {
	return (((12 - 9 * B - 6 * C) * ax +
		 (-18 + 12 * B + 6 * C)) * ax * ax +
		(6 - 2 * B)) / 6;
    }
    else if (ax < 2)
    {
	return ((((-B - 6 * C) * ax +
		  (6 * B + 30 * C)) * ax +
		 (-12 * B - 48 * C)) * ax +
		(8 * B + 24 * C)) / 6;
    }
    else
    {
	return 0;
    }
}

static double
cubic_kernel (double x)
{
    /* This is the Mitchell-Netravali filter.
     *
     * (0.0, 0.5) would give us the Catmull-Rom spline,
     * but that one seems to be indistinguishable from Lanczos2.
     */
    return general_cubic (x, 1/3.0, 1/3.0);
}

static const filter_info_t filters[] =
{
    { PIXMAN_KERNEL_IMPULSE,	        impulse_kernel,   0.0 },
    { PIXMAN_KERNEL_BOX,	        box_kernel,       1.0 },
    { PIXMAN_KERNEL_LINEAR,	        linear_kernel,    2.0 },
    { PIXMAN_KERNEL_CUBIC,		cubic_kernel,     4.0 },
    { PIXMAN_KERNEL_GAUSSIAN,	        gaussian_kernel,  5.0 },
    { PIXMAN_KERNEL_LANCZOS2,	        lanczos2_kernel,  4.0 },
    { PIXMAN_KERNEL_LANCZOS3,	        lanczos3_kernel,  6.0 },
    { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel,      8.0 },
};

/* This function scales @kernel2 by @scale, then
 * aligns @x1 in @kernel1 with @x2 in @kernel2 and
 * and integrates the product of the kernels across @width.
 *
 * This function assumes that the intervals are within
 * the kernels in question. E.g., the caller must not
 * try to integrate a linear kernel ouside of [-1:1]
 */
static double
integral (pixman_kernel_t kernel1, double x1,
	  pixman_kernel_t kernel2, double scale, double x2,
	  double width)
{
    if (kernel1 == PIXMAN_KERNEL_BOX && kernel2 == PIXMAN_KERNEL_BOX)
    {
	return width;
    }
    /* The LINEAR filter is not differentiable at 0, so if the
     * integration interval crosses zero, break it into two
     * separate integrals.
     */
    else if (kernel1 == PIXMAN_KERNEL_LINEAR && x1 < 0 && x1 + width > 0)
    {
	return
	    integral (kernel1, x1, kernel2, scale, x2, - x1) +
	    integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
    }
    else if (kernel2 == PIXMAN_KERNEL_LINEAR && x2 < 0 && x2 + width > 0)
    {
	return
	    integral (kernel1, x1, kernel2, scale, x2, - x2) +
	    integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2);
    }
    else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
    {
	assert (width == 0.0);
	return filters[kernel2].func (x2 * scale);
    }
    else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
    {
	assert (width == 0.0);
	return filters[kernel1].func (x1);
    }
    else
    {
	/* Integration via Simpson's rule
	 * See http://www.intmath.com/integration/6-simpsons-rule.php
	 * 12 segments (6 cubic approximations) seems to produce best
	 * result for lanczos3.linear, which was the combination that
	 * showed the most errors.  This makes sense as the lanczos3
	 * filter is 6 wide.
	 */
#define N_SEGMENTS 12
#define SAMPLE(a1, a2)							\
	(filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
	
	double s = 0.0;
	double h = width / N_SEGMENTS;
	int i;

	s = SAMPLE (x1, x2);

	for (i = 1; i < N_SEGMENTS; i += 2)
	{
	    double a1 = x1 + h * i;
	    double a2 = x2 + h * i;
	    s += 4 * SAMPLE (a1, a2);
	}

	for (i = 2; i < N_SEGMENTS; i += 2)
	{
	    double a1 = x1 + h * i;
	    double a2 = x2 + h * i;
	    s += 2 * SAMPLE (a1, a2);
	}

	s += SAMPLE (x1 + width, x2 + width);
	
	return h * s * (1.0 / 3.0);
    }
}

static void
create_1d_filter (int              width,
		  pixman_kernel_t  reconstruct,
		  pixman_kernel_t  sample,
		  double           scale,
		  int              n_phases,
		  pixman_fixed_t *pstart,
		  pixman_fixed_t *pend
		  )
{
    pixman_fixed_t *p = pstart;
    double step;
    int i;
    if(width <= 0) return;
    step = 1.0 / n_phases;

    for (i = 0; i < n_phases; ++i)
    {
        double frac = step / 2.0 + i * step;
	pixman_fixed_t new_total;
        int x, x1, x2;
	double total, e;

	/* Sample convolution of reconstruction and sampling
	 * filter. See rounding.txt regarding the rounding
	 * and sample positions.
	 */

	x1 = ceil (frac - width / 2.0 - 0.5);
	x2 = x1 + width;
    assert( p >= pstart && p + (x2 - x1) <= pend ); /* assert validity of the following loop */
	total = 0;
        for (x = x1; x < x2; ++x)
        {
	    double pos = x + 0.5 - frac;
	    double rlow = - filters[reconstruct].width / 2.0;
	    double rhigh = rlow + filters[reconstruct].width;
	    double slow = pos - scale * filters[sample].width / 2.0;
	    double shigh = slow + scale * filters[sample].width;
	    double c = 0.0;
	    double ilow, ihigh;

	    if (rhigh >= slow && rlow <= shigh)
	    {
		ilow = MAX (slow, rlow);
		ihigh = MIN (shigh, rhigh);

		c = integral (reconstruct, ilow,
			      sample, 1.0 / scale, ilow - pos,
			      ihigh - ilow);
	    }

            *p = (pixman_fixed_t)floor (c * 65536.0 + 0.5);
	    total += *p;
	    p++;
        }

	/* Normalize, with error diffusion */
	p -= width;
	assert(p >= pstart && p + (x2 - x1) <= pend); /* assert validity of the following loop */

    total = 65536.0 / total;
    new_total = 0;
	e = 0.0;
	for (x = x1; x < x2; ++x)
	{
	    double v = (*p) * total + e;
	    pixman_fixed_t t = floor (v + 0.5);

	    e = v - t;
	    new_total += t;
	    *p++ = t;
	}

	/* pixman_fixed_e's worth of error may remain; put it
	 * at the first sample, since that is the only one that
	 * hasn't had any error diffused into it.
	 */

	assert(p - width >= pstart && p - width < pend); /* assert... */
	*(p - width) += pixman_fixed_1 - new_total;
    }
}


static int
filter_width (pixman_kernel_t reconstruct, pixman_kernel_t sample, double size)
{
    return ceil (filters[reconstruct].width + size * filters[sample].width);
}

#ifdef PIXMAN_GNUPLOT

/* If enable-gnuplot is configured, then you can pipe the output of a
 * pixman-using program to gnuplot and get a continuously-updated plot
 * of the horizontal filter. This works well with demos/scale to test
 * the filter generation.
 *
 * The plot is all the different subposition filters shuffled
 * together. This is misleading in a few cases:
 *
 *  IMPULSE.BOX - goes up and down as the subfilters have different
 *		  numbers of non-zero samples
 *  IMPULSE.TRIANGLE - somewhat crooked for the same reason
 *  1-wide filters - looks triangular, but a 1-wide box would be more
 *		     accurate
 */
static void
gnuplot_filter (int width, int n_phases, const pixman_fixed_t* p)
{
    double step;
    int i, j;
    int first;

    step = 1.0 / n_phases;

    printf ("set style line 1 lc rgb '#0060ad' lt 1 lw 0.5 pt 7 pi 1 ps 0.5\n");
    printf ("plot [x=%g:%g] '-' with linespoints ls 1\n", -width*0.5, width*0.5);
    /* Print a point at the origin so that y==0 line is included: */
    printf ("0 0\n\n");

    /* The position of the first sample of the phase corresponding to
     * frac is given by:
     * 
     *     ceil (frac - width / 2.0 - 0.5) + 0.5 - frac
     * 
     * We have to find the frac that minimizes this expression.
     * 
     * For odd widths, we have
     * 
     *     ceil (frac - width / 2.0 - 0.5) + 0.5 - frac
     *   = ceil (frac) + K - frac
     *   = 1 + K - frac
     * 
     * for some K, so this is minimized when frac is maximized and
     * strictly growing with frac. So for odd widths, we can simply
     * start at the last phase and go backwards.
     * 
     * For even widths, we have
     * 
     *     ceil (frac - width / 2.0 - 0.5) + 0.5 - frac
     *   = ceil (frac - 0.5) + K - frac
     * 
     * The graph for this function (ignoring K) looks like this:
     * 
     *        0.5
     *           |    |\ 
     *           |    | \ 
     *           |    |  \ 
     *         0 |    |   \ 
     *           |\   |
     *           | \  |
     *           |  \ |
     *      -0.5 |   \|
     *   ---------------------------------
     *           0    0.5   1
     * 
     * So in this case we need to start with the phase whose frac is
     * less than, but as close as possible to 0.5, then go backwards
     * until we hit the first phase, then wrap around to the last
     * phase and continue backwards.
     * 
     * Which phase is as close as possible 0.5? The locations of the
     * sampling point corresponding to the kth phase is given by
     * 1/(2 * n_phases) + k / n_phases:
     * 
     *         1/(2 * n_phases) + k / n_phases = 0.5
     *  
     * from which it follows that
     * 
     *         k = (n_phases - 1) / 2
     * 
     * rounded down is the phase in question.
     */
    if (width & 1)
	first = n_phases - 1;
    else
	first = (n_phases - 1) / 2;

    for (j = 0; j < width; ++j)
    {
	for (i = 0; i < n_phases; ++i)
	{
	    int phase = first - i;
	    double frac, pos;

	    if (phase < 0)
		phase = n_phases + phase;

	    frac = step / 2.0 + phase * step;
	    pos = ceil (frac - width / 2.0 - 0.5) + 0.5 - frac + j;

	    printf ("%g %g\n",
		    pos,
		    pixman_fixed_to_double (*(p + phase * width + j)));
	}
    }

    printf ("e\n");
    fflush (stdout);
}

#endif

/* Create the parameter list for a SEPARABLE_CONVOLUTION filter
 * with the given kernels and scale parameters
 */
PIXMAN_EXPORT pixman_fixed_t *
pixman_filter_create_separable_convolution (int             *n_values,
					    pixman_fixed_t   scale_x,
					    pixman_fixed_t   scale_y,
					    pixman_kernel_t  reconstruct_x,
					    pixman_kernel_t  reconstruct_y,
					    pixman_kernel_t  sample_x,
					    pixman_kernel_t  sample_y,
					    int              subsample_bits_x,
					    int	             subsample_bits_y)
{
    double sx = fabs (pixman_fixed_to_double (scale_x));
    double sy = fabs (pixman_fixed_to_double (scale_y));
    pixman_fixed_t *params;
    int subsample_x, subsample_y;
    int width, height;

    width = filter_width (reconstruct_x, sample_x, sx);
    subsample_x = (1 << subsample_bits_x);

    height = filter_width (reconstruct_y, sample_y, sy);
    subsample_y = (1 << subsample_bits_y);

    *n_values = 4 + width * subsample_x + height * subsample_y;
    
    params = malloc (*n_values * sizeof (pixman_fixed_t));
    if (!params)
	return NULL;

    params[0] = pixman_int_to_fixed (width);
    params[1] = pixman_int_to_fixed (height);
    params[2] = pixman_int_to_fixed (subsample_bits_x);
    params[3] = pixman_int_to_fixed (subsample_bits_y);

    {
        pixman_fixed_t
            *xparams = params+4,
            *yparams = xparams + width*subsample_x,
            *endparams = params + *n_values;
        create_1d_filter(width, reconstruct_x, sample_x, sx, subsample_x,
                         xparams, yparams);
        create_1d_filter(height, reconstruct_y, sample_y, sy, subsample_y,
                         yparams, endparams);
    }

#ifdef PIXMAN_GNUPLOT
    gnuplot_filter(width, subsample_x, params + 4);
#endif

    return params;
}