Line data Source code
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
46 0 : impulse_kernel (double x)
47 : {
48 0 : return (x == 0.0)? 1.0 : 0.0;
49 : }
50 :
51 : static double
52 0 : box_kernel (double x)
53 : {
54 0 : return 1;
55 : }
56 :
57 : static double
58 0 : linear_kernel (double x)
59 : {
60 0 : return 1 - fabs (x);
61 : }
62 :
63 : static double
64 0 : gaussian_kernel (double x)
65 : {
66 : #define SQRT2 (1.4142135623730950488016887242096980785696718753769480)
67 : #define SIGMA (SQRT2 / 2.0)
68 :
69 0 : return exp (- x * x / (2 * SIGMA * SIGMA)) / (SIGMA * sqrt (2.0 * M_PI));
70 : }
71 :
72 : static double
73 0 : sinc (double x)
74 : {
75 0 : if (x == 0.0)
76 0 : return 1.0;
77 : else
78 0 : return sin (M_PI * x) / (M_PI * x);
79 : }
80 :
81 : static double
82 0 : lanczos (double x, int n)
83 : {
84 0 : return sinc (x) * sinc (x * (1.0 / n));
85 : }
86 :
87 : static double
88 0 : lanczos2_kernel (double x)
89 : {
90 0 : return lanczos (x, 2);
91 : }
92 :
93 : static double
94 0 : lanczos3_kernel (double x)
95 : {
96 0 : return lanczos (x, 3);
97 : }
98 :
99 : static double
100 0 : nice_kernel (double x)
101 : {
102 0 : return lanczos3_kernel (x * 0.75);
103 : }
104 :
105 : static double
106 0 : general_cubic (double x, double B, double C)
107 : {
108 0 : double ax = fabs(x);
109 :
110 0 : if (ax < 1)
111 : {
112 0 : return ((12 - 9 * B - 6 * C) * ax * ax * ax +
113 0 : (-18 + 12 * B + 6 * C) * ax * ax + (6 - 2 * B)) / 6;
114 : }
115 0 : else if (ax >= 1 && ax < 2)
116 : {
117 0 : return ((-B - 6 * C) * ax * ax * ax +
118 0 : (6 * B + 30 * C) * ax * ax + (-12 * B - 48 * C) *
119 0 : ax + (8 * B + 24 * C)) / 6;
120 : }
121 : else
122 : {
123 0 : return 0;
124 : }
125 : }
126 :
127 : static double
128 0 : 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 0 : 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
159 0 : 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 0 : if (x1 < 0 && x1 + width > 0)
169 : {
170 : return
171 0 : integral (kernel1, x1, kernel2, scale, x2, - x1) +
172 0 : integral (kernel1, 0, kernel2, scale, x2 - x1, width + x1);
173 : }
174 0 : else if (x2 < 0 && x2 + width > 0)
175 : {
176 : return
177 0 : integral (kernel1, x1, kernel2, scale, x2, - x2) +
178 0 : integral (kernel1, x1 - x2, kernel2, scale, 0, width + x2);
179 : }
180 0 : else if (kernel1 == PIXMAN_KERNEL_IMPULSE)
181 : {
182 0 : assert (width == 0.0);
183 0 : return filters[kernel2].func (x2 * scale);
184 : }
185 0 : else if (kernel2 == PIXMAN_KERNEL_IMPULSE)
186 : {
187 0 : assert (width == 0.0);
188 0 : 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 0 : double s = 0.0;
198 0 : double h = width / (double)N_SEGMENTS;
199 : int i;
200 :
201 0 : s = SAMPLE (x1, x2);
202 :
203 0 : for (i = 1; i < N_SEGMENTS; i += 2)
204 : {
205 0 : double a1 = x1 + h * i;
206 0 : double a2 = x2 + h * i;
207 :
208 0 : s += 2 * SAMPLE (a1, a2);
209 :
210 0 : if (i >= 2 && i < N_SEGMENTS - 1)
211 0 : s += 4 * SAMPLE (a1, a2);
212 : }
213 :
214 0 : s += SAMPLE (x1 + width, x2 + width);
215 :
216 0 : return h * s * (1.0 / 3.0);
217 : }
218 : }
219 :
220 : static pixman_fixed_t *
221 0 : 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 0 : size = scale * filters[sample].width + filters[reconstruct].width;
233 0 : *width = ceil (size);
234 :
235 0 : p = params = malloc (*width * n_phases * sizeof (pixman_fixed_t));
236 0 : if (!params)
237 0 : return NULL;
238 :
239 0 : step = 1.0 / n_phases;
240 :
241 0 : for (i = 0; i < n_phases; ++i)
242 : {
243 0 : 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 0 : x1 = ceil (frac - *width / 2.0 - 0.5);
254 0 : x2 = x1 + *width;
255 :
256 0 : total = 0;
257 0 : for (x = x1; x < x2; ++x)
258 : {
259 0 : double pos = x + 0.5 - frac;
260 0 : double rlow = - filters[reconstruct].width / 2.0;
261 0 : double rhigh = rlow + filters[reconstruct].width;
262 0 : double slow = pos - scale * filters[sample].width / 2.0;
263 0 : double shigh = slow + scale * filters[sample].width;
264 0 : double c = 0.0;
265 : double ilow, ihigh;
266 :
267 0 : if (rhigh >= slow && rlow <= shigh)
268 : {
269 0 : ilow = MAX (slow, rlow);
270 0 : ihigh = MIN (shigh, rhigh);
271 :
272 0 : c = integral (reconstruct, ilow,
273 : sample, 1.0 / scale, ilow - pos,
274 : ihigh - ilow);
275 : }
276 :
277 0 : total += c;
278 0 : *p++ = (pixman_fixed_t)(c * 65535.0 + 0.5);
279 : }
280 :
281 : /* Normalize */
282 0 : p -= *width;
283 0 : total = 1 / total;
284 0 : new_total = 0;
285 0 : for (x = x1; x < x2; ++x)
286 : {
287 0 : pixman_fixed_t t = (*p) * total + 0.5;
288 :
289 0 : new_total += t;
290 0 : *p++ = t;
291 : }
292 :
293 0 : if (new_total != pixman_fixed_1)
294 0 : *(p - *width / 2) += (pixman_fixed_1 - new_total);
295 : }
296 :
297 0 : 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 *
304 0 : 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 0 : double sx = fabs (pixman_fixed_to_double (scale_x));
315 0 : double sy = fabs (pixman_fixed_to_double (scale_y));
316 0 : pixman_fixed_t *horz = NULL, *vert = NULL, *params = NULL;
317 : int subsample_x, subsample_y;
318 : int width, height;
319 :
320 0 : subsample_x = (1 << subsample_bits_x);
321 0 : subsample_y = (1 << subsample_bits_y);
322 :
323 0 : horz = create_1d_filter (&width, reconstruct_x, sample_x, sx, subsample_x);
324 0 : vert = create_1d_filter (&height, reconstruct_y, sample_y, sy, subsample_y);
325 :
326 0 : if (!horz || !vert)
327 : goto out;
328 :
329 0 : *n_values = 4 + width * subsample_x + height * subsample_y;
330 :
331 0 : params = malloc (*n_values * sizeof (pixman_fixed_t));
332 0 : if (!params)
333 0 : goto out;
334 :
335 0 : params[0] = pixman_int_to_fixed (width);
336 0 : params[1] = pixman_int_to_fixed (height);
337 0 : params[2] = pixman_int_to_fixed (subsample_bits_x);
338 0 : params[3] = pixman_int_to_fixed (subsample_bits_y);
339 :
340 0 : memcpy (params + 4, horz,
341 0 : width * subsample_x * sizeof (pixman_fixed_t));
342 0 : memcpy (params + 4 + width * subsample_x, vert,
343 0 : height * subsample_y * sizeof (pixman_fixed_t));
344 :
345 : out:
346 0 : free (horz);
347 0 : free (vert);
348 :
349 0 : return params;
350 : }
|
