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