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