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