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., 51 Franklin St, Fifth Floor, Boston,
  20  *  MA 02110-1301, USA.
  21  */
  22
  23 #include <config.h>
  24
  25 #include <stdlib.h>
  26 #include <stdio.h>
  27 #include <string.h>
  28 #include <X11/Xlib.h>
  29 #include <math.h>
  30 #include <assert.h>
  31
  32 #include "wraster.h"
  33 #include "scale.h"
  34 #include "wr_i18n.h"
  35
  36
  37 /*
  38  *----------------------------------------------------------------------
  39  * RScaleImage--
  40  *      Creates a scaled copy of an image.
  41  *
  42  * Returns:
  43  *      The new scaled image.
  44  *
  45  *----------------------------------------------------------------------
  46  */
  47 RImage *RScaleImage(RImage * image, unsigned new_width, unsigned new_height)
  48 {
  49         int ox;
  50         int px, py;
  51         register int x, y, t;
  52         int dx, dy;
  53         unsigned char *s;
  54         unsigned char *d;
  55         RImage *img;
  56
  57         if (image == NULL)
  58                 return NULL;
  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         /*
 192          * The original code did this:
 193          *   if (x != 0) ...
 194          * This code is unsafe, it should be:
 195          *   if (fabs(x) > EPSILON) ...
 196          *
 197          * But the call to fabs is already done in the *ONLY* function
 198          * that call sinc: 'Lanczos3_filter'
 199          *
 200          * The goal was to avoid a Divide-by-0 error, now we also
 201          * avoid a +/-inf result too
 202          */
 203         x *= WM_PI;
 204         if (x > 1.0E-9)
 205                 return (sin(x) / x);
 206         return (1.0);
 207 }
 208
 209 #define Lanczos3_support        (3.0)
 210
 211 static double Lanczos3_filter(double t)
 212 {
 213         if (t < 0)
 214                 t = -t;
 215         if (t < 3.0)
 216                 return (sinc(t) * sinc(t / 3.0));
 217         return (0.0);
 218 }
 219
 220 #define Mitchell_support        (2.0)
 221
 222 #define B       (1.0 / 3.0)
 223 #define C       (1.0 / 3.0)
 224
 225 static double Mitchell_filter(double t)
 226 {
 227         double tt;
 228
 229         tt = t * t;
 230         if (t < 0)
 231                 t = -t;
 232         if (t < 1.0) {
 233                 t = (((12.0 - 9.0 * B - 6.0 * C) * (t * tt))
 234                      + ((-18.0 + 12.0 * B + 6.0 * C) * tt)
 235                      + (6.0 - 2 * B));
 236                 return (t / 6.0);
 237         } else if (t < 2.0) {
 238                 t = (((-1.0 * B - 6.0 * C) * (t * tt))
 239                      + ((6.0 * B + 30.0 * C) * tt)
 240                      + ((-12.0 * B - 48.0 * C) * t)
 241                      + (8.0 * B + 24 * C));
 242                 return (t / 6.0);
 243         }
 244         return (0.0);
 245 }
 246
 247 static double (*filterf)(double) = Mitchell_filter;
 248 static double fwidth = Mitchell_support;
 249
 250 void wraster_change_filter(RScalingFilter type)
 251 {
 252         switch (type) {
 253         case RBoxFilter:
 254                 filterf = box_filter;
 255                 fwidth = box_support;
 256                 break;
 257         case RTriangleFilter:
 258                 filterf = triangle_filter;
 259                 fwidth = triangle_support;
 260                 break;
 261         case RBellFilter:
 262                 filterf = bell_filter;
 263                 fwidth = bell_support;
 264                 break;
 265         case RBSplineFilter:
 266                 filterf = B_spline_filter;
 267                 fwidth = B_spline_support;
 268                 break;
 269         case RLanczos3Filter:
 270                 filterf = Lanczos3_filter;
 271                 fwidth = Lanczos3_support;
 272                 break;
 273         default:
 274         case RMitchellFilter:
 275                 filterf = Mitchell_filter;
 276                 fwidth = Mitchell_support;
 277                 break;
 278         }
 279 }
 280
 281 /*
 282  *      image rescaling routine
 283  */
 284
 285 typedef struct {
 286         int pixel;
 287         double weight;
 288 } CONTRIB;
 289
 290 typedef struct {
 291         int n;                  /* number of contributors */
 292         CONTRIB *p;             /* pointer to list of contributions */
 293 } CLIST;
 294
 295 /* clamp the input to the specified range */
 296 #define CLAMP(v,l,h)    ((v)<(l) ? (l) : (v) > (h) ? (h) : v)
 297
 298 /* return of calloc is not checked if NULL in the function below! */
 299 RImage *RSmoothScaleImage(RImage * src, unsigned new_width, unsigned new_height)
 300 {
 301         CLIST *contrib;                 /* array of contribution lists */
 302         RImage *tmp;            /* intermediate image */
 303         double xscale, yscale;  /* zoom scale factors */
 304         int i, j, k;            /* loop variables */
 305         int n;                  /* pixel number */
 306         double center, left, right;     /* filter calculation variables */
 307         double width, fscale;   /* filter calculation variables */
 308         double rweight, gweight, bweight;
 309         RImage *dst;
 310         unsigned char *p;
 311         unsigned char *sp;
 312         int sch = src->format == RRGBAFormat ? 4 : 3;
 313
 314         dst = RCreateImage(new_width, new_height, False);
 315
 316         /* create intermediate image to hold horizontal zoom */
 317         tmp = RCreateImage(dst->width, src->height, False);
 318         xscale = (double)new_width / (double)src->width;
 319         yscale = (double)new_height / (double)src->height;
 320
 321         /* pre-calculate filter contributions for a row */
 322         contrib = (CLIST *) calloc(new_width, sizeof(CLIST));
 323         if (xscale < 1.0) {
 324                 width = fwidth / xscale;
 325                 fscale = 1.0 / xscale;
 326                 for (i = 0; i < new_width; ++i) {
 327                         contrib[i].n = 0;
 328                         contrib[i].p = (CONTRIB *) calloc((int) ceil(width * 2 + 1), sizeof(CONTRIB));
 329                         center = (double)i / xscale;
 330                         left = ceil(center - width);
 331                         right = floor(center + width);
 332                         for (j = left; j <= right; ++j) {
 333                                 rweight = center - (double)j;
 334                                 rweight = (*filterf) (rweight / fscale) / fscale;
 335                                 if (j < 0) {
 336                                         n = -j;
 337                                 } else if (j >= src->width) {
 338                                         n = (src->width - j) + src->width - 1;
 339                                 } else {
 340                                         n = j;
 341                                 }
 342                                 k = contrib[i].n++;
 343                                 contrib[i].p[k].pixel = n * sch;
 344                                 contrib[i].p[k].weight = rweight;
 345                         }
 346                 }
 347         } else {
 348
 349                 for (i = 0; i < new_width; ++i) {
 350                         contrib[i].n = 0;
 351                         contrib[i].p = (CONTRIB *) calloc((int) ceil(fwidth * 2 + 1), sizeof(CONTRIB));
 352                         center = (double)i / xscale;
 353                         left = ceil(center - fwidth);
 354                         right = floor(center + fwidth);
 355                         for (j = left; j <= right; ++j) {
 356                                 rweight = center - (double)j;
 357                                 rweight = (*filterf) (rweight);
 358                                 if (j < 0) {
 359                                         n = -j;
 360                                 } else if (j >= src->width) {
 361                                         n = (src->width - j) + src->width - 1;
 362                                 } else {
 363                                         n = j;
 364                                 }
 365                                 k = contrib[i].n++;
 366                                 contrib[i].p[k].pixel = n * sch;
 367                                 contrib[i].p[k].weight = rweight;
 368                         }
 369                 }
 370         }
 371
 372         /* apply filter to zoom horizontally from src to tmp */
 373         p = tmp->data;
 374
 375         for (k = 0; k < tmp->height; ++k) {
 376                 CONTRIB *pp;
 377
 378                 sp = src->data + src->width * k * sch;
 379
 380                 for (i = 0; i < tmp->width; ++i) {
 381                         rweight = gweight = bweight = 0.0;
 382
 383                         pp = contrib[i].p;
 384
 385                         for (j = 0; j < contrib[i].n; ++j) {
 386                                 rweight += sp[pp[j].pixel] * pp[j].weight;
 387                                 gweight += sp[pp[j].pixel + 1] * pp[j].weight;
 388                                 bweight += sp[pp[j].pixel + 2] * pp[j].weight;
 389                         }
 390                         *p++ = CLAMP(rweight, 0, 255);
 391                         *p++ = CLAMP(gweight, 0, 255);
 392                         *p++ = CLAMP(bweight, 0, 255);
 393                 }
 394         }
 395
 396         /* free the memory allocated for horizontal filter weights */
 397         for (i = 0; i < new_width; ++i) {
 398                 free(contrib[i].p);
 399         }
 400         free(contrib);
 401
 402         /* pre-calculate filter contributions for a column */
 403         contrib = (CLIST *) calloc(dst->height, sizeof(CLIST));
 404         if (yscale < 1.0) {
 405                 width = fwidth / yscale;
 406                 fscale = 1.0 / yscale;
 407                 for (i = 0; i < dst->height; ++i) {
 408                         contrib[i].n = 0;
 409                         contrib[i].p = (CONTRIB *) calloc((int) ceil(width * 2 + 1), sizeof(CONTRIB));
 410                         center = (double)i / yscale;
 411                         left = ceil(center - width);
 412                         right = floor(center + width);
 413                         for (j = left; j <= right; ++j) {
 414                                 rweight = center - (double)j;
 415                                 rweight = (*filterf) (rweight / fscale) / fscale;
 416                                 if (j < 0) {
 417                                         n = -j;
 418                                 } else if (j >= tmp->height) {
 419                                         n = (tmp->height - j) + tmp->height - 1;
 420                                 } else {
 421                                         n = j;
 422                                 }
 423                                 k = contrib[i].n++;
 424                                 contrib[i].p[k].pixel = n * 3;
 425                                 contrib[i].p[k].weight = rweight;
 426                         }
 427                 }
 428         } else {
 429                 for (i = 0; i < dst->height; ++i) {
 430                         contrib[i].n = 0;
 431                         contrib[i].p = (CONTRIB *) calloc((int) ceil(fwidth * 2 + 1), sizeof(CONTRIB));
 432                         center = (double)i / yscale;
 433                         left = ceil(center - fwidth);
 434                         right = floor(center + fwidth);
 435                         for (j = left; j <= right; ++j) {
 436                                 rweight = center - (double)j;
 437                                 rweight = (*filterf) (rweight);
 438                                 if (j < 0) {
 439                                         n = -j;
 440                                 } else if (j >= tmp->height) {
 441                                         n = (tmp->height - j) + tmp->height - 1;
 442                                 } else {
 443                                         n = j;
 444                                 }
 445                                 k = contrib[i].n++;
 446                                 contrib[i].p[k].pixel = n * 3;
 447                                 contrib[i].p[k].weight = rweight;
 448                         }
 449                 }
 450         }
 451
 452         /* apply filter to zoom vertically from tmp to dst */
 453         sp = malloc(tmp->height * 3);
 454
 455         for (k = 0; k < new_width; ++k) {
 456                 CONTRIB *pp;
 457
 458                 p = dst->data + k * 3;
 459
 460                 /* copy a column into a row */
 461                 {
 462                         int i;
 463                         unsigned char *p, *d;
 464
 465                         d = sp;
 466                         for (i = tmp->height, p = tmp->data + k * 3; i-- > 0; p += tmp->width * 3) {
 467                                 *d++ = *p;
 468                                 *d++ = *(p + 1);
 469                                 *d++ = *(p + 2);
 470                         }
 471                 }
 472                 for (i = 0; i < new_height; ++i) {
 473                         rweight = gweight = bweight = 0.0;
 474
 475                         pp = contrib[i].p;
 476
 477                         for (j = 0; j < contrib[i].n; ++j) {
 478                                 rweight += sp[pp[j].pixel] * pp[j].weight;
 479                                 gweight += sp[pp[j].pixel + 1] * pp[j].weight;
 480                                 bweight += sp[pp[j].pixel + 2] * pp[j].weight;
 481                         }
 482                         *p = CLAMP(rweight, 0, 255);
 483                         *(p + 1) = CLAMP(gweight, 0, 255);
 484                         *(p + 2) = CLAMP(bweight, 0, 255);
 485                         p += new_width * 3;
 486                 }
 487         }
 488         free(sp);
 489
 490         /* free the memory allocated for vertical filter weights */
 491         for (i = 0; i < dst->height; ++i) {
 492                 free(contrib[i].p);
 493         }
 494         free(contrib);
 495
 496         RReleaseImage(tmp);
 497
 498         return dst;
 499 }