1/*
2 * Copyright 2012, Red Hat, Inc.
3 * Copyright 2012, Soren Sandmann
4 *
5 * Permission is hereby granted, free of charge, to any person obtaining a
6 * copy of this software and associated documentation files (the "Software"),
7 * to deal in the Software without restriction, including without limitation
8 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
9 * and/or sell copies of the Software, and to permit persons to whom the
10 * Software is furnished to do so, subject to the following conditions:
11 *
12 * The above copyright notice and this permission notice (including the next
13 * paragraph) shall be included in all copies or substantial portions of the
14 * Software.
15 *
16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
19 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
21 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
22 * DEALINGS IN THE SOFTWARE.
23 *
24 * Author: Soren Sandmann <soren.sandmann@gmail.com>
25 */
26#include <string.h>
27#include <stdlib.h>
28#include <stdio.h>
29#include <math.h>
30#include <assert.h>
31#include <config.h>
32#include "pixman-private.h"
33
34typedef double (* kernel_func_t) (double x);
35
36typedef struct
37{
38    pixman_kernel_t	kernel;
39    kernel_func_t	func;
40    double		width;
41} filter_info_t;
42
43static double
44impulse_kernel (double x)
45{
46    return (x == 0.0)? 1.0 : 0.0;
47}
48
49static double
50box_kernel (double x)
51{
52    return 1;
53}
54
55static double
56linear_kernel (double x)
57{
58    return 1 - fabs (x);
59}
60
61static double
62gaussian_kernel (double x)
63{
64#define SQRT2 (1.4142135623730950488016887242096980785696718753769480)
65#define SIGMA (SQRT2 / 2.0)
66
67    return exp (- x * x / (2 * SIGMA * SIGMA)) / (SIGMA * sqrt (2.0 * M_PI));
68}
69
70static double
71sinc (double x)
72{
73    if (x == 0.0)
74	return 1.0;
75    else
76	return sin (M_PI * x) / (M_PI * x);
77}
78
79static double
80lanczos (double x, int n)
81{
82    return sinc (x) * sinc (x * (1.0 / n));
83}
84
85static double
86lanczos2_kernel (double x)
87{
88    return lanczos (x, 2);
89}
90
91static double
92lanczos3_kernel (double x)
93{
94    return lanczos (x, 3);
95}
96
97static double
98nice_kernel (double x)
99{
100    return lanczos3_kernel (x * 0.75);
101}
102
103static double
104general_cubic (double x, double B, double C)
105{
106    double ax = fabs(x);
107
108    if (ax < 1)
109    {
110	return ((12 - 9 * B - 6 * C) * ax * ax * ax +
111		(-18 + 12 * B + 6 * C) * ax * ax + (6 - 2 * B)) / 6;
112    }
113    else if (ax >= 1 && ax < 2)
114    {
115	return ((-B - 6 * C) * ax * ax * ax +
116		(6 * B + 30 * C) * ax * ax + (-12 * B - 48 * C) *
117		ax + (8 * B + 24 * C)) / 6;
118    }
119    else
120    {
121	return 0;
122    }
123}
124
125static double
126cubic_kernel (double x)
127{
128    /* This is the Mitchell-Netravali filter.
129     *
130     * (0.0, 0.5) would give us the Catmull-Rom spline,
131     * but that one seems to be indistinguishable from Lanczos2.
132     */
133    return general_cubic (x, 1/3.0, 1/3.0);
134}
135
136static const filter_info_t filters[] =
137{
138    { PIXMAN_KERNEL_IMPULSE,	        impulse_kernel,   0.0 },
139    { PIXMAN_KERNEL_BOX,	        box_kernel,       1.0 },
140    { PIXMAN_KERNEL_LINEAR,	        linear_kernel,    2.0 },
141    { PIXMAN_KERNEL_CUBIC,		cubic_kernel,     4.0 },
142    { PIXMAN_KERNEL_GAUSSIAN,	        gaussian_kernel,  6 * SIGMA },
143    { PIXMAN_KERNEL_LANCZOS2,	        lanczos2_kernel,  4.0 },
144    { PIXMAN_KERNEL_LANCZOS3,	        lanczos3_kernel,  6.0 },
145    { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel,      8.0 },
146};
147
148/* This function scales @kernel2 by @scale, then
149 * aligns @x1 in @kernel1 with @x2 in @kernel2 and
150 * and integrates the product of the kernels across @width.
151 *
152 * This function assumes that the intervals are within
153 * the kernels in question. E.g., the caller must not
154 * try to integrate a linear kernel ouside of [-1:1]
155 */
156static double
157integral (pixman_kernel_t kernel1, double x1,
158	  pixman_kernel_t kernel2, double scale, double x2,
159	  double width)
160{
161    /* If the integration interval crosses zero, break it into
162     * two separate integrals. This ensures that filters such
163     * as LINEAR that are not differentiable at 0 will still
164     * integrate properly.
165     */
166    if (x1 < 0 && x1 + width > 0)
167    {
168	return
169	    integral (kernel1, x1, kernel2, scale, x2, - x1) +
170	    integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
171    }
172    else if (x2 < 0 && x2 + width > 0)
173    {
174	return
175	    integral (kernel1, x1, kernel2, scale, x2, - x2) +
176	    integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2);
177    }
178    else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
179    {
180	assert (width == 0.0);
181	return filters[kernel2].func (x2 * scale);
182    }
183    else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
184    {
185	assert (width == 0.0);
186	return filters[kernel1].func (x1);
187    }
188    else
189    {
190	/* Integration via Simpson's rule */
191#define N_SEGMENTS 128
192#define SAMPLE(a1, a2)							\
193	(filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
194
195	double s = 0.0;
196	double h = width / (double)N_SEGMENTS;
197	int i;
198
199	s = SAMPLE (x1, x2);
200
201	for (i = 1; i < N_SEGMENTS; i += 2)
202	{
203	    double a1 = x1 + h * i;
204	    double a2 = x2 + h * i;
205
206	    s += 2 * SAMPLE (a1, a2);
207
208	    if (i >= 2 && i < N_SEGMENTS - 1)
209		s += 4 * SAMPLE (a1, a2);
210	}
211
212	s += SAMPLE (x1 + width, x2 + width);
213
214	return h * s * (1.0 / 3.0);
215    }
216}
217
218static pixman_fixed_t *
219create_1d_filter (int             *width,
220		  pixman_kernel_t  reconstruct,
221		  pixman_kernel_t  sample,
222		  double           scale,
223		  int              n_phases)
224{
225    pixman_fixed_t *params, *p;
226    double step;
227    double size;
228    int i;
229
230    size = scale * filters[sample].width + filters[reconstruct].width;
231    *width = ceil (size);
232
233    p = params = malloc (*width * n_phases * sizeof (pixman_fixed_t));
234    if (!params)
235        return NULL;
236
237    step = 1.0 / n_phases;
238
239    for (i = 0; i < n_phases; ++i)
240    {
241        double frac = step / 2.0 + i * step;
242	pixman_fixed_t new_total;
243        int x, x1, x2;
244	double total;
245
246	/* Sample convolution of reconstruction and sampling
247	 * filter. See rounding.txt regarding the rounding
248	 * and sample positions.
249	 */
250
251	x1 = ceil (frac - *width / 2.0 - 0.5);
252        x2 = x1 + *width;
253
254	total = 0;
255        for (x = x1; x < x2; ++x)
256        {
257	    double pos = x + 0.5 - frac;
258	    double rlow = - filters[reconstruct].width / 2.0;
259	    double rhigh = rlow + filters[reconstruct].width;
260	    double slow = pos - scale * filters[sample].width / 2.0;
261	    double shigh = slow + scale * filters[sample].width;
262	    double c = 0.0;
263	    double ilow, ihigh;
264
265	    if (rhigh >= slow && rlow <= shigh)
266	    {
267		ilow = MAX (slow, rlow);
268		ihigh = MIN (shigh, rhigh);
269
270		c = integral (reconstruct, ilow,
271			      sample, 1.0 / scale, ilow - pos,
272			      ihigh - ilow);
273	    }
274
275	    total += c;
276            *p++ = (pixman_fixed_t)(c * 65535.0 + 0.5);
277        }
278
279	/* Normalize */
280	p -= *width;
281        total = 1 / total;
282        new_total = 0;
283	for (x = x1; x < x2; ++x)
284	{
285	    pixman_fixed_t t = (*p) * total + 0.5;
286
287	    new_total += t;
288	    *p++ = t;
289	}
290
291	if (new_total != pixman_fixed_1)
292	    *(p - *width / 2) += (pixman_fixed_1 - new_total);
293    }
294
295    return params;
296}
297
298/* Create the parameter list for a SEPARABLE_CONVOLUTION filter
299 * with the given kernels and scale parameters
300 */
301PIXMAN_EXPORT pixman_fixed_t *
302pixman_filter_create_separable_convolution (int             *n_values,
303					    pixman_fixed_t   scale_x,
304					    pixman_fixed_t   scale_y,
305					    pixman_kernel_t  reconstruct_x,
306					    pixman_kernel_t  reconstruct_y,
307					    pixman_kernel_t  sample_x,
308					    pixman_kernel_t  sample_y,
309					    int              subsample_bits_x,
310					    int	             subsample_bits_y)
311{
312    double sx = fabs (pixman_fixed_to_double (scale_x));
313    double sy = fabs (pixman_fixed_to_double (scale_y));
314    pixman_fixed_t *horz = NULL, *vert = NULL, *params = NULL;
315    int subsample_x, subsample_y;
316    int width, height;
317
318    subsample_x = (1 << subsample_bits_x);
319    subsample_y = (1 << subsample_bits_y);
320
321    horz = create_1d_filter (&width, reconstruct_x, sample_x, sx, subsample_x);
322    vert = create_1d_filter (&height, reconstruct_y, sample_y, sy, subsample_y);
323
324    if (!horz || !vert)
325        goto out;
326
327    *n_values = 4 + width * subsample_x + height * subsample_y;
328
329    params = malloc (*n_values * sizeof (pixman_fixed_t));
330    if (!params)
331        goto out;
332
333    params[0] = pixman_int_to_fixed (width);
334    params[1] = pixman_int_to_fixed (height);
335    params[2] = pixman_int_to_fixed (subsample_bits_x);
336    params[3] = pixman_int_to_fixed (subsample_bits_y);
337
338    memcpy (params + 4, horz,
339	    width * subsample_x * sizeof (pixman_fixed_t));
340    memcpy (params + 4 + width * subsample_x, vert,
341	    height * subsample_y * sizeof (pixman_fixed_t));
342
343out:
344    free (horz);
345    free (vert);
346
347    return params;
348}
349