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., 675 Mass Ave, Cambridge, MA 02139, USA.
  20  */
  21
  22 #include <config.h>
  23
  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  *----------------------------------------------------------------------
  42  * RScaleImage--
  43  *      Creates a scaled copy of an image.
  44  *
  45  * Returns:
  46  *      The new scaled image.
  47  *
  48  *----------------------------------------------------------------------
  49  */
  50 #ifndef broken_code
  51 RImage*
  52 RScaleImage(RImage *image, unsigned new_width, unsigned new_height)
  53 {
  54     int ox;
  55     int px, py;
  56     register int x, y, t;
  57     int dx, dy;
  58     unsigned char *s;
  59     unsigned char *d;
  60     RImage *img;
  61
  62     assert(new_width >= 0 && new_height >= 0);
  63
  64     if (new_width == image->width && new_height == image->height)
  65         return RCloneImage(image);
  66
  67     img = RCreateImage(new_width, new_height, image->format==RRGBAFormat);
  68
  69     if (!img)
  70         return NULL;
  71
  72     /* fixed point math idea taken from Imlib by
  73      * Carsten Haitzler (Rasterman) */
  74     dx = (image->width<<16)/new_width;
  75     dy = (image->height<<16)/new_height;
  76
  77     py = 0;
  78
  79     d = img->data;
  80
  81     if (image->format == RRGBAFormat) {
  82         for (y=0; y<new_height; y++) {
  83             t = image->width*(py>>16);
  84
  85             s = image->data+(t<<2); /* image->data+t*4 */
  86
  87             ox = 0;
  88             px = 0;
  89             for (x=0; x<new_width; x++) {
  90                 px += dx;
  91
  92                 *(d++) = *(s);
  93                 *(d++) = *(s+1);
  94                 *(d++) = *(s+2);
  95                 *(d++) = *(s+3);
  96
  97                 t = (px - ox)>>16;
  98                 ox += t<<16;
  99
 100                 s += t<<2; /* t*4 */
 101             }
 102             py += dy;
 103         }
 104     } else {
 105         for (y=0; y<new_height; y++) {
 106             t = image->width*(py>>16);
 107
 108             s = image->data+(t<<1)+t; /* image->data+t*3 */
 109
 110             ox = 0;
 111             px = 0;
 112             for (x=0; x<new_width; x++) {
 113                 px += dx;
 114
 115                 *(d++) = *(s);
 116                 *(d++) = *(s+1);
 117                 *(d++) = *(s+2);
 118
 119                 t = (px - ox)>>16;
 120                 ox += t<<16;
 121
 122                 s += (t<<1)+t; /* t*3 */
 123             }
 124             py += dy;
 125         }
 126     }
 127
 128     return img;
 129 }
 130
 131 #else
 132 RImage*
 133 RScaleImage(RImage *src, unsigned new_width, unsigned new_height)
 134 {
 135     int ddy, ee;
 136     int h2;
 137     int yd;
 138     int xd, xs;
 139     RImage *dst;
 140     int e, xd2;
 141     unsigned char *sr, *sg, *sb, *sa;
 142     unsigned char *dr, *dg, *db, *da;
 143     int ys = 0;
 144
 145
 146
 147     dst = RCreateImage(new_width, new_height, src->data[3]!=NULL);
 148
 149     ddy = src->height/2;
 150     ee = (ddy/2) - dst->height;
 151     h2 = new_height/2;
 152
 153     xd = dst->width;
 154     xs = src->width/2;
 155     e = (src->width/2)-xd;
 156     xd2 = xd/2;
 157
 158     sr = src->data[0];
 159     sg = src->data[1];
 160     sb = src->data[2];
 161     sa = src->data[3];
 162
 163     dr = dst->data[0];
 164     dg = dst->data[1];
 165     db = dst->data[2];
 166     da = dst->data[3];
 167
 168     if (sa == NULL) {
 169         for (yd = 0; yd < new_height; yd++) {
 170             int x;
 171
 172             sr = src->data[0] + ys * src->width;
 173             sg = src->data[1] + ys * src->width;
 174             sb = src->data[2] + ys * src->width;
 175
 176             for (x = 0; x < xd; x++) {
 177                 *(dr++) = *sr;
 178                 *(dg++) = *sg;
 179                 *(db++) = *sb;
 180
 181                 while (e >= 0) {
 182                     sr++;
 183                     sg++;
 184                     sb++;
 185                     e -= xd2;
 186                 }
 187                 e += xs;
 188             }
 189             while (ee >= 0) {
 190                 ys++;
 191                 ee -= h2;
 192             }
 193             ee += ddy;
 194         }
 195     } else {
 196         for (yd = 0; yd < new_height; yd++) {
 197             int x;
 198
 199             sr = src->data[0] + ys * src->width;
 200             sg = src->data[1] + ys * src->width;
 201             sb = src->data[2] + ys * src->width;
 202             sa = src->data[3] + ys * src->width;
 203
 204             for (x = 0; x < xd; x++) {
 205                 *(dr++) = *sr;
 206                 *(dg++) = *sg;
 207                 *(db++) = *sb;
 208                 *(da++) = *sa;
 209
 210                 while (e >= 0) {
 211                     sr++;
 212                     sg++;
 213                     sb++;
 214                     sa++;
 215                     e -= xd2;
 216                 }
 217                 e += xs;
 218             }
 219             while (ee >= 0) {
 220                 ys++;
 221                 ee -= h2;
 222             }
 223             ee += ddy;
 224         }
 225     }
 226
 227     return dst;
 228 }
 229 #endif
 230
 231
 232
 233 /*
 234  * Filtered Image Rescaling code copy/pasted from
 235  * Graphics Gems III
 236  * Public Domain 1991 by Dale Schumacher
 237  */
 238
 239
 240 /*
 241  *      filter function definitions
 242  */
 243 #if 0
 244 #define filter_support          (1.0)
 245
 246 static double
 247 filter(t)
 248 double t;
 249 {
 250     /* f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1 */
 251     if(t < 0.0) t = -t;
 252     if(t < 1.0) return((2.0 * t - 3.0) * t * t + 1.0);
 253     return(0.0);
 254 }
 255 #endif
 256 #define box_support             (0.5)
 257
 258 static double
 259 box_filter(t)
 260 double t;
 261 {
 262     if((t > -0.5) && (t <= 0.5)) return(1.0);
 263     return(0.0);
 264 }
 265
 266 #define triangle_support        (1.0)
 267
 268 static double
 269 triangle_filter(t)
 270 double t;
 271 {
 272     if(t < 0.0) t = -t;
 273     if(t < 1.0) return(1.0 - t);
 274     return(0.0);
 275 }
 276
 277 #define bell_support            (1.5)
 278
 279 static double
 280 bell_filter(t)          /* box (*) box (*) box */
 281 double t;
 282 {
 283     if(t < 0) t = -t;
 284     if(t < .5) return(.75 - (t * t));
 285     if(t < 1.5) {
 286         t = (t - 1.5);
 287         return(.5 * (t * t));
 288     }
 289     return(0.0);
 290 }
 291
 292 #define B_spline_support        (2.0)
 293
 294 static double
 295 B_spline_filter(t)      /* box (*) box (*) box (*) box */
 296 double t;
 297 {
 298     double tt;
 299
 300     if(t < 0) t = -t;
 301     if(t < 1) {
 302         tt = t * t;
 303         return((.5 * tt * t) - tt + (2.0 / 3.0));
 304     } else if(t < 2) {
 305         t = 2 - t;
 306         return((1.0 / 6.0) * (t * t * t));
 307     }
 308     return(0.0);
 309 }
 310
 311 static double
 312 sinc(x)
 313 double x;
 314 {
 315     x *= PI;
 316     if(x != 0) return(sin(x) / x);
 317     return(1.0);
 318 }
 319
 320 #define Lanczos3_support        (3.0)
 321
 322 static double
 323 Lanczos3_filter(t)
 324 double t;
 325 {
 326     if(t < 0) t = -t;
 327     if(t < 3.0) return(sinc(t) * sinc(t/3.0));
 328     return(0.0);
 329 }
 330
 331 #define Mitchell_support        (2.0)
 332
 333 #define B       (1.0 / 3.0)
 334 #define C       (1.0 / 3.0)
 335
 336 static double
 337 Mitchell_filter(t)
 338 double t;
 339 {
 340     double tt;
 341
 342     tt = t * t;
 343     if(t < 0) t = -t;
 344     if(t < 1.0) {
 345         t = (((12.0 - 9.0 * B - 6.0 * C) * (t * tt))
 346              + ((-18.0 + 12.0 * B + 6.0 * C) * tt)
 347              + (6.0 - 2 * B));
 348         return(t / 6.0);
 349     } else if(t < 2.0) {
 350         t = (((-1.0 * B - 6.0 * C) * (t * tt))
 351              + ((6.0 * B + 30.0 * C) * tt)
 352              + ((-12.0 * B - 48.0 * C) * t)
 353              + (8.0 * B + 24 * C));
 354         return(t / 6.0);
 355     }
 356     return(0.0);
 357 }
 358
 359 static double (*filterf)() = Mitchell_filter;
 360 static double fwidth = Mitchell_support;
 361
 362 void
 363 _wraster_change_filter(int type)
 364 {
 365     switch (type) {
 366     case RBoxFilter:
 367         filterf = box_filter;
 368         fwidth = box_support;
 369         break;
 370     case RTriangleFilter:
 371         filterf = triangle_filter;
 372         fwidth = triangle_support;
 373         break;
 374     case RBellFilter:
 375         filterf = bell_filter;
 376         fwidth = bell_support;
 377         break;
 378     case RBSplineFilter:
 379         filterf = B_spline_filter;
 380         fwidth = B_spline_support;
 381         break;
 382     case RLanczos3Filter:
 383         filterf = Lanczos3_filter;
 384         fwidth = Lanczos3_support;
 385         break;
 386     default:
 387     case RMitchellFilter:
 388         filterf = Mitchell_filter;
 389         fwidth = Mitchell_support;
 390         break;
 391     }
 392 }
 393
 394
 395 /*
 396  *      image rescaling routine
 397  */
 398
 399 typedef struct {
 400     int pixel;
 401     double      weight;
 402 } CONTRIB;
 403
 404 typedef struct {
 405     int n;              /* number of contributors */
 406     CONTRIB     *p;             /* pointer to list of contributions */
 407 } CLIST;
 408
 409 CLIST   *contrib;               /* array of contribution lists */
 410
 411 /* clamp the input to the specified range */
 412 #define CLAMP(v,l,h)    ((v)<(l) ? (l) : (v) > (h) ? (h) : v)
 413
 414
 415 /* return of calloc is not checked if NULL in the function below! */
 416 RImage*
 417 RSmoothScaleImage(RImage *src, unsigned new_width, unsigned new_height)
 418 {
 419     RImage *tmp;                       /* intermediate image */
 420     double xscale, yscale;             /* zoom scale factors */
 421     int i, j, k;                       /* loop variables */
 422     int n;                             /* pixel number */
 423     double center, left, right;        /* filter calculation variables */
 424     double width, fscale;              /* filter calculation variables */
 425     double rweight, gweight, bweight;
 426     RImage *dst;
 427     unsigned char *p;
 428     unsigned char *sp;
 429     int sch = src->format == RRGBAFormat ? 4 : 3;
 430
 431
 432     dst = RCreateImage(new_width, new_height, False);
 433
 434     /* create intermediate image to hold horizontal zoom */
 435     tmp = RCreateImage(dst->width, src->height, False);
 436     xscale = (double)new_width / (double)src->width;
 437     yscale = (double)new_height / (double)src->height;
 438
 439     /* pre-calculate filter contributions for a row */
 440     contrib = (CLIST *)calloc(new_width, sizeof(CLIST));
 441     if (xscale < 1.0) {
 442         width = fwidth / xscale;
 443         fscale = 1.0 / xscale;
 444         for (i = 0; i < new_width; ++i) {
 445             contrib[i].n = 0;
 446             contrib[i].p = (CONTRIB *)calloc((int)(width * 2 + 1),
 447                                              sizeof(CONTRIB));
 448             center = (double) i / xscale;
 449             left = ceil(center - width);
 450             right = floor(center + width);
 451             for(j = left; j <= right; ++j) {
 452                 rweight = center - (double) j;
 453                 rweight = (*filterf)(rweight / fscale) / fscale;
 454                 if(j < 0) {
 455                     n = -j;
 456                 } else if(j >= src->width) {
 457                     n = (src->width - j) + src->width - 1;
 458                 } else {
 459                     n = j;
 460                 }
 461                 k = contrib[i].n++;
 462                 contrib[i].p[k].pixel = n*sch;
 463                 contrib[i].p[k].weight = rweight;
 464             }
 465         }
 466     } else {
 467
 468         for(i = 0; i < new_width; ++i) {
 469             contrib[i].n = 0;
 470             contrib[i].p = (CONTRIB *)calloc((int) (fwidth * 2 + 1),
 471                                              sizeof(CONTRIB));
 472             center = (double) i / xscale;
 473             left = ceil(center - fwidth);
 474             right = floor(center + fwidth);
 475             for(j = left; j <= right; ++j) {
 476                 rweight = center - (double) j;
 477                 rweight = (*filterf)(rweight);
 478                 if(j < 0) {
 479                     n = -j;
 480                 } else if(j >= src->width) {
 481                     n = (src->width - j) + src->width - 1;
 482                 } else {
 483                     n = j;
 484                 }
 485                 k = contrib[i].n++;
 486                 contrib[i].p[k].pixel = n*sch;
 487                 contrib[i].p[k].weight = rweight;
 488             }
 489         }
 490     }
 491
 492     /* apply filter to zoom horizontally from src to tmp */
 493     p = tmp->data;
 494
 495
 496     for(k = 0; k < tmp->height; ++k) {
 497         CONTRIB *pp;
 498
 499         sp = src->data + src->width*k*sch;
 500
 501         for(i = 0; i < tmp->width; ++i) {
 502             rweight = gweight = bweight = 0.0;
 503
 504             pp = contrib[i].p;
 505
 506             for(j = 0; j < contrib[i].n; ++j) {
 507                 rweight += sp[pp[j].pixel] * pp[j].weight;
 508                 gweight += sp[pp[j].pixel+1] * pp[j].weight;
 509                 bweight += sp[pp[j].pixel+2] * pp[j].weight;
 510             }
 511             *p++ = CLAMP(rweight, 0, 255);
 512             *p++ = CLAMP(gweight, 0, 255);
 513             *p++ = CLAMP(bweight, 0, 255);
 514         }
 515     }
 516
 517     /* free the memory allocated for horizontal filter weights */
 518     for(i = 0; i < tmp->width; ++i) {
 519         free(contrib[i].p);
 520     }
 521     free(contrib);
 522
 523     /* pre-calculate filter contributions for a column */
 524     contrib = (CLIST *)calloc(dst->height, sizeof(CLIST));
 525     if(yscale < 1.0) {
 526         width = fwidth / yscale;
 527         fscale = 1.0 / yscale;
 528         for(i = 0; i < dst->height; ++i) {
 529             contrib[i].n = 0;
 530             contrib[i].p = (CONTRIB *)calloc((int) (width * 2 + 1),
 531                                              sizeof(CONTRIB));
 532             center = (double) i / yscale;
 533             left = ceil(center - width);
 534             right = floor(center + width);
 535             for(j = left; j <= right; ++j) {
 536                 rweight = center - (double) j;
 537                 rweight = (*filterf)(rweight / fscale) / fscale;
 538                 if(j < 0) {
 539                     n = -j;
 540                 } else if(j >= tmp->height) {
 541                     n = (tmp->height - j) + tmp->height - 1;
 542                 } else {
 543                     n = j;
 544                 }
 545                 k = contrib[i].n++;
 546                 contrib[i].p[k].pixel = n*3;
 547                 contrib[i].p[k].weight = rweight;
 548             }
 549         }
 550     } else {
 551         for(i = 0; i < dst->height; ++i) {
 552             contrib[i].n = 0;
 553             contrib[i].p = (CONTRIB *)calloc((int) (fwidth * 2 + 1),
 554                                              sizeof(CONTRIB));
 555             center = (double) i / yscale;
 556             left = ceil(center - fwidth);
 557             right = floor(center + fwidth);
 558             for(j = left; j <= right; ++j) {
 559                 rweight = center - (double) j;
 560                 rweight = (*filterf)(rweight);
 561                 if(j < 0) {
 562                     n = -j;
 563                 } else if(j >= tmp->height) {
 564                     n = (tmp->height - j) + tmp->height - 1;
 565                 } else {
 566                     n = j;
 567                 }
 568                 k = contrib[i].n++;
 569                 contrib[i].p[k].pixel = n*3;
 570                 contrib[i].p[k].weight = rweight;
 571             }
 572         }
 573     }
 574
 575     /* apply filter to zoom vertically from tmp to dst */
 576     sp = malloc(tmp->height*3);
 577
 578     for(k = 0; k < new_width; ++k) {
 579         CONTRIB *pp;
 580
 581         p = dst->data + k*3;
 582
 583         /* copy a column into a row */
 584         {
 585             int i;
 586             unsigned char *p, *d;
 587
 588             d = sp;
 589             for(i = tmp->height, p = tmp->data + k*3; i-- > 0;
 590                 p += tmp->width*3) {
 591                 *d++ = *p;
 592                 *d++ = *(p+1);
 593                 *d++ = *(p+2);
 594             }
 595         }
 596         for(i = 0; i < new_height; ++i) {
 597             rweight = gweight = bweight = 0.0;
 598
 599             pp = contrib[i].p;
 600
 601             for(j = 0; j < contrib[i].n; ++j) {
 602                 rweight += sp[pp[j].pixel] * pp[j].weight;
 603                 gweight += sp[pp[j].pixel+1] * pp[j].weight;
 604                 bweight += sp[pp[j].pixel+2] * pp[j].weight;
 605             }
 606             *p = CLAMP(rweight, 0, 255);
 607             *(p+1) = CLAMP(gweight, 0, 255);
 608             *(p+2) = CLAMP(bweight, 0, 255);
 609             p += new_width*3;
 610         }
 611     }
 612     free(sp);
 613
 614     /* free the memory allocated for vertical filter weights */
 615     for(i = 0; i < dst->height; ++i) {
 616         free(contrib[i].p);
 617     }
 618     free(contrib);
 619
 620     RReleaseImage(tmp);
 621
 622     return dst;
 623 }
 624