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