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 #ifdef HAVE_CONFIG_H
32 #include <config.h>
33 #endif
34 #include "pixman-private.h"
35
36 typedef double (* kernel_func_t) (double x);
37
38 typedef struct
39 {
40 pixman_kernel_t kernel;
41 kernel_func_t func;
42 double width;
43 } filter_info_t;
44
45 static double
impulse_kernel(double x)46 impulse_kernel (double x)
47 {
48 return (x == 0.0)? 1.0 : 0.0;
49 }
50
51 static double
box_kernel(double x)52 box_kernel (double x)
53 {
54 return 1;
55 }
56
57 static double
linear_kernel(double x)58 linear_kernel (double x)
59 {
60 return 1 - fabs (x);
61 }
62
63 static double
gaussian_kernel(double x)64 gaussian_kernel (double x)
65 {
66 #define SQRT2 (1.4142135623730950488016887242096980785696718753769480)
67 #define SIGMA (SQRT2 / 2.0)
68
69 return exp (- x * x / (2 * SIGMA * SIGMA)) / (SIGMA * sqrt (2.0 * M_PI));
70 }
71
72 static double
sinc(double x)73 sinc (double x)
74 {
75 if (x == 0.0)
76 return 1.0;
77 else
78 return sin (M_PI * x) / (M_PI * x);
79 }
80
81 static double
lanczos(double x,int n)82 lanczos (double x, int n)
83 {
84 return sinc (x) * sinc (x * (1.0 / n));
85 }
86
87 static double
lanczos2_kernel(double x)88 lanczos2_kernel (double x)
89 {
90 return lanczos (x, 2);
91 }
92
93 static double
lanczos3_kernel(double x)94 lanczos3_kernel (double x)
95 {
96 return lanczos (x, 3);
97 }
98
99 static double
nice_kernel(double x)100 nice_kernel (double x)
101 {
102 return lanczos3_kernel (x * 0.75);
103 }
104
105 static double
general_cubic(double x,double B,double C)106 general_cubic (double x, double B, double C)
107 {
108 double ax = fabs(x);
109
110 if (ax < 1)
111 {
112 return ((12 - 9 * B - 6 * C) * ax * ax * ax +
113 (-18 + 12 * B + 6 * C) * ax * ax + (6 - 2 * B)) / 6;
114 }
115 else if (ax >= 1 && ax < 2)
116 {
117 return ((-B - 6 * C) * ax * ax * ax +
118 (6 * B + 30 * C) * ax * ax + (-12 * B - 48 * C) *
119 ax + (8 * B + 24 * C)) / 6;
120 }
121 else
122 {
123 return 0;
124 }
125 }
126
127 static double
cubic_kernel(double x)128 cubic_kernel (double x)
129 {
130 /* This is the Mitchell-Netravali filter.
131 *
132 * (0.0, 0.5) would give us the Catmull-Rom spline,
133 * but that one seems to be indistinguishable from Lanczos2.
134 */
135 return general_cubic (x, 1/3.0, 1/3.0);
136 }
137
138 static const filter_info_t filters[] =
139 {
140 { PIXMAN_KERNEL_IMPULSE, impulse_kernel, 0.0 },
141 { PIXMAN_KERNEL_BOX, box_kernel, 1.0 },
142 { PIXMAN_KERNEL_LINEAR, linear_kernel, 2.0 },
143 { PIXMAN_KERNEL_CUBIC, cubic_kernel, 4.0 },
144 { PIXMAN_KERNEL_GAUSSIAN, gaussian_kernel, 6 * SIGMA },
145 { PIXMAN_KERNEL_LANCZOS2, lanczos2_kernel, 4.0 },
146 { PIXMAN_KERNEL_LANCZOS3, lanczos3_kernel, 6.0 },
147 { PIXMAN_KERNEL_LANCZOS3_STRETCHED, nice_kernel, 8.0 },
148 };
149
150 /* This function scales @kernel2 by @scale, then
151 * aligns @x1 in @kernel1 with @x2 in @kernel2 and
152 * and integrates the product of the kernels across @width.
153 *
154 * This function assumes that the intervals are within
155 * the kernels in question. E.g., the caller must not
156 * try to integrate a linear kernel ouside of [-1:1]
157 */
158 static double
integral(pixman_kernel_t kernel1,double x1,pixman_kernel_t kernel2,double scale,double x2,double width)159 integral (pixman_kernel_t kernel1, double x1,
160 pixman_kernel_t kernel2, double scale, double x2,
161 double width)
162 {
163 /* If the integration interval crosses zero, break it into
164 * two separate integrals. This ensures that filters such
165 * as LINEAR that are not differentiable at 0 will still
166 * integrate properly.
167 */
168 if (x1 < 0 && x1 + width > 0)
169 {
170 return
171 integral (kernel1, x1, kernel2, scale, x2, - x1) +
172 integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
173 }
174 else if (x2 < 0 && x2 + width > 0)
175 {
176 return
177 integral (kernel1, x1, kernel2, scale, x2, - x2) +
178 integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2);
179 }
180 else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
181 {
182 assert (width == 0.0);
183 return filters[kernel2].func (x2 * scale);
184 }
185 else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
186 {
187 assert (width == 0.0);
188 return filters[kernel1].func (x1);
189 }
190 else
191 {
192 /* Integration via Simpson's rule */
193 #define N_SEGMENTS 128
194 #define SAMPLE(a1, a2) \
195 (filters[kernel1].func ((a1)) * filters[kernel2].func ((a2) * scale))
196
197 double s = 0.0;
198 double h = width / (double)N_SEGMENTS;
199 int i;
200
201 s = SAMPLE (x1, x2);
202
203 for (i = 1; i < N_SEGMENTS; i += 2)
204 {
205 double a1 = x1 + h * i;
206 double a2 = x2 + h * i;
207
208 s += 2 * SAMPLE (a1, a2);
209
210 if (i >= 2 && i < N_SEGMENTS - 1)
211 s += 4 * SAMPLE (a1, a2);
212 }
213
214 s += SAMPLE (x1 + width, x2 + width);
215
216 return h * s * (1.0 / 3.0);
217 }
218 }
219
220 static pixman_fixed_t *
create_1d_filter(int * width,pixman_kernel_t reconstruct,pixman_kernel_t sample,double scale,int n_phases)221 create_1d_filter (int *width,
222 pixman_kernel_t reconstruct,
223 pixman_kernel_t sample,
224 double scale,
225 int n_phases)
226 {
227 pixman_fixed_t *params, *p;
228 double step;
229 double size;
230 int i;
231
232 size = scale * filters[sample].width + filters[reconstruct].width;
233 *width = ceil (size);
234
235 p = params = malloc (*width * n_phases * sizeof (pixman_fixed_t));
236 if (!params)
237 return NULL;
238
239 step = 1.0 / n_phases;
240
241 for (i = 0; i < n_phases; ++i)
242 {
243 double frac = step / 2.0 + i * step;
244 pixman_fixed_t new_total;
245 int x, x1, x2;
246 double total;
247
248 /* Sample convolution of reconstruction and sampling
249 * filter. See rounding.txt regarding the rounding
250 * and sample positions.
251 */
252
253 x1 = ceil (frac - *width / 2.0 - 0.5);
254 x2 = x1 + *width;
255
256 total = 0;
257 for (x = x1; x < x2; ++x)
258 {
259 double pos = x + 0.5 - frac;
260 double rlow = - filters[reconstruct].width / 2.0;
261 double rhigh = rlow + filters[reconstruct].width;
262 double slow = pos - scale * filters[sample].width / 2.0;
263 double shigh = slow + scale * filters[sample].width;
264 double c = 0.0;
265 double ilow, ihigh;
266
267 if (rhigh >= slow && rlow <= shigh)
268 {
269 ilow = MAX (slow, rlow);
270 ihigh = MIN (shigh, rhigh);
271
272 c = integral (reconstruct, ilow,
273 sample, 1.0 / scale, ilow - pos,
274 ihigh - ilow);
275 }
276
277 total += c;
278 *p++ = (pixman_fixed_t)(c * 65535.0 + 0.5);
279 }
280
281 /* Normalize */
282 p -= *width;
283 total = 1 / total;
284 new_total = 0;
285 for (x = x1; x < x2; ++x)
286 {
287 pixman_fixed_t t = (*p) * total + 0.5;
288
289 new_total += t;
290 *p++ = t;
291 }
292
293 if (new_total != pixman_fixed_1)
294 *(p - *width / 2) += (pixman_fixed_1 - new_total);
295 }
296
297 return params;
298 }
299
300 /* Create the parameter list for a SEPARABLE_CONVOLUTION filter
301 * with the given kernels and scale parameters
302 */
303 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)304 pixman_filter_create_separable_convolution (int *n_values,
305 pixman_fixed_t scale_x,
306 pixman_fixed_t scale_y,
307 pixman_kernel_t reconstruct_x,
308 pixman_kernel_t reconstruct_y,
309 pixman_kernel_t sample_x,
310 pixman_kernel_t sample_y,
311 int subsample_bits_x,
312 int subsample_bits_y)
313 {
314 double sx = fabs (pixman_fixed_to_double (scale_x));
315 double sy = fabs (pixman_fixed_to_double (scale_y));
316 pixman_fixed_t *horz = NULL, *vert = NULL, *params = NULL;
317 int subsample_x, subsample_y;
318 int width, height;
319
320 subsample_x = (1 << subsample_bits_x);
321 subsample_y = (1 << subsample_bits_y);
322
323 horz = create_1d_filter (&width, reconstruct_x, sample_x, sx, subsample_x);
324 vert = create_1d_filter (&height, reconstruct_y, sample_y, sy, subsample_y);
325
326 if (!horz || !vert)
327 goto out;
328
329 *n_values = 4 + width * subsample_x + height * subsample_y;
330
331 params = malloc (*n_values * sizeof (pixman_fixed_t));
332 if (!params)
333 goto out;
334
335 params[0] = pixman_int_to_fixed (width);
336 params[1] = pixman_int_to_fixed (height);
337 params[2] = pixman_int_to_fixed (subsample_bits_x);
338 params[3] = pixman_int_to_fixed (subsample_bits_y);
339
340 memcpy (params + 4, horz,
341 width * subsample_x * sizeof (pixman_fixed_t));
342 memcpy (params + 4 + width * subsample_x, vert,
343 height * subsample_y * sizeof (pixman_fixed_t));
344
345 out:
346 free (horz);
347 free (vert);
348
349 return params;
350 }
351