wrlib/scale.c

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