~ubuntu-branches/ubuntu/jaunty/gimp/jaunty-security

« back to all changes in this revision

Viewing changes to app/paint-funcs/scale-funcs.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel Holbach
  • Date: 2007-05-02 16:33:03 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20070502163303-bvzhjzbpw8qglc4y
Tags: 2.3.16-1ubuntu1
* Resynchronized with Debian, remaining Ubuntu changes:
  - debian/rules: i18n magic.
* debian/control.in:
  - Maintainer: Ubuntu Core Developers <ubuntu-devel@lists.ubuntu.com>
* debian/patches/02_help-message.patch,
  debian/patches/03_gimp.desktop.in.in.patch,
  debian/patches/10_dont_show_wizard.patch: updated.
* debian/patches/04_composite-signedness.patch,
  debian/patches/05_add-letter-spacing.patch: dropped, used upstream.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* GIMP - The GNU Image Manipulation Program
 
2
 * Copyright (C) 1995 Spencer Kimball and Peter Mattis
 
3
 *
 
4
 * This program is free software; you can redistribute it and/or modify
 
5
 * it under the terms of the GNU General Public License as published by
 
6
 * the Free Software Foundation; either version 2 of the License, or
 
7
 * (at your option) any later version.
 
8
 *
 
9
 * This program is distributed in the hope that it will be useful,
 
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
 * GNU General Public License for more details.
 
13
 *
 
14
 * You should have received a copy of the GNU General Public License
 
15
 * along with this program; if not, write to the Free Software
 
16
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 
17
 */
 
18
 
 
19
#include "config.h"
 
20
 
 
21
#include <string.h>
 
22
 
 
23
#include <glib-object.h>
 
24
 
 
25
#include "libgimpmath/gimpmath.h"
 
26
 
 
27
#include "paint-funcs-types.h"
 
28
 
 
29
#include "base/pixel-region.h"
 
30
 
 
31
#include "scale-funcs.h"
 
32
 
 
33
 
 
34
#define EPSILON          (0.0001)  /* arbitary small number for avoiding zero */
 
35
 
 
36
static void  scale_region_no_resample (PixelRegion           *srcPR,
 
37
                                       PixelRegion           *destPR);
 
38
static void  scale_region_lanczos     (PixelRegion           *srcPR,
 
39
                                       PixelRegion           *dstPR,
 
40
                                       GimpProgressFunc       progress_callback,
 
41
                                       gpointer               progress_data);
 
42
 
 
43
static void  expand_line              (gdouble               *dest,
 
44
                                       const gdouble         *src,
 
45
                                       gint                   bytes,
 
46
                                       gint                   old_width,
 
47
                                       gint                   width,
 
48
                                       GimpInterpolationType  interp);
 
49
static void  shrink_line              (gdouble               *dest,
 
50
                                       const gdouble         *src,
 
51
                                       gint                   bytes,
 
52
                                       gint                   old_width,
 
53
                                       gint                   width,
 
54
                                       GimpInterpolationType  interp);
 
55
 
 
56
 
 
57
/* Catmull-Rom spline - not bad
 
58
  * basic intro http://www.mvps.org/directx/articles/catmull/
 
59
  * This formula will calculate an interpolated point between pt1 and pt2
 
60
  * dx=0 returns pt1; dx=1 returns pt2
 
61
  */
 
62
 
 
63
static inline gdouble
 
64
cubic_spline_fit (gdouble dx,
 
65
                  gint    pt0,
 
66
                  gint    pt1,
 
67
                  gint    pt2,
 
68
                  gint    pt3)
 
69
{
 
70
 
 
71
  return (gdouble) ((( ( - pt0 + 3 * pt1 - 3 * pt2 + pt3 ) * dx +
 
72
      ( 2 * pt0 - 5 * pt1 + 4 * pt2 - pt3 ) ) * dx +
 
73
      ( - pt0 + pt2 ) ) * dx + (pt1 + pt1) ) / 2.0;
 
74
}
 
75
 
 
76
 
 
77
 
 
78
/*
 
79
 * non-interpolating scale_region.  [adam]
 
80
 */
 
81
static void
 
82
scale_region_no_resample (PixelRegion *srcPR,
 
83
                          PixelRegion *destPR)
 
84
{
 
85
  const gint  width       = destPR->w;
 
86
  const gint  height      = destPR->h;
 
87
  const gint  orig_width  = srcPR->w;
 
88
  const gint  orig_height = srcPR->h;
 
89
  const gint  bytes       = srcPR->bytes;
 
90
  gint       *x_src_offsets;
 
91
  gint       *y_src_offsets;
 
92
  gint       *offset;
 
93
  guchar     *src;
 
94
  guchar     *dest;
 
95
  gint        last_src_y;
 
96
  gint        row_bytes;
 
97
  gint        x, y;
 
98
  gint        b;
 
99
 
 
100
  /*  the data pointers...  */
 
101
  x_src_offsets = g_new (gint, width * bytes);
 
102
  y_src_offsets = g_new (gint, height);
 
103
  src  = g_new (guchar, orig_width * bytes);
 
104
  dest = g_new (guchar, width * bytes);
 
105
 
 
106
  /*  pre-calc the scale tables  */
 
107
  offset = x_src_offsets;
 
108
  for (x = 0; x < width; x++)
 
109
    for (b = 0; b < bytes; b++)
 
110
      *offset++ = b + bytes * ((x * orig_width + orig_width / 2) / width);
 
111
 
 
112
  offset = y_src_offsets;
 
113
  for (y = 0; y < height; y++)
 
114
    *offset++ = (y * orig_height + orig_height / 2) / height;
 
115
 
 
116
  /*  do the scaling  */
 
117
  row_bytes = width * bytes;
 
118
  last_src_y = -1;
 
119
 
 
120
  for (y = 0; y < height; y++)
 
121
    {
 
122
      /* if the source of this line was the same as the source
 
123
       *  of the last line, there's no point in re-rescaling.
 
124
       */
 
125
      if (y_src_offsets[y] != last_src_y)
 
126
        {
 
127
          pixel_region_get_row (srcPR, 0, y_src_offsets[y], orig_width, src, 1);
 
128
 
 
129
          for (x = 0; x < row_bytes ; x++)
 
130
            dest[x] = src[x_src_offsets[x]];
 
131
 
 
132
          last_src_y = y_src_offsets[y];
 
133
        }
 
134
 
 
135
      pixel_region_set_row (destPR, 0, y, width, dest);
 
136
    }
 
137
 
 
138
  g_free (x_src_offsets);
 
139
  g_free (y_src_offsets);
 
140
  g_free (src);
 
141
  g_free (dest);
 
142
}
 
143
 
 
144
 
 
145
static void
 
146
get_premultiplied_double_row (PixelRegion *srcPR,
 
147
                              gint         x,
 
148
                              gint         y,
 
149
                              gint         w,
 
150
                              gdouble     *row,
 
151
                              guchar      *tmp_src,
 
152
                              gint         n)
 
153
{
 
154
  const gint bytes = srcPR->bytes;
 
155
  gint       b;
 
156
 
 
157
  pixel_region_get_row (srcPR, x, y, w, tmp_src, n);
 
158
 
 
159
  if (pixel_region_has_alpha (srcPR))
 
160
    {
 
161
      /* premultiply the alpha into the double array */
 
162
      const gint  alpha = bytes - 1;
 
163
      gdouble    *irow  = row;
 
164
 
 
165
      for (x = 0; x < w; x++)
 
166
        {
 
167
          gdouble  mod_alpha = tmp_src[alpha] / 255.0;
 
168
 
 
169
          for (b = 0; b < alpha; b++)
 
170
            irow[b] = mod_alpha * tmp_src[b];
 
171
 
 
172
          irow[b] = tmp_src[alpha];
 
173
          irow += bytes;
 
174
          tmp_src += bytes;
 
175
        }
 
176
    }
 
177
  else /* no alpha */
 
178
    {
 
179
      for (x = 0; x < w * bytes; x++)
 
180
        row[x] = tmp_src[x];
 
181
    }
 
182
 
 
183
  /* set the off edge pixels to their nearest neighbor */
 
184
  for (b = 0; b < 2 * bytes; b++)
 
185
    row[b - 2 * bytes] = row[b % bytes];
 
186
 
 
187
  for (b = 0; b < 2 * bytes; b++)
 
188
    row[b + w * bytes] = row[(w - 1) * bytes + b % bytes];
 
189
}
 
190
 
 
191
 
 
192
static void
 
193
expand_line (gdouble               *dest,
 
194
             const gdouble         *src,
 
195
             gint                   bpp,
 
196
             gint                   old_width,
 
197
             gint                   width,
 
198
             GimpInterpolationType  interp)
 
199
{
 
200
  const gdouble *s;
 
201
  const gdouble  ratio = (gdouble) old_width / (gdouble) width; /* ie reverse scaling_factor */
 
202
  gint           x, b;
 
203
  gint           src_col;
 
204
  gdouble        frac;
 
205
 
 
206
  /* we can overflow src's boundaries, so we expect our caller to have
 
207
  allocated extra space for us to do so safely (see scale_region ()) */
 
208
 
 
209
  switch (interp)
 
210
    {
 
211
      /* -0.5 is because cubic() interpolates a position between 2nd and 3rd data points
 
212
       * we are assigning to 2nd in dest, hence mean shift of +0.5
 
213
       * +1, -1 ensures we dont (int) a negative; first src col only.
 
214
       */
 
215
    case GIMP_INTERPOLATION_CUBIC:
 
216
      for (x = 0; x < width; x++)
 
217
        {
 
218
          gdouble xr = x * ratio - 0.5;
 
219
 
 
220
          if (xr < 0)
 
221
            src_col = (gint) (xr + 1) - 1;
 
222
          else
 
223
            src_col = (gint) xr;
 
224
 
 
225
          frac = xr - src_col;
 
226
          s = &src[src_col * bpp];
 
227
 
 
228
          for (b = 0; b < bpp; b++)
 
229
            dest[b] = cubic_spline_fit (frac, s[b - bpp], s[b], s[b + bpp],
 
230
                                        s[b + bpp * 2]);
 
231
 
 
232
          dest += bpp;
 
233
        }
 
234
 
 
235
      break;
 
236
 
 
237
      /* -0.5 corrects the drift from averaging between adjacent points and assigning to dest[b]
 
238
       * +1, -1 ensures we dont (int) a negative; first src col only.
 
239
       */
 
240
    case GIMP_INTERPOLATION_LINEAR:
 
241
      for (x = 0; x < width; x++)
 
242
        {
 
243
          gdouble xr = (x * ratio + 1 - 0.5) - 1;
 
244
 
 
245
          src_col = (gint) xr;
 
246
          frac = xr - src_col;
 
247
          s = &src[src_col * bpp];
 
248
 
 
249
          for (b = 0; b < bpp; b++)
 
250
            dest[b] = ((s[b + bpp] - s[b]) * frac + s[b]);
 
251
 
 
252
          dest += bpp;
 
253
        }
 
254
      break;
 
255
 
 
256
    case GIMP_INTERPOLATION_NONE:
 
257
    case GIMP_INTERPOLATION_LANCZOS:
 
258
      g_assert_not_reached ();
 
259
      break;
 
260
 
 
261
    default:
 
262
      break;
 
263
    }
 
264
}
 
265
 
 
266
 
 
267
 
 
268
static void
 
269
shrink_line (gdouble               *dest,
 
270
             const gdouble         *src,
 
271
             gint                   bytes,
 
272
             gint                   old_width,
 
273
             gint                   width,
 
274
             GimpInterpolationType  interp)
 
275
{
 
276
  const gdouble *srcp;
 
277
  gdouble       *destp;
 
278
  gdouble        accum[4];
 
279
  gdouble        slice;
 
280
  const gdouble  avg_ratio = (gdouble) width / old_width;
 
281
  const gdouble  inv_width = 1.0 / width;
 
282
  gint           slicepos;      /* slice position relative to width */
 
283
  gint           x;
 
284
  gint           b;
 
285
 
 
286
#if 0
 
287
  g_printerr ("shrink_line bytes=%d old_width=%d width=%d interp=%d "
 
288
              "avg_ratio=%f\n",
 
289
              bytes, old_width, width, interp, avg_ratio);
 
290
#endif
 
291
 
 
292
  g_return_if_fail (bytes <= 4);
 
293
 
 
294
  /* This algorithm calculates the weighted average of pixel data that
 
295
     each output pixel must receive, taking into account that it always
 
296
     scales down, i.e. there's always more than one input pixel per each
 
297
     output pixel.  */
 
298
 
 
299
  srcp = src;
 
300
  destp = dest;
 
301
 
 
302
  slicepos = 0;
 
303
 
 
304
  /* Initialize accum to the first pixel slice.  As there is no partial
 
305
     pixel at start, that value is 0.  The source data is interleaved, so
 
306
     we maintain BYTES accumulators at the same time to deal with that
 
307
     many channels simultaneously.  */
 
308
  for (b = 0; b < bytes; b++)
 
309
    accum[b] = 0.0;
 
310
 
 
311
  for (x = 0; x < width; x++)
 
312
    {
 
313
      /* Accumulate whole pixels.  */
 
314
      do
 
315
        {
 
316
          for (b = 0; b < bytes; b++)
 
317
            accum[b] += *srcp++;
 
318
 
 
319
          slicepos += width;
 
320
        }
 
321
      while (slicepos < old_width);
 
322
      slicepos -= old_width;
 
323
 
 
324
      if (! (slicepos < width))
 
325
        g_warning ("Assertion (slicepos < width) failed. Please report.");
 
326
 
 
327
      if (slicepos == 0)
 
328
        {
 
329
          /* Simplest case: we have reached a whole pixel boundary.  Store
 
330
             the average value per channel and reset the accumulators for
 
331
             the next round.
 
332
 
 
333
             The main reason to treat this case separately is to avoid an
 
334
             access to out-of-bounds memory for the first pixel.  */
 
335
          for (b = 0; b < bytes; b++)
 
336
            {
 
337
              *destp++ = accum[b] * avg_ratio;
 
338
              accum[b] = 0.0;
 
339
            }
 
340
        }
 
341
      else
 
342
        {
 
343
          for (b = 0; b < bytes; b++)
 
344
            {
 
345
              /* We have accumulated a whole pixel per channel where just a
 
346
                 slice of it was needed.  Subtract now the previous pixel's
 
347
                 extra slice.  */
 
348
              slice = srcp[- bytes + b] * slicepos * inv_width;
 
349
              *destp++ = (accum[b] - slice) * avg_ratio;
 
350
 
 
351
              /* That slice is the initial value for the next round.  */
 
352
              accum[b] = slice;
 
353
            }
 
354
        }
 
355
    }
 
356
 
 
357
  /* Sanity check: srcp should point to the next-to-last position, and
 
358
     slicepos should be zero.  */
 
359
  if (! (srcp - src == old_width * bytes && slicepos == 0))
 
360
    g_warning ("Assertion (srcp - src == old_width * bytes && slicepos == 0)"
 
361
               " failed. Please report.");
 
362
}
 
363
 
 
364
static inline void
 
365
rotate_pointers (guchar  **p,
 
366
                 guint32   n)
 
367
{
 
368
  guchar  *tmp = p[0];
 
369
  guint32  i;
 
370
 
 
371
  for (i = 0; i < n-1; i++)
 
372
    p[i] = p[i+1];
 
373
 
 
374
  p[i] = tmp;
 
375
}
 
376
 
 
377
static void
 
378
get_scaled_row (gdouble              **src,
 
379
                gint                   y,
 
380
                gint                   new_width,
 
381
                PixelRegion           *srcPR,
 
382
                gdouble               *row,
 
383
                guchar                *src_tmp,
 
384
                GimpInterpolationType  interpolation_type)
 
385
{
 
386
  /* get the necesary lines from the source image, scale them,
 
387
     and put them into src[] */
 
388
 
 
389
  rotate_pointers ((gpointer) src, 4);
 
390
 
 
391
  if (y < 0)
 
392
    y = 0;
 
393
 
 
394
  if (y < srcPR->h)
 
395
    {
 
396
      get_premultiplied_double_row (srcPR, 0, y, srcPR->w, row, src_tmp, 1);
 
397
 
 
398
      if (new_width > srcPR->w)
 
399
        expand_line (src[3], row, srcPR->bytes,
 
400
                     srcPR->w, new_width, interpolation_type);
 
401
      else if (srcPR->w > new_width)
 
402
        shrink_line (src[3], row, srcPR->bytes,
 
403
                     srcPR->w, new_width, interpolation_type);
 
404
      else /* no scailing needed */
 
405
        memcpy (src[3], row, sizeof (gdouble) * new_width * srcPR->bytes);
 
406
    }
 
407
  else
 
408
    {
 
409
      memcpy (src[3], src[2], sizeof (gdouble) * new_width * srcPR->bytes);
 
410
    }
 
411
}
 
412
 
 
413
void
 
414
scale_region (PixelRegion           *srcPR,
 
415
              PixelRegion           *destPR,
 
416
              GimpInterpolationType  interpolation,
 
417
              GimpProgressFunc       progress_callback,
 
418
              gpointer               progress_data)
 
419
{
 
420
  gdouble *src[4];
 
421
  guchar  *src_tmp;
 
422
  guchar  *dest;
 
423
  gdouble *row, *accum;
 
424
  gint     bytes, b;
 
425
  gint     width, height;
 
426
  gint     orig_width, orig_height;
 
427
  gdouble  y_ratio;
 
428
  gint     i;
 
429
  gint     old_y = -4;
 
430
  gint     new_y;
 
431
  gint     x, y;
 
432
 
 
433
  if (interpolation == GIMP_INTERPOLATION_NONE)
 
434
    {
 
435
      scale_region_no_resample (srcPR, destPR);
 
436
      return;
 
437
    }
 
438
 
 
439
  orig_width = srcPR->w;
 
440
  orig_height = srcPR->h;
 
441
 
 
442
  width = destPR->w;
 
443
  height = destPR->h;
 
444
 
 
445
  if (interpolation == GIMP_INTERPOLATION_LANCZOS && orig_height <= height)
 
446
    {
 
447
      scale_region_lanczos (srcPR, destPR, progress_callback, progress_data);
 
448
      return;
 
449
    }
 
450
 
 
451
#if 0
 
452
  g_printerr ("scale_region: (%d x %d) -> (%d x %d)\n",
 
453
              orig_width, orig_height, width, height);
 
454
#endif
 
455
 
 
456
  /*  find the ratios of old y to new y  */
 
457
  y_ratio = (gdouble) orig_height / (gdouble) height;
 
458
 
 
459
  bytes = destPR->bytes;
 
460
 
 
461
  /*  the data pointers...  */
 
462
  for (i = 0; i < 4; i++)
 
463
    src[i] = g_new (gdouble, width * bytes);
 
464
 
 
465
  dest = g_new (guchar, width * bytes);
 
466
 
 
467
  src_tmp = g_new (guchar, orig_width * bytes);
 
468
 
 
469
  /* offset the row pointer by 2*bytes so the range of the array
 
470
     is [-2*bytes] to [(orig_width + 2)*bytes] */
 
471
  row = g_new (gdouble, (orig_width + 2 * 2) * bytes);
 
472
  row += bytes * 2;
 
473
 
 
474
  accum = g_new (gdouble, width * bytes);
 
475
 
 
476
  /*  Scale the selected region  */
 
477
 
 
478
  for (y = 0; y < height; y++)
 
479
    {
 
480
      if (progress_callback && (y % 64 == 0))
 
481
        progress_callback (0, height, y, progress_data);
 
482
 
 
483
      if (height < orig_height)
 
484
        {
 
485
          const gdouble inv_ratio = 1.0 / y_ratio;
 
486
          gint          max;
 
487
          gdouble       frac;
 
488
 
 
489
          if (y == 0) /* load the first row if this is the first time through */
 
490
            get_scaled_row (&src[0], 0, width, srcPR, row,
 
491
                            src_tmp,
 
492
                            interpolation);
 
493
 
 
494
          new_y = (int) (y * y_ratio);
 
495
          frac = 1.0 - (y * y_ratio - new_y);
 
496
 
 
497
          for (x = 0; x < width * bytes; x++)
 
498
            accum[x] = src[3][x] * frac;
 
499
 
 
500
          max = (int) ((y + 1) * y_ratio) - new_y - 1;
 
501
 
 
502
          get_scaled_row (&src[0], ++new_y, width, srcPR, row,
 
503
                          src_tmp,
 
504
                          interpolation);
 
505
 
 
506
          while (max > 0)
 
507
            {
 
508
              for (x = 0; x < width * bytes; x++)
 
509
                accum[x] += src[3][x];
 
510
 
 
511
              get_scaled_row (&src[0], ++new_y, width, srcPR, row,
 
512
                              src_tmp,
 
513
                              interpolation);
 
514
              max--;
 
515
            }
 
516
 
 
517
          frac = (y + 1) * y_ratio - ((int) ((y + 1) * y_ratio));
 
518
 
 
519
          for (x = 0; x < width * bytes; x++)
 
520
            {
 
521
              accum[x] += frac * src[3][x];
 
522
              accum[x] *= inv_ratio;
 
523
            }
 
524
        }
 
525
      else if (height > orig_height)
 
526
        {
 
527
          new_y = floor (y * y_ratio - 0.5);
 
528
 
 
529
          while (old_y <= new_y)
 
530
            {
 
531
              /* get the necesary lines from the source image, scale them,
 
532
                 and put them into src[] */
 
533
              get_scaled_row (&src[0], old_y + 2, width, srcPR, row,
 
534
                              src_tmp,
 
535
                              interpolation);
 
536
              old_y++;
 
537
            }
 
538
 
 
539
          switch (interpolation)
 
540
            {
 
541
            case GIMP_INTERPOLATION_CUBIC:
 
542
              {
 
543
                gdouble p0, p1, p2, p3;
 
544
                gdouble dy = (y * y_ratio - 0.5) - new_y;
 
545
 
 
546
                p0 = cubic_spline_fit (dy, 1, 0, 0, 0);
 
547
                p1 = cubic_spline_fit (dy, 0, 1, 0, 0);
 
548
                p2 = cubic_spline_fit (dy, 0, 0, 1, 0);
 
549
                p3 = cubic_spline_fit (dy, 0, 0, 0, 1);
 
550
 
 
551
                for (x = 0; x < width * bytes; x++)
 
552
                  accum[x] = (p0 * src[0][x] + p1 * src[1][x] +
 
553
                              p2 * src[2][x] + p3 * src[3][x]);
 
554
              }
 
555
 
 
556
              break;
 
557
 
 
558
            case GIMP_INTERPOLATION_LINEAR:
 
559
              {
 
560
                gdouble idy = (y * y_ratio - 0.5) - new_y;
 
561
                gdouble dy  = 1.0 - idy;
 
562
 
 
563
                for (x = 0; x < width * bytes; x++)
 
564
                  accum[x] = dy * src[1][x] + idy * src[2][x];
 
565
              }
 
566
 
 
567
              break;
 
568
 
 
569
            case GIMP_INTERPOLATION_NONE:
 
570
            case GIMP_INTERPOLATION_LANCZOS:
 
571
              g_assert_not_reached ();
 
572
              break;
 
573
            }
 
574
        }
 
575
      else /* height == orig_height */
 
576
        {
 
577
          get_scaled_row (&src[0], y, width, srcPR, row,
 
578
                          src_tmp,
 
579
                          interpolation);
 
580
          memcpy (accum, src[3], sizeof (gdouble) * width * bytes);
 
581
        }
 
582
 
 
583
      if (pixel_region_has_alpha (srcPR))
 
584
        {
 
585
          /* unmultiply the alpha */
 
586
          gdouble  inv_alpha;
 
587
          gdouble *p     = accum;
 
588
          gint     alpha = bytes - 1;
 
589
          gint     result;
 
590
          guchar  *d     = dest;
 
591
 
 
592
          for (x = 0; x < width; x++)
 
593
            {
 
594
              if (p[alpha] > 0.001)
 
595
                {
 
596
                  inv_alpha = 255.0 / p[alpha];
 
597
 
 
598
                  for (b = 0; b < alpha; b++)
 
599
                    {
 
600
                      result = RINT (inv_alpha * p[b]);
 
601
 
 
602
                      if (result < 0)
 
603
                        d[b] = 0;
 
604
                      else if (result > 255)
 
605
                        d[b] = 255;
 
606
                      else
 
607
                        d[b] = result;
 
608
                    }
 
609
 
 
610
                  result = RINT (p[alpha]);
 
611
 
 
612
                  if (result > 255)
 
613
                    d[alpha] = 255;
 
614
                  else
 
615
                    d[alpha] = result;
 
616
                }
 
617
              else /* alpha <= 0 */
 
618
                {
 
619
                  for (b = 0; b <= alpha; b++)
 
620
                    d[b] = 0;
 
621
                }
 
622
 
 
623
              d += bytes;
 
624
              p += bytes;
 
625
            }
 
626
        }
 
627
      else
 
628
        {
 
629
          gint w = width * bytes;
 
630
 
 
631
          for (x = 0; x < w; x++)
 
632
            {
 
633
              if (accum[x] < 0.0)
 
634
                dest[x] = 0;
 
635
              else if (accum[x] > 255.0)
 
636
                dest[x] = 255;
 
637
              else
 
638
                dest[x] = RINT (accum[x]);
 
639
            }
 
640
        }
 
641
 
 
642
      pixel_region_set_row (destPR, 0, y, width, dest);
 
643
    }
 
644
 
 
645
  /*  free up temporary arrays  */
 
646
  g_free (accum);
 
647
 
 
648
  for (i = 0; i < 4; i++)
 
649
    g_free (src[i]);
 
650
 
 
651
  g_free (src_tmp);
 
652
  g_free (dest);
 
653
 
 
654
  row -= 2 * bytes;
 
655
  g_free (row);
 
656
}
 
657
 
 
658
void
 
659
subsample_region (PixelRegion *srcPR,
 
660
                  PixelRegion *destPR,
 
661
                  gint         subsample)
 
662
{
 
663
  const gint     width       = destPR->w;
 
664
  const gint     height      = destPR->h;
 
665
  const gint     orig_width  = srcPR->w / subsample;
 
666
  const gint     orig_height = srcPR->h / subsample;
 
667
  const gdouble  x_ratio     = (gdouble) orig_width / (gdouble) width;
 
668
  const gdouble  y_ratio     = (gdouble) orig_height / (gdouble) height;
 
669
  const gint     bytes       = destPR->bytes;
 
670
  const gint     destwidth   = destPR->rowstride;
 
671
  guchar        *src,  *s;
 
672
  guchar        *dest, *d;
 
673
  gdouble       *row,  *r;
 
674
  gint           src_row, src_col;
 
675
  gdouble        x_sum, y_sum;
 
676
  gdouble        x_last, y_last;
 
677
  gdouble       *x_frac, y_frac, tot_frac;
 
678
  gint           i, j;
 
679
  gint           b;
 
680
  gint           frac;
 
681
  gint           advance_dest;
 
682
 
 
683
#if 0
 
684
  g_printerr ("subsample_region: (%d x %d) -> (%d x %d)\n",
 
685
              orig_width, orig_height, width, height);
 
686
#endif
 
687
 
 
688
  /*  the data pointers...  */
 
689
  src  = g_new (guchar, orig_width * bytes);
 
690
  dest = destPR->data;
 
691
 
 
692
  /*  allocate an array to help with the calculations  */
 
693
  row    = g_new (gdouble, width * bytes);
 
694
  x_frac = g_new (gdouble, width + orig_width);
 
695
 
 
696
  /*  initialize the pre-calculated pixel fraction array  */
 
697
  src_col = 0;
 
698
  x_sum = (gdouble) src_col;
 
699
  x_last = x_sum;
 
700
 
 
701
  for (i = 0; i < width + orig_width; i++)
 
702
    {
 
703
      if (x_sum + x_ratio <= (src_col + 1 + EPSILON))
 
704
        {
 
705
          x_sum += x_ratio;
 
706
          x_frac[i] = x_sum - x_last;
 
707
        }
 
708
      else
 
709
        {
 
710
          src_col ++;
 
711
          x_frac[i] = src_col - x_last;
 
712
        }
 
713
 
 
714
      x_last += x_frac[i];
 
715
    }
 
716
 
 
717
  /*  clear the "row" array  */
 
718
  memset (row, 0, sizeof (gdouble) * width * bytes);
 
719
 
 
720
  /*  counters...  */
 
721
  src_row = 0;
 
722
  y_sum = (gdouble) src_row;
 
723
  y_last = y_sum;
 
724
 
 
725
  pixel_region_get_row (srcPR,
 
726
                        srcPR->x, srcPR->y + src_row * subsample,
 
727
                        orig_width * subsample,
 
728
                        src, subsample);
 
729
 
 
730
  /*  Scale the selected region  */
 
731
  for (i = 0; i < height; )
 
732
    {
 
733
      src_col = 0;
 
734
      x_sum = (gdouble) src_col;
 
735
 
 
736
      /* determine the fraction of the src pixel we are using for y */
 
737
      if (y_sum + y_ratio <= (src_row + 1 + EPSILON))
 
738
        {
 
739
          y_sum += y_ratio;
 
740
          y_frac = y_sum - y_last;
 
741
          advance_dest = TRUE;
 
742
        }
 
743
      else
 
744
        {
 
745
          src_row ++;
 
746
          y_frac = src_row - y_last;
 
747
          advance_dest = FALSE;
 
748
        }
 
749
 
 
750
      y_last += y_frac;
 
751
 
 
752
      s = src;
 
753
      r = row;
 
754
 
 
755
      frac = 0;
 
756
      j = width;
 
757
 
 
758
      while (j)
 
759
        {
 
760
          tot_frac = x_frac[frac++] * y_frac;
 
761
 
 
762
          for (b = 0; b < bytes; b++)
 
763
            r[b] += s[b] * tot_frac;
 
764
 
 
765
          /*  increment the destination  */
 
766
          if (x_sum + x_ratio <= (src_col + 1 + EPSILON))
 
767
            {
 
768
              r += bytes;
 
769
              x_sum += x_ratio;
 
770
              j--;
 
771
            }
 
772
 
 
773
          /* increment the source */
 
774
          else
 
775
            {
 
776
              s += bytes;
 
777
              src_col++;
 
778
            }
 
779
        }
 
780
 
 
781
      if (advance_dest)
 
782
        {
 
783
          tot_frac = 1.0 / (x_ratio * y_ratio);
 
784
 
 
785
          /*  copy "row" to "dest"  */
 
786
          d = dest;
 
787
          r = row;
 
788
 
 
789
          j = width;
 
790
          while (j--)
 
791
            {
 
792
              b = bytes;
 
793
 
 
794
              while (b--)
 
795
                *d++ = (guchar) (*r++ * tot_frac + 0.5);
 
796
            }
 
797
 
 
798
          dest += destwidth;
 
799
 
 
800
          /*  clear the "row" array  */
 
801
          memset (row, 0, sizeof (gdouble) * destwidth);
 
802
 
 
803
          i++;
 
804
        }
 
805
      else
 
806
        {
 
807
          pixel_region_get_row (srcPR,
 
808
                                srcPR->x,
 
809
                                srcPR->y + src_row * subsample,
 
810
                                orig_width * subsample,
 
811
                                src, subsample);
 
812
        }
 
813
    }
 
814
 
 
815
  /*  free up temporary arrays  */
 
816
  g_free (row);
 
817
  g_free (x_frac);
 
818
  g_free (src);
 
819
}
 
820
 
 
821
 
 
822
/* Lanczos */
 
823
static inline gdouble
 
824
sinc (gdouble x)
 
825
{
 
826
  gdouble y = x * G_PI;
 
827
 
 
828
  if (ABS (x) < LANCZOS_MIN)
 
829
    return 1.0;
 
830
 
 
831
  return sin (y) / y;
 
832
}
 
833
 
 
834
static inline gdouble
 
835
lanczos_sum (guchar         *ptr,
 
836
             const gdouble  *kernel,  /* 1-D kernel of transform coeffs */
 
837
             gint            u,
 
838
             gint            bytes,
 
839
             gint            byte)
 
840
{
 
841
  gdouble sum = 0;
 
842
  gint    i;
 
843
 
 
844
  for (i = 0; i < LANCZOS_WIDTH2 ; i++)
 
845
    sum += kernel[i] * ptr[ (u + i - LANCZOS_WIDTH) * bytes + byte ];
 
846
 
 
847
  return sum;
 
848
}
 
849
 
 
850
static inline gdouble
 
851
lanczos_sum_mul (guchar         *ptr,
 
852
                 const gdouble  *kernel,    /* 1-D kernel of transform coeffs */
 
853
                 gint            u,
 
854
                 gint            bytes,
 
855
                 gint            byte,
 
856
                 gint            alpha)
 
857
{
 
858
  gdouble sum = 0;
 
859
  gint    i;
 
860
 
 
861
  for (i = 0; i < LANCZOS_WIDTH2 ; i++ )
 
862
    sum += kernel[i] * ptr[ (u + i - LANCZOS_WIDTH) * bytes + byte ]
 
863
                 * ptr[ (u + i - LANCZOS_WIDTH) * bytes + alpha];
 
864
 
 
865
  return sum;
 
866
}
 
867
 
 
868
static gboolean
 
869
inv_lin_trans (const gdouble *t,
 
870
               gdouble       *it)
 
871
{
 
872
  gdouble d = (t[0] * t[4]) - (t[1] * t[3]);  /* determinant */
 
873
 
 
874
  if (fabs(d) < EPSILON )
 
875
    return FALSE;
 
876
 
 
877
  it[0] =    t[4] / d;
 
878
  it[1] =   -t[1] / d;
 
879
  it[2] = (( t[1] * t[5]) - (t[2] * t[4])) / d;
 
880
  it[3] =   -t[3] / d;
 
881
  it[4] =    t[0] / d;
 
882
  it[5] = (( t[2] * t[3]) - (t[0] * t[5])) / d;
 
883
 
 
884
  return TRUE;
 
885
}
 
886
 
 
887
 
 
888
/*
 
889
 * allocate and fill lookup table of Lanczos windowed sinc function
 
890
 * use gfloat since errors due to granularity of array far exceed data precision
 
891
 */
 
892
gfloat *
 
893
create_lanczos_lookup (void)
 
894
{
 
895
  const gdouble dx = LANCZOS_WIDTH / (gdouble) (LANCZOS_SAMPLES - 1);
 
896
 
 
897
  gfloat  *lookup = g_new (gfloat, LANCZOS_SAMPLES);
 
898
  gdouble  x      = 0.0;
 
899
  gint     i;
 
900
 
 
901
  for (i = 0; i < LANCZOS_SAMPLES; i++)
 
902
    {
 
903
      lookup[i] = ((ABS (x) < LANCZOS_WIDTH) ?
 
904
                   (sinc (x) * sinc (x / LANCZOS_WIDTH)) : 0.0);
 
905
      x += dx;
 
906
    }
 
907
 
 
908
  return lookup;
 
909
}
 
910
 
 
911
static void
 
912
scale_region_lanczos (PixelRegion      *srcPR,
 
913
                      PixelRegion      *dstPR,
 
914
                      GimpProgressFunc  progress_callback,
 
915
                      gpointer          progress_data)
 
916
 
 
917
{
 
918
  gfloat        *kernel_lookup  = NULL;             /* Lanczos lookup table                */
 
919
  gdouble        x_kernel[LANCZOS_WIDTH2],    /* 1-D kernels of Lanczos window coeffs */
 
920
                 y_kernel[LANCZOS_WIDTH2];
 
921
 
 
922
  gdouble        newval;                      /* new interpolated RGB value          */
 
923
 
 
924
  guchar        *win_buf = NULL;              /* Sliding window buffer               */
 
925
  guchar        *win_ptr[LANCZOS_WIDTH2];     /* Ponters to sliding window rows      */
 
926
 
 
927
  guchar        *dst_buf = NULL;              /* Pointer to destination image data   */
 
928
 
 
929
  gint           x, y;                        /* Position in destination image       */
 
930
  gint           i, byte;                     /* loop vars  */
 
931
  gint           row;
 
932
 
 
933
  gdouble        trans[6], itrans[6];         /* Scale transformations               */
 
934
 
 
935
  const gint     dst_width     = dstPR->w;
 
936
  const gint     dst_height    = dstPR->h;
 
937
  const gint     bytes         = dstPR->bytes;
 
938
  const gint     src_width     = srcPR->w;
 
939
  const gint     src_height    = srcPR->h;
 
940
 
 
941
  const gint     src_row_span  = src_width * bytes;
 
942
  const gint     dst_row_span  = dst_width * bytes;
 
943
  const gint     win_row_span  = (src_width + LANCZOS_WIDTH2) * bytes;
 
944
 
 
945
  const gdouble  scale_x       = dst_width / (gdouble) src_width;
 
946
  const gdouble  scale_y       = dst_height / (gdouble) src_height;
 
947
 
 
948
  for (i = 0; i < 6; i++)
 
949
    trans[i] = 0.0;
 
950
 
 
951
  trans[0] = scale_x;
 
952
  trans[4] = scale_y;
 
953
 
 
954
  if (! inv_lin_trans (trans, itrans))
 
955
    {
 
956
      g_warning ("transformation matrix is not invertible");
 
957
      return;
 
958
    }
 
959
 
 
960
  /* allocate buffer for destination row */
 
961
  dst_buf = g_new0 (guchar, dst_row_span);
 
962
 
 
963
  /* if no scaling needed copy data */
 
964
  if (dst_width == src_width && dst_height == src_height)
 
965
    {
 
966
       for (i = 0 ; i < src_height ; i++)
 
967
         {
 
968
           pixel_region_get_row (srcPR, 0, i, src_width, dst_buf, 1);
 
969
           pixel_region_set_row (dstPR, 0, i, dst_width, dst_buf);
 
970
         }
 
971
       g_free(dst_buf);
 
972
       return;
 
973
    }
 
974
 
 
975
  /* allocate and fill kernel_lookup lookup table */
 
976
  kernel_lookup = create_lanczos_lookup ();
 
977
 
 
978
  /* allocate buffer for source rows */
 
979
  win_buf = g_new0 (guchar, win_row_span * LANCZOS_WIDTH2);
 
980
 
 
981
  /* Set the window pointers */
 
982
  for ( i = 0 ; i < LANCZOS_WIDTH2 ; i++ )
 
983
    win_ptr[i] = win_buf + (win_row_span * i) + LANCZOS_WIDTH * bytes;
 
984
 
 
985
  /* fill the data for the first loop */
 
986
  for ( i = 0 ; i <= LANCZOS_WIDTH && i < src_height ; i++)
 
987
    pixel_region_get_row (srcPR, 0, i, src_width, win_ptr[i + LANCZOS_WIDTH], 1);
 
988
 
 
989
  for (row = y = 0; y < dst_height; y++)
 
990
    {
 
991
      if (progress_callback && (y % 64 == 0))
 
992
        progress_callback (0, dst_height, y, progress_data);
 
993
 
 
994
      pixel_region_get_row (dstPR, 0, y, dst_width, dst_buf, 1);
 
995
      for  (x = 0; x < dst_width; x++)
 
996
        {
 
997
          gdouble  dsrc_x ,dsrc_y;        /* corresponding scaled position in source image */
 
998
          gint     int_src_x, int_src_y;  /* integer part of coordinates in source image   */
 
999
          gint     x_shift, y_shift;      /* index into Lanczos lookup                     */
 
1000
          gdouble  kx_sum, ky_sum;        /* sums of Lanczos kernel coeffs                 */
 
1001
 
 
1002
          /* -0.5 corrects image drift.due to average offset used in lookup */
 
1003
          dsrc_x = x / scale_x - 0.5;
 
1004
          dsrc_y = y / scale_y - 0.5;
 
1005
 
 
1006
          /* avoid (int) on negative*/
 
1007
          if (dsrc_x > 0)
 
1008
            int_src_x = (gint) (dsrc_x);
 
1009
          else
 
1010
            int_src_x = (gint) (dsrc_x + 1) - 1;
 
1011
 
 
1012
          if (dsrc_y > 0)
 
1013
            int_src_y = (gint) (dsrc_y);
 
1014
          else
 
1015
            int_src_y = (gint) (dsrc_y + 1) - 1;
 
1016
 
 
1017
          /* calc lookup offsets for non-interger remainders */
 
1018
          x_shift = (gint) ((dsrc_x - int_src_x) * LANCZOS_SPP + 0.5);
 
1019
          y_shift = (gint) ((dsrc_y - int_src_y) * LANCZOS_SPP + 0.5);
 
1020
 
 
1021
          /*  Fill x_kernel[] and y_kernel[] with lanczos coeffs
 
1022
           *
 
1023
           *  kernel_lookup = Is a lookup table that contains half of the symetrical windowed-sinc func.
 
1024
           *
 
1025
           *  x_shift, y_shift = shift from kernel center due to fractional part
 
1026
           *           of interpollation
 
1027
           *
 
1028
           *  The for-loop creates two 1-D kernels for convolution.
 
1029
           *    - If the center position +/- LANCZOS_WIDTH is out of
 
1030
           *      the source image coordinates set the value to 0.0
 
1031
           *      FIXME => partial kernel. Define a more rigourous border mode.
 
1032
           *    - If the kernel index is out of range set value to 0.0
 
1033
           *      ( caused by offset coeff. obselete??)
 
1034
           */
 
1035
          kx_sum = ky_sum = 0.0;
 
1036
 
 
1037
          for (i = LANCZOS_WIDTH; i >= -LANCZOS_WIDTH; i--)
 
1038
            {
 
1039
              gint pos = i * LANCZOS_SPP;
 
1040
 
 
1041
              if ( int_src_x + i >= 0 && int_src_x + i < src_width)
 
1042
                kx_sum += x_kernel[LANCZOS_WIDTH + i] = kernel_lookup[ABS (x_shift - pos)];
 
1043
              else
 
1044
                x_kernel[LANCZOS_WIDTH + i] = 0.0;
 
1045
 
 
1046
              if ( int_src_y + i >= 0 && int_src_y + i < src_height)
 
1047
                ky_sum += y_kernel[LANCZOS_WIDTH + i] = kernel_lookup[ABS (y_shift - pos)];
 
1048
              else
 
1049
                y_kernel[LANCZOS_WIDTH + i] = 0.0;
 
1050
            }
 
1051
 
 
1052
          /* normalise the kernel arrays */
 
1053
          for (i = -LANCZOS_WIDTH; i <= LANCZOS_WIDTH; i++)
 
1054
          {
 
1055
             x_kernel[LANCZOS_WIDTH + i] /= kx_sum;
 
1056
             y_kernel[LANCZOS_WIDTH + i] /= ky_sum;
 
1057
          }
 
1058
 
 
1059
          /*
 
1060
            Scaling up
 
1061
            New determined source row is > than last read row
 
1062
            rotate the pointers and get next source row from region.
 
1063
            If no more source rows are available fill buffer with 0
 
1064
            ( Probably not necessary because multipliers should be 0).
 
1065
          */
 
1066
          for ( ; row < int_src_y ; )
 
1067
            {
 
1068
              row++;
 
1069
              rotate_pointers (win_ptr, LANCZOS_WIDTH2);
 
1070
              if ( row + LANCZOS_WIDTH < src_height)
 
1071
                pixel_region_get_row (srcPR, 0,
 
1072
                                      row + LANCZOS_WIDTH, src_width,
 
1073
                                      win_ptr[LANCZOS_WIDTH2 - 1], 1);
 
1074
              else
 
1075
                memset (win_ptr[LANCZOS_WIDTH2 - 1], 0,
 
1076
                        sizeof (guchar) * src_row_span);
 
1077
            }
 
1078
           /*
 
1079
              Scaling down
 
1080
            */
 
1081
            for ( ; row > int_src_y ; )
 
1082
            {
 
1083
              row--;
 
1084
              for ( i = 0 ; i < LANCZOS_WIDTH2 - 1 ; i++ )
 
1085
                 rotate_pointers (win_ptr, LANCZOS_WIDTH2);
 
1086
              if ( row >=  0)
 
1087
                pixel_region_get_row (srcPR, 0,
 
1088
                                      row, src_width,
 
1089
                                      win_ptr[0], 1);
 
1090
              else
 
1091
                memset (win_ptr[0], 0,
 
1092
                        sizeof (guchar) * src_row_span);
 
1093
 
 
1094
            }
 
1095
 
 
1096
 
 
1097
           if (pixel_region_has_alpha (srcPR))
 
1098
             {
 
1099
               const gint alpha = bytes - 1;
 
1100
               gint       byte;
 
1101
               gdouble    arecip;
 
1102
               gdouble    aval;
 
1103
 
 
1104
               aval = 0.0;
 
1105
               for (i = 0; i < LANCZOS_WIDTH2 ; i++ )
 
1106
                 aval += y_kernel[i] * lanczos_sum (win_ptr[i], x_kernel,
 
1107
                         int_src_x, bytes, alpha);
 
1108
 
 
1109
               if (aval <= 0.0)
 
1110
                 {
 
1111
                   arecip = 0.0;
 
1112
                   dst_buf[x * bytes + alpha] = 0;
 
1113
                 }
 
1114
               else if (aval > 255.0)
 
1115
                 {
 
1116
                   arecip = 1.0 / aval;
 
1117
                   dst_buf[x * bytes + alpha] = 255;
 
1118
                 }
 
1119
               else
 
1120
                 {
 
1121
                   arecip = 1.0 / aval;
 
1122
                   dst_buf[x * bytes + alpha] = RINT (aval);
 
1123
                 }
 
1124
 
 
1125
               for (byte = 0; byte < alpha; byte++)
 
1126
                 {
 
1127
                   newval = 0.0;
 
1128
                   for (i = 0; i < LANCZOS_WIDTH2; i++ )
 
1129
                     newval += y_kernel[i] * lanczos_sum_mul (win_ptr[i], x_kernel,
 
1130
                               int_src_x, bytes, byte, alpha);
 
1131
                   newval *= arecip;
 
1132
                   dst_buf[x * bytes + byte] = CLAMP (newval, 0, 255);
 
1133
                 }
 
1134
             }
 
1135
           else
 
1136
             {
 
1137
               for (byte = 0; byte < bytes; byte++)
 
1138
                 {
 
1139
                   /* Calculate new value */
 
1140
                   newval = 0.0;
 
1141
                   for (i = 0; i < LANCZOS_WIDTH2; i++ )
 
1142
                     newval += y_kernel[i] * lanczos_sum (win_ptr[i], x_kernel,
 
1143
                               int_src_x, bytes, byte);
 
1144
                   dst_buf[x * bytes + byte] = CLAMP ((gint) newval, 0, 255);
 
1145
                 }
 
1146
             }
 
1147
        }
 
1148
 
 
1149
      pixel_region_set_row (dstPR, 0, y , dst_width, dst_buf);
 
1150
    }
 
1151
 
 
1152
  g_free (dst_buf);
 
1153
  g_free (win_buf);
 
1154
  g_free (kernel_lookup);
 
1155
}