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(double t)
 238 {
 239         /* f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1 */
 240         if (t < 0.0)
 241                 t = -t;
 242         if (t < 1.0)
 243                 return ((2.0 * t - 3.0) * t * t + 1.0);
 244         return (0.0);
 245 }
 246 #endif
 247 #define box_support             (0.5)
 248
 249 static double box_filter(double t)
 250 {
 251         if ((t > -0.5) && (t <= 0.5))
 252                 return (1.0);
 253         return (0.0);
 254 }
 255
 256 #define triangle_support        (1.0)
 257
 258 static double triangle_filter(double t)
 259 {
 260         if (t < 0.0)
 261                 t = -t;
 262         if (t < 1.0)
 263                 return (1.0 - t);
 264         return (0.0);
 265 }
 266
 267 #define bell_support            (1.5)
 268
 269 static double bell_filter(double t)  /* box (*) box (*) box */
 270 {
 271         if (t < 0)
 272                 t = -t;
 273         if (t < .5)
 274                 return (.75 - (t * t));
 275         if (t < 1.5) {
 276                 t = (t - 1.5);
 277                 return (.5 * (t * t));
 278         }
 279         return (0.0);
 280 }
 281
 282 #define B_spline_support        (2.0)
 283
 284 static double B_spline_filter(double t)   /* box (*) box (*) box (*) box */
 285 {
 286         double tt;
 287
 288         if (t < 0)
 289                 t = -t;
 290         if (t < 1) {
 291                 tt = t * t;
 292                 return ((.5 * tt * t) - tt + (2.0 / 3.0));
 293         } else if (t < 2) {
 294                 t = 2 - t;
 295                 return ((1.0 / 6.0) * (t * t * t));
 296         }
 297         return (0.0);
 298 }
 299
 300 static double sinc(double x)
 301 {
 302         x *= PI;
 303         if (x != 0)
 304                 return (sin(x) / x);
 305         return (1.0);
 306 }
 307
 308 #define Lanczos3_support        (3.0)
 309
 310 static double Lanczos3_filter(double t)
 311 {
 312         if (t < 0)
 313                 t = -t;
 314         if (t < 3.0)
 315                 return (sinc(t) * sinc(t / 3.0));
 316         return (0.0);
 317 }
 318
 319 #define Mitchell_support        (2.0)
 320
 321 #define B       (1.0 / 3.0)
 322 #define C       (1.0 / 3.0)
 323
 324 static double Mitchell_filter(double t)
 325 {
 326         double tt;
 327
 328         tt = t * t;
 329         if (t < 0)
 330                 t = -t;
 331         if (t < 1.0) {
 332                 t = (((12.0 - 9.0 * B - 6.0 * C) * (t * tt))
 333                      + ((-18.0 + 12.0 * B + 6.0 * C) * tt)
 334                      + (6.0 - 2 * B));
 335                 return (t / 6.0);
 336         } else if (t < 2.0) {
 337                 t = (((-1.0 * B - 6.0 * C) * (t * tt))
 338                      + ((6.0 * B + 30.0 * C) * tt)
 339                      + ((-12.0 * B - 48.0 * C) * t)
 340                      + (8.0 * B + 24 * C));
 341                 return (t / 6.0);
 342         }
 343         return (0.0);
 344 }
 345
 346 static double (*filterf)(double) = Mitchell_filter;
 347 static double fwidth = Mitchell_support;
 348
 349 void _wraster_change_filter(int type)
 350 {
 351         switch (type) {
 352         case RBoxFilter:
 353                 filterf = box_filter;
 354                 fwidth = box_support;
 355                 break;
 356         case RTriangleFilter:
 357                 filterf = triangle_filter;
 358                 fwidth = triangle_support;
 359                 break;
 360         case RBellFilter:
 361                 filterf = bell_filter;
 362                 fwidth = bell_support;
 363                 break;
 364         case RBSplineFilter:
 365                 filterf = B_spline_filter;
 366                 fwidth = B_spline_support;
 367                 break;
 368         case RLanczos3Filter:
 369                 filterf = Lanczos3_filter;
 370                 fwidth = Lanczos3_support;
 371                 break;
 372         default:
 373         case RMitchellFilter:
 374                 filterf = Mitchell_filter;
 375                 fwidth = Mitchell_support;
 376                 break;
 377         }
 378 }
 379
 380 /*
 381  *      image rescaling routine
 382  */
 383
 384 typedef struct {
 385         int pixel;
 386         double weight;
 387 } CONTRIB;
 388
 389 typedef struct {
 390         int n;                  /* number of contributors */
 391         CONTRIB *p;             /* pointer to list of contributions */
 392 } CLIST;
 393
 394 CLIST *contrib;                 /* array of contribution lists */
 395
 396 /* clamp the input to the specified range */
 397 #define CLAMP(v,l,h)    ((v)<(l) ? (l) : (v) > (h) ? (h) : v)
 398
 399 /* return of calloc is not checked if NULL in the function below! */
 400 RImage *RSmoothScaleImage(RImage * src, unsigned new_width, unsigned new_height)
 401 {
 402         RImage *tmp;            /* intermediate image */
 403         double xscale, yscale;  /* zoom scale factors */
 404         int i, j, k;            /* loop variables */
 405         int n;                  /* pixel number */
 406         double center, left, right;     /* filter calculation variables */
 407         double width, fscale;   /* filter calculation variables */
 408         double rweight, gweight, bweight;
 409         RImage *dst;
 410         unsigned char *p;
 411         unsigned char *sp;
 412         int sch = src->format == RRGBAFormat ? 4 : 3;
 413
 414         dst = RCreateImage(new_width, new_height, False);
 415
 416         /* create intermediate image to hold horizontal zoom */
 417         tmp = RCreateImage(dst->width, src->height, False);
 418         xscale = (double)new_width / (double)src->width;
 419         yscale = (double)new_height / (double)src->height;
 420
 421         /* pre-calculate filter contributions for a row */
 422         contrib = (CLIST *) calloc(new_width, sizeof(CLIST));
 423         if (xscale < 1.0) {
 424                 width = fwidth / xscale;
 425                 fscale = 1.0 / xscale;
 426                 for (i = 0; i < new_width; ++i) {
 427                         contrib[i].n = 0;
 428                         contrib[i].p = (CONTRIB *) calloc((int)(width * 2 + 1), sizeof(CONTRIB));
 429                         center = (double)i / xscale;
 430                         left = ceil(center - width);
 431                         right = floor(center + width);
 432                         for (j = left; j <= right; ++j) {
 433                                 rweight = center - (double)j;
 434                                 rweight = (*filterf) (rweight / fscale) / fscale;
 435                                 if (j < 0) {
 436                                         n = -j;
 437                                 } else if (j >= src->width) {
 438                                         n = (src->width - j) + src->width - 1;
 439                                 } else {
 440                                         n = j;
 441                                 }
 442                                 k = contrib[i].n++;
 443                                 contrib[i].p[k].pixel = n * sch;
 444                                 contrib[i].p[k].weight = rweight;
 445                         }
 446                 }
 447         } else {
 448
 449                 for (i = 0; i < new_width; ++i) {
 450                         contrib[i].n = 0;
 451                         contrib[i].p = (CONTRIB *) calloc((int)(fwidth * 2 + 1), sizeof(CONTRIB));
 452                         center = (double)i / xscale;
 453                         left = ceil(center - fwidth);
 454                         right = floor(center + fwidth);
 455                         for (j = left; j <= right; ++j) {
 456                                 rweight = center - (double)j;
 457                                 rweight = (*filterf) (rweight);
 458                                 if (j < 0) {
 459                                         n = -j;
 460                                 } else if (j >= src->width) {
 461                                         n = (src->width - j) + src->width - 1;
 462                                 } else {
 463                                         n = j;
 464                                 }
 465                                 k = contrib[i].n++;
 466                                 contrib[i].p[k].pixel = n * sch;
 467                                 contrib[i].p[k].weight = rweight;
 468                         }
 469                 }
 470         }
 471
 472         /* apply filter to zoom horizontally from src to tmp */
 473         p = tmp->data;
 474
 475         for (k = 0; k < tmp->height; ++k) {
 476                 CONTRIB *pp;
 477
 478                 sp = src->data + src->width * k * sch;
 479
 480                 for (i = 0; i < tmp->width; ++i) {
 481                         rweight = gweight = bweight = 0.0;
 482
 483                         pp = contrib[i].p;
 484
 485                         for (j = 0; j < contrib[i].n; ++j) {
 486                                 rweight += sp[pp[j].pixel] * pp[j].weight;
 487                                 gweight += sp[pp[j].pixel + 1] * pp[j].weight;
 488                                 bweight += sp[pp[j].pixel + 2] * pp[j].weight;
 489                         }
 490                         *p++ = CLAMP(rweight, 0, 255);
 491                         *p++ = CLAMP(gweight, 0, 255);
 492                         *p++ = CLAMP(bweight, 0, 255);
 493                 }
 494         }
 495
 496         /* free the memory allocated for horizontal filter weights */
 497         for (i = 0; i < tmp->width; ++i) {
 498                 free(contrib[i].p);
 499         }
 500         free(contrib);
 501
 502         /* pre-calculate filter contributions for a column */
 503         contrib = (CLIST *) calloc(dst->height, sizeof(CLIST));
 504         if (yscale < 1.0) {
 505                 width = fwidth / yscale;
 506                 fscale = 1.0 / yscale;
 507                 for (i = 0; i < dst->height; ++i) {
 508                         contrib[i].n = 0;
 509                         contrib[i].p = (CONTRIB *) calloc((int)(width * 2 + 1), sizeof(CONTRIB));
 510                         center = (double)i / yscale;
 511                         left = ceil(center - width);
 512                         right = floor(center + width);
 513                         for (j = left; j <= right; ++j) {
 514                                 rweight = center - (double)j;
 515                                 rweight = (*filterf) (rweight / fscale) / fscale;
 516                                 if (j < 0) {
 517                                         n = -j;
 518                                 } else if (j >= tmp->height) {
 519                                         n = (tmp->height - j) + tmp->height - 1;
 520                                 } else {
 521                                         n = j;
 522                                 }
 523                                 k = contrib[i].n++;
 524                                 contrib[i].p[k].pixel = n * 3;
 525                                 contrib[i].p[k].weight = rweight;
 526                         }
 527                 }
 528         } else {
 529                 for (i = 0; i < dst->height; ++i) {
 530                         contrib[i].n = 0;
 531                         contrib[i].p = (CONTRIB *) calloc((int)(fwidth * 2 + 1), sizeof(CONTRIB));
 532                         center = (double)i / yscale;
 533                         left = ceil(center - fwidth);
 534                         right = floor(center + fwidth);
 535                         for (j = left; j <= right; ++j) {
 536                                 rweight = center - (double)j;
 537                                 rweight = (*filterf) (rweight);
 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         }
 551
 552         /* apply filter to zoom vertically from tmp to dst */
 553         sp = malloc(tmp->height * 3);
 554
 555         for (k = 0; k < new_width; ++k) {
 556                 CONTRIB *pp;
 557
 558                 p = dst->data + k * 3;
 559
 560                 /* copy a column into a row */
 561                 {
 562                         int i;
 563                         unsigned char *p, *d;
 564
 565                         d = sp;
 566                         for (i = tmp->height, p = tmp->data + k * 3; i-- > 0; p += tmp->width * 3) {
 567                                 *d++ = *p;
 568                                 *d++ = *(p + 1);
 569                                 *d++ = *(p + 2);
 570                         }
 571                 }
 572                 for (i = 0; i < new_height; ++i) {
 573                         rweight = gweight = bweight = 0.0;
 574
 575                         pp = contrib[i].p;
 576
 577                         for (j = 0; j < contrib[i].n; ++j) {
 578                                 rweight += sp[pp[j].pixel] * pp[j].weight;
 579                                 gweight += sp[pp[j].pixel + 1] * pp[j].weight;
 580                                 bweight += sp[pp[j].pixel + 2] * pp[j].weight;
 581                         }
 582                         *p = CLAMP(rweight, 0, 255);
 583                         *(p + 1) = CLAMP(gweight, 0, 255);
 584                         *(p + 2) = CLAMP(bweight, 0, 255);
 585                         p += new_width * 3;
 586                 }
 587         }
 588         free(sp);
 589
 590         /* free the memory allocated for vertical filter weights */
 591         for (i = 0; i < dst->height; ++i) {
 592                 free(contrib[i].p);
 593         }
 594         free(contrib);
 595
 596         RReleaseImage(tmp);
 597
 598         return dst;
 599 }