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