1 /* scale.c - image scaling
2 *
3 * Raster graphics library
4 *
5 * Copyright (c) 1997-2003 Alfredo K. Kojima
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Library General Public
9 * License as published by the Free Software Foundation; either
10 * version 2 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Library General Public License for more details.
16 *
17 * You should have received a copy of the GNU Library General Public
18 * License along with this library; if not, write to the Free
19 * Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
20 * MA 02110-1301, USA.
21 */
22
23 #include <config.h>
24
25 #include <stdlib.h>
26 #include <stdio.h>
27 #include <string.h>
28 #include <X11/Xlib.h>
29 #include <math.h>
30 #include <assert.h>
31
32 #include "wraster.h"
33 #include "scale.h"
34
35 /*
36 *----------------------------------------------------------------------
37 * RScaleImage--
38 * Creates a scaled copy of an image.
39 *
40 * Returns:
41 * The new scaled image.
42 *
43 *----------------------------------------------------------------------
44 */
RScaleImage(RImage * image,unsigned new_width,unsigned new_height)45 RImage *RScaleImage(RImage * image, unsigned new_width, unsigned new_height)
46 {
47 int ox;
48 int px, py;
49 register int x, y, t;
50 int dx, dy;
51 unsigned char *s;
52 unsigned char *d;
53 RImage *img;
54
55 if (image == NULL)
56 return NULL;
57
58 if (new_width == image->width && new_height == image->height)
59 return RCloneImage(image);
60
61 img = RCreateImage(new_width, new_height, image->format == RRGBAFormat);
62
63 if (!img)
64 return NULL;
65
66 /* fixed point math idea taken from Imlib by
67 * Carsten Haitzler (Rasterman) */
68 dx = (image->width << 16) / new_width;
69 dy = (image->height << 16) / new_height;
70
71 py = 0;
72
73 d = img->data;
74
75 if (image->format == RRGBAFormat) {
76 for (y = 0; y < new_height; y++) {
77 t = image->width * (py >> 16);
78
79 s = image->data + (t << 2); /* image->data+t*4 */
80
81 ox = 0;
82 px = 0;
83 for (x = 0; x < new_width; x++) {
84 px += dx;
85
86 *(d++) = *(s);
87 *(d++) = *(s + 1);
88 *(d++) = *(s + 2);
89 *(d++) = *(s + 3);
90
91 t = (px - ox) >> 16;
92 ox += t << 16;
93
94 s += t << 2; /* t*4 */
95 }
96 py += dy;
97 }
98 } else {
99 for (y = 0; y < new_height; y++) {
100 t = image->width * (py >> 16);
101
102 s = image->data + (t << 1) + t; /* image->data+t*3 */
103
104 ox = 0;
105 px = 0;
106 for (x = 0; x < new_width; x++) {
107 px += dx;
108
109 *(d++) = *(s);
110 *(d++) = *(s + 1);
111 *(d++) = *(s + 2);
112
113 t = (px - ox) >> 16;
114 ox += t << 16;
115
116 s += (t << 1) + t; /* t*3 */
117 }
118 py += dy;
119 }
120 }
121
122 return img;
123 }
124
125 /*
126 * Filtered Image Rescaling code copy/pasted from
127 * Graphics Gems III
128 * Public Domain 1991 by Dale Schumacher
129 */
130
131 /*
132 * filter function definitions
133 */
134 #define box_support (0.5)
135
box_filter(double t)136 static double box_filter(double t)
137 {
138 if ((t > -0.5) && (t <= 0.5))
139 return (1.0);
140 return (0.0);
141 }
142
143 #define triangle_support (1.0)
144
triangle_filter(double t)145 static double triangle_filter(double t)
146 {
147 if (t < 0.0)
148 t = -t;
149 if (t < 1.0)
150 return (1.0 - t);
151 return (0.0);
152 }
153
154 #define bell_support (1.5)
155
bell_filter(double t)156 static double bell_filter(double t) /* box (*) box (*) box */
157 {
158 if (t < 0)
159 t = -t;
160 if (t < .5)
161 return (.75 - (t * t));
162 if (t < 1.5) {
163 t = (t - 1.5);
164 return (.5 * (t * t));
165 }
166 return (0.0);
167 }
168
169 #define B_spline_support (2.0)
170
B_spline_filter(double t)171 static double B_spline_filter(double t) /* box (*) box (*) box (*) box */
172 {
173 double tt;
174
175 if (t < 0)
176 t = -t;
177 if (t < 1) {
178 tt = t * t;
179 return ((.5 * tt * t) - tt + (2.0 / 3.0));
180 } else if (t < 2) {
181 t = 2 - t;
182 return ((1.0 / 6.0) * (t * t * t));
183 }
184 return (0.0);
185 }
186
sinc(double x)187 static double sinc(double x)
188 {
189 /*
190 * The original code did this:
191 * if (x != 0) ...
192 * This code is unsafe, it should be:
193 * if (fabs(x) > EPSILON) ...
194 *
195 * But the call to fabs is already done in the *ONLY* function
196 * that call sinc: 'Lanczos3_filter'
197 *
198 * The goal was to avoid a Divide-by-0 error, now we also
199 * avoid a +/-inf result too
200 */
201 x *= WM_PI;
202 if (x > 1.0E-9)
203 return (sin(x) / x);
204 return (1.0);
205 }
206
207 #define Lanczos3_support (3.0)
208
Lanczos3_filter(double t)209 static double Lanczos3_filter(double t)
210 {
211 if (t < 0)
212 t = -t;
213 if (t < 3.0)
214 return (sinc(t) * sinc(t / 3.0));
215 return (0.0);
216 }
217
218 #define Mitchell_support (2.0)
219
220 #define B (1.0 / 3.0)
221 #define C (1.0 / 3.0)
222
Mitchell_filter(double t)223 static double Mitchell_filter(double t)
224 {
225 double tt;
226
227 tt = t * t;
228 if (t < 0)
229 t = -t;
230 if (t < 1.0) {
231 t = (((12.0 - 9.0 * B - 6.0 * C) * (t * tt))
232 + ((-18.0 + 12.0 * B + 6.0 * C) * tt)
233 + (6.0 - 2 * B));
234 return (t / 6.0);
235 } else if (t < 2.0) {
236 t = (((-1.0 * B - 6.0 * C) * (t * tt))
237 + ((6.0 * B + 30.0 * C) * tt)
238 + ((-12.0 * B - 48.0 * C) * t)
239 + (8.0 * B + 24 * C));
240 return (t / 6.0);
241 }
242 return (0.0);
243 }
244
245 static double (*filterf)(double) = Mitchell_filter;
246 static double fwidth = Mitchell_support;
247
wraster_change_filter(RScalingFilter type)248 void wraster_change_filter(RScalingFilter type)
249 {
250 switch (type) {
251 case RBoxFilter:
252 filterf = box_filter;
253 fwidth = box_support;
254 break;
255 case RTriangleFilter:
256 filterf = triangle_filter;
257 fwidth = triangle_support;
258 break;
259 case RBellFilter:
260 filterf = bell_filter;
261 fwidth = bell_support;
262 break;
263 case RBSplineFilter:
264 filterf = B_spline_filter;
265 fwidth = B_spline_support;
266 break;
267 case RLanczos3Filter:
268 filterf = Lanczos3_filter;
269 fwidth = Lanczos3_support;
270 break;
271 default:
272 case RMitchellFilter:
273 filterf = Mitchell_filter;
274 fwidth = Mitchell_support;
275 break;
276 }
277 }
278
279 /*
280 * image rescaling routine
281 */
282
283 typedef struct {
284 int pixel;
285 double weight;
286 } CONTRIB;
287
288 typedef struct {
289 int n; /* number of contributors */
290 CONTRIB *p; /* pointer to list of contributions */
291 } CLIST;
292
293 /* clamp the input to the specified range */
294 #define CLAMP(v,l,h) ((v)<(l) ? (l) : (v) > (h) ? (h) : v)
295
296 /* return of calloc is not checked if NULL in the function below! */
RSmoothScaleImage(RImage * src,unsigned new_width,unsigned new_height)297 RImage *RSmoothScaleImage(RImage * src, unsigned new_width, unsigned new_height)
298 {
299 CLIST *contrib; /* array of contribution lists */
300 RImage *tmp; /* intermediate image */
301 double xscale, yscale; /* zoom scale factors */
302 int i, j, k; /* loop variables */
303 int n; /* pixel number */
304 double center, left, right; /* filter calculation variables */
305 double width, fscale; /* filter calculation variables */
306 double rweight, gweight, bweight;
307 RImage *dst;
308 unsigned char *p;
309 unsigned char *sp;
310 int sch = src->format == RRGBAFormat ? 4 : 3;
311
312 dst = RCreateImage(new_width, new_height, False);
313
314 /* create intermediate image to hold horizontal zoom */
315 tmp = RCreateImage(dst->width, src->height, False);
316 xscale = (double)new_width / (double)src->width;
317 yscale = (double)new_height / (double)src->height;
318
319 /* pre-calculate filter contributions for a row */
320 contrib = (CLIST *) calloc(new_width, sizeof(CLIST));
321 if (xscale < 1.0) {
322 width = fwidth / xscale;
323 fscale = 1.0 / xscale;
324 for (i = 0; i < new_width; ++i) {
325 contrib[i].n = 0;
326 contrib[i].p = (CONTRIB *) calloc((int) ceil(width * 2 + 1), sizeof(CONTRIB));
327 center = (double)i / xscale;
328 left = ceil(center - width);
329 right = floor(center + width);
330 for (j = left; j <= right; ++j) {
331 rweight = center - (double)j;
332 rweight = (*filterf) (rweight / fscale) / fscale;
333 if (j < 0) {
334 n = -j;
335 } else if (j >= src->width) {
336 n = (src->width - j) + src->width - 1;
337 } else {
338 n = j;
339 }
340 k = contrib[i].n++;
341 contrib[i].p[k].pixel = n * sch;
342 contrib[i].p[k].weight = rweight;
343 }
344 }
345 } else {
346
347 for (i = 0; i < new_width; ++i) {
348 contrib[i].n = 0;
349 contrib[i].p = (CONTRIB *) calloc((int) ceil(fwidth * 2 + 1), sizeof(CONTRIB));
350 center = (double)i / xscale;
351 left = ceil(center - fwidth);
352 right = floor(center + fwidth);
353 for (j = left; j <= right; ++j) {
354 rweight = center - (double)j;
355 rweight = (*filterf) (rweight);
356 if (j < 0) {
357 n = -j;
358 } else if (j >= src->width) {
359 n = (src->width - j) + src->width - 1;
360 } else {
361 n = j;
362 }
363 k = contrib[i].n++;
364 contrib[i].p[k].pixel = n * sch;
365 contrib[i].p[k].weight = rweight;
366 }
367 }
368 }
369
370 /* apply filter to zoom horizontally from src to tmp */
371 p = tmp->data;
372
373 for (k = 0; k < tmp->height; ++k) {
374 CONTRIB *pp;
375
376 sp = src->data + src->width * k * sch;
377
378 for (i = 0; i < tmp->width; ++i) {
379 rweight = gweight = bweight = 0.0;
380
381 pp = contrib[i].p;
382
383 for (j = 0; j < contrib[i].n; ++j) {
384 rweight += sp[pp[j].pixel] * pp[j].weight;
385 gweight += sp[pp[j].pixel + 1] * pp[j].weight;
386 bweight += sp[pp[j].pixel + 2] * pp[j].weight;
387 }
388 *p++ = CLAMP(rweight, 0, 255);
389 *p++ = CLAMP(gweight, 0, 255);
390 *p++ = CLAMP(bweight, 0, 255);
391 }
392 }
393
394 /* free the memory allocated for horizontal filter weights */
395 for (i = 0; i < new_width; ++i) {
396 free(contrib[i].p);
397 }
398 free(contrib);
399
400 /* pre-calculate filter contributions for a column */
401 contrib = (CLIST *) calloc(dst->height, sizeof(CLIST));
402 if (yscale < 1.0) {
403 width = fwidth / yscale;
404 fscale = 1.0 / yscale;
405 for (i = 0; i < dst->height; ++i) {
406 contrib[i].n = 0;
407 contrib[i].p = (CONTRIB *) calloc((int) ceil(width * 2 + 1), sizeof(CONTRIB));
408 center = (double)i / yscale;
409 left = ceil(center - width);
410 right = floor(center + width);
411 for (j = left; j <= right; ++j) {
412 rweight = center - (double)j;
413 rweight = (*filterf) (rweight / fscale) / fscale;
414 if (j < 0) {
415 n = -j;
416 } else if (j >= tmp->height) {
417 n = (tmp->height - j) + tmp->height - 1;
418 } else {
419 n = j;
420 }
421 k = contrib[i].n++;
422 contrib[i].p[k].pixel = n * 3;
423 contrib[i].p[k].weight = rweight;
424 }
425 }
426 } else {
427 for (i = 0; i < dst->height; ++i) {
428 contrib[i].n = 0;
429 contrib[i].p = (CONTRIB *) calloc((int) ceil(fwidth * 2 + 1), sizeof(CONTRIB));
430 center = (double)i / yscale;
431 left = ceil(center - fwidth);
432 right = floor(center + fwidth);
433 for (j = left; j <= right; ++j) {
434 rweight = center - (double)j;
435 rweight = (*filterf) (rweight);
436 if (j < 0) {
437 n = -j;
438 } else if (j >= tmp->height) {
439 n = (tmp->height - j) + tmp->height - 1;
440 } else {
441 n = j;
442 }
443 k = contrib[i].n++;
444 contrib[i].p[k].pixel = n * 3;
445 contrib[i].p[k].weight = rweight;
446 }
447 }
448 }
449
450 /* apply filter to zoom vertically from tmp to dst */
451 sp = malloc(tmp->height * 3);
452
453 for (k = 0; k < new_width; ++k) {
454 CONTRIB *pp;
455
456 p = dst->data + k * 3;
457
458 /* copy a column into a row */
459 {
460 int i;
461 unsigned char *p, *d;
462
463 d = sp;
464 for (i = tmp->height, p = tmp->data + k * 3; i-- > 0; p += tmp->width * 3) {
465 *d++ = *p;
466 *d++ = *(p + 1);
467 *d++ = *(p + 2);
468 }
469 }
470 for (i = 0; i < new_height; ++i) {
471 rweight = gweight = bweight = 0.0;
472
473 pp = contrib[i].p;
474
475 for (j = 0; j < contrib[i].n; ++j) {
476 rweight += sp[pp[j].pixel] * pp[j].weight;
477 gweight += sp[pp[j].pixel + 1] * pp[j].weight;
478 bweight += sp[pp[j].pixel + 2] * pp[j].weight;
479 }
480 *p = CLAMP(rweight, 0, 255);
481 *(p + 1) = CLAMP(gweight, 0, 255);
482 *(p + 2) = CLAMP(bweight, 0, 255);
483 p += new_width * 3;
484 }
485 }
486 free(sp);
487
488 /* free the memory allocated for vertical filter weights */
489 for (i = 0; i < dst->height; ++i) {
490 free(contrib[i].p);
491 }
492 free(contrib);
493
494 RReleaseImage(tmp);
495
496 return dst;
497 }
498