wrlib/scale.c

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