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