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