Diff of /resample.c [000000] .. [2a97a8] Maximize Restore

  Switch to unified view

a b/resample.c
1
/* Panorama_Tools -   Generate, Edit and Convert Panoramic Images
2
   Copyright (C) 1998,1999 - Helmut Dersch  der@fh-furtwangen.de
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, or (at your option)
7
   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., 675 Mass Ave, Cambridge, MA 02139, USA.  */
17
18
/*------------------------------------------------------------*/
19
20
21
22
23
24
// Program specific includes
25
26
#include "filter.h"           
27
28
29
// Standard C includes
30
31
#include <stdio.h>
32
#include <stdarg.h>
33
#include <stdlib.h>
34
#include <string.h>
35
#include <math.h>
36
37
38
//            This file uses functions of type
39
//    resample( unsigned char *dst,   unsigned char **rgb,
40
//                            register double Dx, 
41
//                            register double Dy,
42
//                            int color, int SamplesPerPixel);
43
//
44
// dst - output pixel
45
// rgb - input pixels, may be Lab as well.
46
// Dx  - offset of output pixel position in x-direction
47
// Dy  - offset of output pixel position in y-direction
48
// color = 0: all rgb colors; color = 1,2,3: one of r,g,b
49
// BytesPerPixel = 3,4. Using color != 0, any value should (?) work.
50
51
52
53
// Arrays used for Gamma correction
54
55
56
PTGamma glu; // Lookup table
57
58
59
60
// Some locally needed math functions
61
62
static double     sinc        ( double x );
63
static double     cubic01     ( double x );
64
static double     cubic12     ( double x ); 
65
66
67
68
// Interpolators
69
70
static void nn( unsigned char *dst, unsigned char **rgb, 
71
      register double Dx, register double Dy, int color, int SamplesPerPixel);
72
      
73
static void bil( unsigned char *dst, unsigned char **rgb,  
74
      register double Dx, register double Dy, int color, int SamplesPerPixel);
75
      
76
static void poly3( unsigned char *dst, unsigned char **rgb,  
77
      register double Dx, register double Dy, int color, int SamplesPerPixel);
78
      
79
static void spline16( unsigned char *dst, unsigned char **rgb,  
80
      register double Dx, register double Dy, int color, int SamplesPerPixel);
81
      
82
static void spline36( unsigned char *dst, unsigned char **rgb,  
83
      register double Dx, register double Dy, int color, int SamplesPerPixel);
84
      
85
static void spline64( unsigned char *dst, unsigned char **rgb,  
86
      register double Dx, register double Dy, int color, int SamplesPerPixel);
87
88
static void sinc256( unsigned char *dst, unsigned char **rgb,  
89
      register double Dx, register double Dy, int color, int SamplesPerPixel);
90
91
static void sinc1024( unsigned char *dst, unsigned char **rgb,  
92
      register double Dx, register double Dy, int color, int SamplesPerPixel);
93
94
static void nn_16( unsigned char *dst, unsigned char **rgb, 
95
      register double Dx, register double Dy, int color, int SamplesPerPixel);
96
      
97
static void bil_16( unsigned char *dst, unsigned char **rgb,  
98
      register double Dx, register double Dy, int color, int SamplesPerPixel);
99
      
100
static void poly3_16( unsigned char *dst, unsigned char **rgb,  
101
      register double Dx, register double Dy, int color, int SamplesPerPixel);
102
      
103
static void spline16_16( unsigned char *dst, unsigned char **rgb,  
104
      register double Dx, register double Dy, int color, int SamplesPerPixel);
105
      
106
static void spline36_16( unsigned char *dst, unsigned char **rgb,  
107
      register double Dx, register double Dy, int color, int SamplesPerPixel);
108
      
109
static void spline64_16( unsigned char *dst, unsigned char **rgb,  
110
      register double Dx, register double Dy, int color, int SamplesPerPixel);
111
112
static void sinc256_16( unsigned char *dst, unsigned char **rgb,  
113
      register double Dx, register double Dy, int color, int SamplesPerPixel);
114
115
static void sinc1024_16( unsigned char *dst, unsigned char **rgb,  
116
      register double Dx, register double Dy, int color, int SamplesPerPixel);
117
118
119
120
121
122
123
// Various interpolators; a[] is array of coeeficients; 0 <= x < 1
124
125
126
#define   NNEIGHBOR(x, a , NDIM )                                         \
127
          a[0] = 1.0; 
128
129
130
#define   BILINEAR(x, a, NDIM )                                           \
131
          a[1] = x;                                                       \
132
          a[0] = 1.0 - x; 
133
          
134
135
// Unused; has been replaced by 'CUBIC'.
136
137
#define   POLY3( x, a , NDIM )                                            \
138
          a[3] = (  x * x - 1.0) * x / 6.0;                               \
139
          a[2] = ( (1.0 - x) * x / 2.0 + 1.0) * x;                        \
140
          a[1] = ( ( 1.0/2.0  * x - 1.0 ) * x - 1.0/2.0 ) * x + 1.0;      \
141
          a[0] = ( ( -1.0/6.0 * x + 1.0/2.0 ) * x - 1.0/3.0 ) * x ;
142
143
144
145
#define       SPLINE16( x, a, NDIM )                                          \
146
          a[3] = ( ( 1.0/3.0  * x - 1.0/5.0 ) * x -   2.0/15.0 ) * x;     \
147
          a[2] = ( ( 6.0/5.0 - x     ) * x +   4.0/5.0 ) * x;             \
148
          a[1] = ( ( x - 9.0/5.0 ) * x -   1.0/5.0     ) * x + 1.0;       \
149
          a[0] = ( ( -1.0/3.0 * x + 4.0/5.0     ) * x -   7.0/15.0 ) * x ;
150
151
152
#define       CUBIC( x, a, NDIM )                                         \
153
          a[3] = cubic12( 2.0 - x);                                   \
154
          a[2] = cubic01( 1.0 - x);                                   \
155
          a[1] = cubic01( x );                                        \
156
          a[0] = cubic12( x + 1.0);                                   \
157
158
159
160
161
#define       SPLINE36( x, a , NDIM )                                                     \
162
  a[5] = ( ( -  1.0/11.0  * x +  12.0/ 209.0 ) * x +   7.0/ 209.0  ) * x;             \
163
  a[4] = ( (    6.0/11.0  * x -  72.0/ 209.0 ) * x -  42.0/ 209.0  ) * x;             \
164
  a[3] = ( ( - 13.0/11.0  * x + 288.0/ 209.0 ) * x + 168.0/ 209.0  ) * x;             \
165
  a[2] = ( (   13.0/11.0  * x - 453.0/ 209.0 ) * x -   3.0/ 209.0  ) * x + 1.0;       \
166
  a[1] = ( ( -  6.0/11.0  * x + 270.0/ 209.0 ) * x - 156.0/ 209.0  ) * x;             \
167
  a[0] = ( (    1.0/11.0  * x -  45.0/ 209.0 ) * x +  26.0/ 209.0  ) * x;
168
169
170
171
#define       SPLINE64( x, a , NDIM )                                                     \
172
  a[7] = ((  1.0/41.0 * x -   45.0/2911.0) * x -   26.0/2911.0) * x;                  \
173
  a[6] = ((- 6.0/41.0 * x +  270.0/2911.0) * x +  156.0/2911.0) * x;                  \
174
  a[5] = (( 24.0/41.0 * x - 1080.0/2911.0) * x -  624.0/2911.0) * x;                  \
175
  a[4] = ((-49.0/41.0 * x + 4050.0/2911.0) * x + 2340.0/2911.0) * x;                  \
176
  a[3] = (( 49.0/41.0 * x - 6387.0/2911.0) * x -    3.0/2911.0) * x + 1.0;            \
177
  a[2] = ((-24.0/41.0 * x + 4032.0/2911.0) * x - 2328.0/2911.0) * x;                  \
178
  a[1] = ((  6.0/41.0 * x - 1008.0/2911.0) * x +  582.0/2911.0) * x;                  \
179
  a[0] = ((- 1.0/41.0 * x +  168.0/2911.0) * x -   97.0/2911.0) * x;                  
180
181
182
#define       SINC( x, a, NDIM )                                      \
183
  {                                                               \
184
      register int idx;                                           \
185
      register double xadd;                                       \
186
      for( idx = 0, xadd = NDIM / 2 - 1.0 + x;                    \
187
           idx < NDIM / 2;                                        \
188
           xadd-=1.0)                                             \
189
      {                                                           \
190
          a[idx++] = sinc( xadd ) * sinc( xadd / ( NDIM / 2 ));   \
191
      }                                                           \
192
      for( xadd = 1.0 - x;                                        \
193
           idx < NDIM;                                            \
194
           xadd+=1.0)                                             \
195
      {                                                           \
196
          a[idx++] = sinc( xadd ) * sinc( xadd / ( NDIM / 2 ));   \
197
      }                                                           \
198
  }                                                               \
199
      
200
201
202
203
204
205
206
207
208
// Set up the arrays for gamma correction
209
210
int SetUpGamma( double pgamma, int psize)
211
{
212
  int i;
213
  double gnorm, xg, rgamma = 1.0/pgamma;
214
215
  if( psize == 1 )
216
  {
217
      glu.ChannelSize     =   256;
218
      glu.ChannelStretch  =    16;
219
  }
220
  else if( psize == 2 )
221
  {
222
      glu.ChannelSize     = 65536;
223
      glu.ChannelStretch  =     4;
224
  }
225
  else
226
      return -1;
227
228
  glu.GammaSize = glu.ChannelSize * glu.ChannelStretch;
229
  
230
  glu.DeGamma     = NULL;
231
  glu.Gamma       = NULL;
232
  glu.DeGamma     = (double*)         malloc( glu.ChannelSize * sizeof( double ) );
233
  glu.Gamma       = (unsigned short*) malloc( glu.GammaSize * sizeof( unsigned short) );
234
  
235
  if( glu.DeGamma == NULL || glu.Gamma == NULL )
236
  {
237
      PrintError("Not enough memory");
238
      return -1;
239
  }
240
241
  glu.DeGamma[0] = 0.0;
242
  gnorm = (glu.ChannelSize-1) / pow( glu.ChannelSize-1 , pgamma ) ; 
243
  for(i=1; i<glu.ChannelSize; i++)
244
  {
245
      glu.DeGamma[i] = pow( (double)i , pgamma ) * gnorm;
246
  }
247
248
  glu.Gamma[0] = 0;
249
  gnorm = (glu.ChannelSize-1) /  pow( glu.ChannelSize-1 , rgamma ) ; 
250
  if( psize == 1 )
251
  {
252
      for(i=1; i<glu.GammaSize; i++)
253
      {
254
          xg   = pow(  ((double)i) / glu.ChannelStretch , rgamma ) * gnorm;
255
          DBL_TO_UC( glu.Gamma[i], xg );
256
      }
257
  }
258
  else
259
  {
260
      for(i=1; i<glu.GammaSize; i++)
261
      {
262
          xg   = pow(  ((double)i) / glu.ChannelStretch , rgamma ) * gnorm;
263
          DBL_TO_US( glu.Gamma[i], xg );
264
      }
265
  }
266
  return 0;
267
}
268
269
unsigned short gamma_correct( double pix )
270
{
271
  int k = glu.ChannelStretch * pix;
272
  if( k < 0 )
273
      return 0;
274
  if( k > glu.GammaSize - 1 )
275
      return glu.ChannelSize - 1;
276
  return (glu.Gamma)[ k ] ;
277
}
278
279
280
281
282
/////////// N x N Sampler /////////////////////////////////////////////
283
284
#define RESAMPLE_N( intpol, ndim, psize )                             \
285
  double yr[ndim], yg[ndim], yb[ndim], w[ndim];                       \
286
  register double rd, gd, bd, weight ;                                \
287
  register int k,i;                                                   \
288
  register unsigned psize *r, *ri;                                    \
289
                                                                      \
290
  intpol( Dx, w, ndim )                                               \
291
  if( color == 0 )                                                    \
292
  {                                                                   \
293
      for(k=0; k<ndim; k++)                                           \
294
      {                                                               \
295
          r = ((unsigned psize**)rgb)[k]  + SamplesPerPixel - 3;      \
296
          rd = gd = bd = 0.0;                                         \
297
                                                                      \
298
          for(i=0; i<ndim; i++)                                       \
299
          {                                                           \
300
              weight = w[ i ];                                        \
301
              ri     = r + i * SamplesPerPixel;                       \
302
              rd += glu.DeGamma[(int)*ri++] * weight;                 \
303
              gd += glu.DeGamma[(int)*ri++] * weight;                 \
304
              bd += glu.DeGamma[(int)*ri]   * weight;                 \
305
          }                                                           \
306
          yr[k] = rd; yg[k] = gd; yb[k] = bd;                         \
307
      }                                                               \
308
                                                                      \
309
      intpol( Dy, w, ndim )                                           \
310
      rd = gd = bd = 0.0;                                             \
311
                                                                      \
312
      for(i=0; i<ndim; i++)                                           \
313
      {                                                               \
314
          weight = w[ i ];                                            \
315
          rd += yr[i] * weight;                                       \
316
          gd += yg[i] * weight;                                       \
317
          bd += yb[i] * weight;                                       \
318
      }                                                               \
319
      *((unsigned psize*)dst)++   =   gamma_correct( rd );            \
320
      *((unsigned psize*)dst)++   =   gamma_correct( gd );            \
321
      *((unsigned psize*)dst)     =   gamma_correct( bd );            \
322
  }                                                                   \
323
  else                                                                \
324
  {                                                                   \
325
      color-=1;                                                       \
326
      for(k=0; k<ndim; k++)                                           \
327
      {                                                               \
328
          r = ((unsigned psize**)rgb)[k] + SamplesPerPixel - 3 + color;\
329
          yr[k] =  0.0;                                               \
330
                                                                      \
331
          for(i=0; i<ndim; i++)                                       \
332
          {                                                           \
333
              yr[k] += glu.DeGamma[(int)r[i*SamplesPerPixel]] * w[i];\
334
          }                                                           \
335
      }                                                               \
336
                                                                      \
337
      intpol( Dy, w, ndim )                                           \
338
      rd = 0.0;                                                       \
339
                                                                      \
340
      for(i=0; i<ndim; i++)                                           \
341
      {                                                               \
342
          rd += yr[i] * w[ i ];                                       \
343
      }                                                               \
344
       *((unsigned psize*)dst+color)  =   gamma_correct( rd );        \
345
  }                                                                   \
346
347
348
349
350
static double sinc( double x )
351
{
352
  x *= PI;
353
  if(x != 0.0) 
354
      return(sin(x) / x);
355
  return(1.0);
356
}
357
358
359
// Cubic polynomial with parameter A
360
// A = -1: sharpen; A = - 0.5 homogeneous
361
// make sure x >= 0
362
#define   A   (-0.75)
363
364
// 0 <= x < 1
365
static double cubic01( double x )
366
{
367
  return  (( A + 2.0 )*x - ( A + 3.0 ))*x*x +1.0;
368
}
369
// 1 <= x < 2
370
371
static double cubic12( double x )
372
{
373
  return  (( A * x - 5.0 * A ) * x + 8.0 * A ) * x - 4.0 * A;
374
375
}
376
377
#undef A
378
379
380
381
382
383
// ---------- Sampling functions ----------------------------------
384
385
386
// Nearest neighbor sampling, nowhere used (yet)
387
388
static void nn( unsigned char *dst, unsigned char **rgb, 
389
      register double Dx, register double Dy, int color, int SamplesPerPixel)
390
      {
391
#pragma unused(Dx)
392
#pragma unused(Dy)
393
          RESAMPLE_N( NNEIGHBOR, 1, char) }
394
395
// Bilinear sampling, nowhere used (yet).
396
397
static void bil( unsigned char *dst, unsigned char **rgb,  
398
      register double Dx, register double Dy, int color, int SamplesPerPixel)
399
      {   RESAMPLE_N( BILINEAR, 2, char)  }
400
401
402
// Lowest quality sampler in distribution; since version 1.8b1 changed to closely
403
// resemble Photoshop's bicubic interpolation
404
405
static void poly3( unsigned char *dst, unsigned char **rgb,  
406
      register double Dx, register double Dy, int color, int SamplesPerPixel)
407
      {   RESAMPLE_N( CUBIC, 4, char)     }
408
409
410
// Spline using 16 pixels; smoother and less artefacts than poly3, softer; same speed
411
412
static void spline16( unsigned char *dst, unsigned char **rgb,  
413
      register double Dx, register double Dy, int color, int SamplesPerPixel)
414
      {   RESAMPLE_N( SPLINE16, 4, char)  }
415
416
// Spline using 36 pixels; significantly sharper than both poly3 and spline16,
417
// almost no artefacts
418
419
static void spline36( unsigned char *dst, unsigned char **rgb,  
420
      register double Dx, register double Dy, int color, int SamplesPerPixel)
421
      {   RESAMPLE_N( SPLINE36, 6, char)  }
422
423
// Not used anymore
424
425
static void spline64( unsigned char *dst, unsigned char **rgb,  
426
      register double Dx, register double Dy, int color, int SamplesPerPixel)
427
      {   RESAMPLE_N( SPLINE64, 8, char)  }
428
429
430
// Highest quality sampler since version 1.8b1
431
// Extremely slow, but defintely worth every second.
432
433
static void sinc256( unsigned char *dst, unsigned char **rgb,  
434
      register double Dx, register double Dy, int color, int SamplesPerPixel)
435
      {   RESAMPLE_N( SINC, 16, char)     }
436
      
437
438
// Highest quality sampler since version 1.8b1
439
// Extremely slow, but defintely worth every second.
440
441
static void sinc1024( unsigned char *dst, unsigned char **rgb,  
442
      register double Dx, register double Dy, int color, int SamplesPerPixel)
443
      {   RESAMPLE_N( SINC, 32, char)     }
444
445
446
//--------------- Same as above, for shorts (16 bit channel size-------------------
447
448
// Nearest neighbor sampling, nowhere used (yet)
449
450
static void nn_16( unsigned char *dst, unsigned char **rgb, 
451
      register double Dx, register double Dy, int color, int SamplesPerPixel)
452
      {
453
#pragma unused(Dx)
454
#pragma unused(Dy)
455
          RESAMPLE_N( NNEIGHBOR, 1, short)    }
456
457
// Bilinear sampling, nowhere used (yet).
458
459
static void bil_16( unsigned char *dst, unsigned char **rgb,  
460
      register double Dx, register double Dy, int color, int SamplesPerPixel)
461
      {   RESAMPLE_N( BILINEAR, 2, short)     }
462
463
464
// Lowest quality sampler in distribution; since version 1.8b1 changed to closely
465
// resemble Photoshop's bicubic interpolation
466
467
static void poly3_16( unsigned char *dst, unsigned char **rgb,  
468
      register double Dx, register double Dy, int color, int SamplesPerPixel)
469
      {   RESAMPLE_N( CUBIC, 4, short)    }
470
471
472
// Spline using 16 pixels; smoother and less artefacts than poly3, softer; same speed
473
474
static void spline16_16( unsigned char *dst, unsigned char **rgb,  
475
      register double Dx, register double Dy, int color, int SamplesPerPixel)
476
      {   RESAMPLE_N( SPLINE16, 4, short)     }
477
478
// Spline using 36 pixels; significantly sharper than both poly3 and spline16,
479
// almost no artefacts
480
481
static void spline36_16( unsigned char *dst, unsigned char **rgb,  
482
      register double Dx, register double Dy, int color, int SamplesPerPixel)
483
      {   RESAMPLE_N( SPLINE36, 6, short)     }
484
485
// Not used anymore
486
487
static void spline64_16( unsigned char *dst, unsigned char **rgb,  
488
      register double Dx, register double Dy, int color, int SamplesPerPixel)
489
      {   RESAMPLE_N( SPLINE64, 8, short)     }
490
491
492
// Highest quality sampler since version 1.8b1
493
// Extremely slow, but defintely worth every second.
494
495
static void sinc256_16( unsigned char *dst, unsigned char **rgb,  
496
      register double Dx, register double Dy, int color, int SamplesPerPixel)
497
      {   RESAMPLE_N( SINC, 16, short)    }
498
      
499
500
// Highest quality sampler since version 1.8b1
501
// Extremely slow, but defintely worth every second.
502
503
static void sinc1024_16( unsigned char *dst, unsigned char **rgb,  
504
      register double Dx, register double Dy, int color, int SamplesPerPixel)
505
      {   RESAMPLE_N( SINC, 32, short)    }
506
      
507
508
509
510
//    Main transformation function. Destination image is calculated using transformation
511
//    Function "func". Either all colors (color = 0) or one of rgb (color =1,2,3) are
512
//    determined. If successful, TrPtr->success = 1. Memory for destination image
513
//    must have been allocated and locked!
514
515
void transForm( TrformStr *TrPtr, fDesc *fD, int color){
516
  register int        x, y;       // Loop through destination image
517
  register int        i, k;       // Auxilliary loop variables
518
  int             skip = 0;   // Update progress counter
519
  unsigned char       *dest,*src,*sry;// Source and destination image data
520
  register unsigned char      *sr;    // Source  image data
521
  char*           progressMessage;// Message to be displayed by progress reporter
522
  char                    percent[8]; // Number displayed by Progress reporter
523
  int         valid;      // Is this pixel valid? (i.e. inside source image)
524
  long            coeff;      // pixel coefficient in destination image
525
  long            cy;     // rownum in destimage
526
  int         xc,yc;
527
528
  double          x_d, y_d;   // Cartesian Coordinates of point ("target") in Destination image
529
  double          Dx, Dy;     // Coordinates of target in Source image
530
  int             xs, ys; 
531
532
  unsigned char       **rgb  = NULL, 
533
              *cdata = NULL;  // Image data handed to sampler
534
535
  double          max_x = (double) TrPtr->src->width; // Maximum x values in source image
536
  double          max_y = (double) TrPtr->src->height; // Maximum y values in source image
537
  double          min_x =  -1.0;//0.0; // Minimum x values in source image
538
  double          min_y =  -1.0;//0.0; // Minimum y values in source image
539
540
  int         mix   = TrPtr->src->width - 1; // maximum x-index src
541
  int         mix2;
542
  int         miy   = TrPtr->src->height - 1;// maximum y-index src
543
  int         miy2;
544
545
  // Variables used to convert screen coordinates to cartesian coordinates
546
547
      
548
  double          w2  = (double) TrPtr->dest->width  / 2.0 - 0.5;  // Steve's L
549
  double          h2  = (double) TrPtr->dest->height / 2.0 - 0.5;
550
  double          sw2 = (double) TrPtr->src->width   / 2.0 - 0.5;
551
  double          sh2 = (double) TrPtr->src->height  / 2.0 - 0.5;
552
  
553
  int         BytesPerLine    = TrPtr->src->bytesPerLine;
554
  int         BytesPerPixel, FirstColorByte, SamplesPerPixel, BytesPerSample;
555
556
  int         n, n2;      // How many pixels should be used for interpolation 
557
  intFunc         intp;       // Function used to interpolate
558
  int             lu = 0;     // Use lookup table?
559
  int         wrap_x = FALSE;
560
  double          theGamma;   // gamma handed to SetUpGamma()
561
562
  // Selection rectangle
563
  PTRect          destRect;
564
  if( TrPtr->dest->selection.bottom == 0 && TrPtr->dest->selection.right == 0 ){
565
      destRect.left   = 0;
566
      destRect.right  = TrPtr->dest->width;
567
      destRect.top    = 0;
568
      destRect.bottom = TrPtr->dest->height;
569
  }else{
570
      memcpy( &destRect, &TrPtr->dest->selection, sizeof(PTRect) );
571
  }
572
573
574
  switch( TrPtr->src->bitsPerPixel ){
575
      case 64:    FirstColorByte = 2; BytesPerPixel = 8; SamplesPerPixel = 4; BytesPerSample = 2; break;
576
      case 48:    FirstColorByte = 0; BytesPerPixel = 6; SamplesPerPixel = 3; BytesPerSample = 2; break;
577
      case 32:    FirstColorByte = 1; BytesPerPixel = 4; SamplesPerPixel = 4; BytesPerSample = 1; break;
578
      case 24:    FirstColorByte = 0; BytesPerPixel = 3; SamplesPerPixel = 3; BytesPerSample = 1; break;
579
      case  8:    FirstColorByte = 0; BytesPerPixel = 1; SamplesPerPixel = 1; BytesPerSample = 1; break;
580
      default:    PrintError("Unsupported Pixel Size: %d", TrPtr->src->bitsPerPixel);
581
                  TrPtr->success = 0;
582
                  return;
583
  }
584
  
585
  // Set interpolator etc:
586
  switch( TrPtr->interpolator ){
587
      case _poly3:// Third order polynomial fitting 16 nearest pixels
588
          if( BytesPerSample == 1 ) intp = poly3; else intp = poly3_16;       
589
          n = 4;
590
          break;
591
      case _spline16:// Cubic Spline fitting 16 nearest pixels
592
          if( BytesPerSample == 1 ) intp = spline16; else intp = spline16_16;     
593
          n = 4;
594
          break;
595
      case _spline36: // Cubic Spline fitting 36 nearest pixels
596
          if( BytesPerSample == 1 ) intp = spline36; else intp = spline36_16;     
597
          n = 6;
598
          break;
599
      case _spline64: // Cubic Spline fitting 64 nearest pixels
600
          if( BytesPerSample == 1 ) intp = spline64; else intp = spline64_16; 
601
          n = 8;
602
          break;
603
      case _sinc256:  // sinc windowed to 256 (2*8)^2 pixels
604
          if( BytesPerSample == 1 ) intp = sinc256; else intp = sinc256_16;   
605
          n = 16;
606
          break;
607
      case _sinc1024: // sinc windowed to 1024 (2*16)^2 pixels
608
          if( BytesPerSample == 1 ) intp = sinc1024; else intp = sinc1024_16; 
609
          n = 32;
610
          break;
611
      case _bilinear: // Bilinear fit using 4 nearest points
612
          if( BytesPerSample == 1 ) intp = bil; else intp = bil_16;   
613
          n = 2;
614
          break;
615
      case _nn:// nearest neighbor fit using 4 nearest points
616
          if( BytesPerSample == 1 ) intp = nn; else intp = nn_16; 
617
          n = 1;
618
          break;
619
      default: 
620
          PrintError( "Invalid Interpolator selected" );
621
          TrPtr->success = 0;
622
          return;
623
  }
624
625
  // Set up arrays that hold color data for interpolators
626
627
  rgb     = (unsigned char**) malloc( n * sizeof(unsigned char*) );
628
  cdata   = (unsigned char*)  malloc( n * n * BytesPerPixel * sizeof( unsigned char ) );
629
  
630
  
631
  if( rgb == NULL || cdata == NULL ){
632
      PrintError( "Not enough Memory" );
633
      TrPtr->success = 0;
634
      goto Trform_exit;
635
  }
636
      
637
  n2 = n/2 ;
638
  mix2 = mix +1 - n;
639
  miy2 = miy +1 - n;
640
641
  dest = *TrPtr->dest->data;
642
  src  = *TrPtr->src->data; // is locked
643
644
  if(TrPtr->mode & _show_progress){
645
      switch(color){
646
          case 0: progressMessage = "Image Conversion";   break;
647
          case 1: switch( TrPtr->src->dataformat){
648
                      case _RGB:  progressMessage = "Red Channel"  ; break;
649
                      case _Lab:  progressMessage = "Lightness"    ; break;
650
                  } break;
651
          case 2: switch( TrPtr->src->dataformat){
652
                      case _RGB:  progressMessage = "Green Channel"; break;
653
                      case _Lab:  progressMessage = "Color A"      ; break;
654
                  } break; 
655
          case 3: switch( TrPtr->src->dataformat){
656
                      case _RGB:  progressMessage = "Blue Channel"; break;
657
                      case _Lab:  progressMessage = "Color B"     ; break;
658
                  } break; 
659
          default: progressMessage = "Something is wrong here";
660
      }
661
      Progress( _initProgress, progressMessage );
662
  }
663
664
  if(TrPtr->mode & _wrapX)
665
      wrap_x = TRUE;
666
667
  if( TrPtr->src->dataformat == _RGB )    // Gamma correct only RGB-images
668
      theGamma = TrPtr->gamma;
669
  else
670
      theGamma = 1.0;
671
  
672
  if( SetUpGamma( theGamma, BytesPerSample) != 0 ){
673
      PrintError( "Could not set up lookup table for Gamma Correction" );
674
      TrPtr->success = 0;
675
      goto Trform_exit;
676
  }
677
678
  for(y=destRect.top; y<destRect.bottom; y++){
679
      // Update Progress report and check for cancel every 5 lines.
680
      skip++;
681
      if( skip == 5 ){
682
          if(TrPtr->mode & _show_progress){   
683
              sprintf( percent, "%d", (int) (y * 100)/ TrPtr->dest->height);
684
              if( ! Progress( _setProgress, percent ) ){
685
                  TrPtr->success = 0;
686
                  goto Trform_exit;
687
              }
688
          }else{
689
              if( ! Progress( _idleProgress, 0) ){
690
                  TrPtr->success = 0;
691
                  goto Trform_exit;
692
              }
693
          }
694
          skip = 0;
695
      }
696
      
697
      // y-coordinate in dest image relative to center        
698
      y_d = (double) y - h2 ;
699
      cy  = (y-destRect.top) * TrPtr->dest->bytesPerLine; 
700
      
701
      for(x=destRect.left; x<destRect.right; x++){
702
          // Calculate pixel coefficient in dest image just once
703
704
          coeff = cy  + BytesPerPixel * (x-destRect.left);        
705
706
          // Convert destination screen coordinates to cartesian coordinates.         
707
          x_d = (double) x - w2 ;
708
          
709
          // Get source cartesian coordinates 
710
          fD->func( x_d, y_d , &Dx, &Dy, fD->param);
711
712
          // Convert source cartesian coordinates to screen coordinates 
713
          Dx += sw2;
714
          Dy =  sh2 + Dy ;
715
          
716
717
          // Is the pixel valid, i.e. from within source image?
718
          if( (Dx >= max_x)   || (Dy >= max_y) || (Dx < min_x) || (Dy < min_y)  )
719
              valid = FALSE;
720
          else
721
              valid = TRUE;
722
723
          // Convert only valid pixels
724
          if( valid ){
725
              // Extract integer and fractions of source screen coordinates
726
              xc    =  (int)floor( Dx ) ; Dx -= (double)xc;
727
              yc    =  (int)floor( Dy ) ; Dy -= (double)yc;
728
              
729
              // if alpha channel marks valid portions, set valid 
730
              if(TrPtr->mode & _honor_valid)
731
              switch( FirstColorByte ){
732
                  case 1:{
733
                      int xn = xc, yn = yc;
734
                      if( xn < 0 ) xn = 0; //  -1  crashes Windows
735
                      if( yn < 0 ) yn = 0; //  -1  crashes Windows
736
                      if( src[ yn * BytesPerLine + BytesPerPixel * xn] == 0 )
737
                          valid = FALSE;
738
                      }
739
                      break;
740
                  case 2:{
741
                      int xn = xc, yn = yc;
742
                      if( xn < 0 ) xn = 0; //  -1  crashes Windows
743
                      if( yn < 0 ) yn = 0; //  -1  crashes Windows
744
                      if( *((USHORT*)(src + yn * BytesPerLine + BytesPerPixel * xn)) == 0 )
745
                          valid = FALSE;
746
                      }
747
                      break;
748
                  default: break;
749
              }
750
          }
751
          
752
          if( valid ){    
753
              ys = yc +1 - n2 ; // smallest y-index used for interpolation
754
              xs = xc +1 - n2 ; // smallest x-index used for interpolation
755
                  
756
              // y indices used: yc-(n2-1)....yc+n2
757
              // x indices used: xc-(n2-1)....xc+n2
758
                  
759
              if( ys >= 0 && ys <= miy2 && xs >= 0 && xs <= mix2 ){  // all interpolation pixels inside image
760
                  sry = src + ys * BytesPerLine + xs * BytesPerPixel;
761
                  for(i = 0;  i < n;  i++, sry += BytesPerLine){
762
                      rgb[i] = sry;
763
                  }
764
              }else{ // edge pixels
765
                  if( ys < 0 )
766
                      sry = src;
767
                  else if( ys > miy )
768
                      sry = src + miy * BytesPerLine;
769
                  else
770
                      sry = src + ys  * BytesPerLine;
771
                  
772
                  for(i = 0; i < n; i++){ 
773
                      xs = xc +1 - n2 ; // smallest x-index used for interpolation
774
                      if( wrap_x ){
775
                          while( xs < 0 )  xs += (mix+1);
776
                          while( xs > mix) xs -= (mix+1);
777
                      }
778
                      if( xs < 0 )
779
                           sr = sry;
780
                      else if( xs > mix )
781
                          sr = sry + mix *BytesPerPixel;
782
                      else
783
                          sr = sry + BytesPerPixel * xs;
784
                  
785
                      rgb[i] = cdata + i * n * BytesPerPixel;
786
                      for(k = 0; k < n; k++ ){
787
                          memcpy( &(rgb[i][k * BytesPerPixel]), sr, BytesPerPixel);
788
                          xs++;
789
                          if( wrap_x ){
790
                              while( xs < 0 )  xs += (mix+1);
791
                              while( xs > mix) xs -= (mix+1);
792
                          }
793
                          if( xs < 0 )
794
                              sr = sry;
795
                          else if( xs > mix )
796
                              sr = sry + mix *BytesPerPixel;
797
                          else
798
                              sr = sry + BytesPerPixel * xs;
799
                      }
800
                       ys++;
801
                       if( ys > 0 && ys <= miy )
802
                          sry +=  BytesPerLine; 
803
                  }
804
805
              }
806
              
807
              switch( FirstColorByte ){
808
                  case 1: dest[ coeff++ ] = UCHAR_MAX;         // Set alpha channel
809
                          break;
810
                  case 2: *((USHORT*)(dest + coeff)) = USHRT_MAX; coeff+=2;
811
                          break;
812
                  default: break;
813
              }
814
                      
815
              intp( &(dest[ coeff ]), rgb, Dx, Dy, color, SamplesPerPixel ); 
816
817
          }else{  // not valid
818
              memset( &(dest[ coeff ]), 0 ,BytesPerPixel ); 
819
          }
820
821
      }
822
  }
823
824
//    if(TrPtr->mode & _show_progress){
825
//        Progress( _disposeProgress, percent );
826
//    }   
827
  TrPtr->success = 1;
828
829
830
Trform_exit:
831
  if( rgb )       free( rgb );
832
  if( cdata )         free( cdata );
833
  if( glu.DeGamma )   free( glu.DeGamma );    glu.DeGamma     = NULL;
834
  if( glu.Gamma )     free( glu.Gamma );  glu.Gamma   = NULL;
835
  return;
836
}
837
  
838
839
840
841
#if 0
842
843
// An unused lookup version of sinc256
844
// Somewhat  faster on non-floating point machines like Intel
845
846
static double*   SetUpWeights(  );
847
848
849
// Weigths for sinc function
850
851
#define   NUM_WEIGHTS 256
852
853
static double *wt = NULL;
854
855
856
857
858
#define       SINC256( x, a )                     \
859
  if( wt == NULL ) wt = SetUpWeights( );      \
860
  if( wt != NULL )                            \
861
  {                                           \
862
      int xn = x * NUM_WEIGHTS;               \
863
      a[15]   = wt[ 8*NUM_WEIGHTS - xn ];     \
864
      a[14]   = wt[ 7*NUM_WEIGHTS - xn ];     \
865
      a[13]   = wt[ 6*NUM_WEIGHTS - xn ];     \
866
      a[12]   = wt[ 5*NUM_WEIGHTS - xn ];     \
867
      a[11]   = wt[ 4*NUM_WEIGHTS - xn ];     \
868
      a[10]   = wt[ 3*NUM_WEIGHTS - xn ];     \
869
      a[ 9]   = wt[ 2*NUM_WEIGHTS - xn ];     \
870
      a[ 8]   = wt[ 1*NUM_WEIGHTS - xn ];     \
871
      a[ 7]   = wt[ 0*NUM_WEIGHTS + xn ];     \
872
      a[ 6]   = wt[ 1*NUM_WEIGHTS + xn ];     \
873
      a[ 5]   = wt[ 2*NUM_WEIGHTS + xn ];     \
874
      a[ 4]   = wt[ 3*NUM_WEIGHTS + xn ];     \
875
      a[ 3]   = wt[ 4*NUM_WEIGHTS + xn ];     \
876
      a[ 2]   = wt[ 5*NUM_WEIGHTS + xn ];     \
877
      a[ 1]   = wt[ 6*NUM_WEIGHTS + xn ];     \
878
      a[ 0]   = wt[ 7*NUM_WEIGHTS + xn ];     \
879
  }                                           \
880
881
// Create Weights for A * NUM_WEIGHTS positions
882
883
static double* SetUpWeights(  )
884
{
885
#define A     8
886
  int i,k,id;
887
  double dx = 1.0 / (double)NUM_WEIGHTS;
888
  double *w;
889
  
890
  w = (double*)malloc( A * NUM_WEIGHTS * sizeof(double) );
891
  if( w )
892
  {
893
      for( k=0; k < A ; k++ )
894
      {
895
          id = k * NUM_WEIGHTS;
896
          for( i=0; i<NUM_WEIGHTS; i++ )
897
          {
898
              w[id + i] = sinc8( (double)k + i*dx );
899
          }
900
      }
901
  }
902
  return w;
903
}
904
      
905
#undef A
906
907
908
#endif        
909
910
911
912
913
#if 0
914
/////////////// Results of calc for poly3  ////////////////////////////////////////////////////////////   
915
916
Equations:
917
918
1:  y0 = -a3 + a2 - a1 + a0 
919
2:  y1 = a0
920
3:  y2 = a3 + a2 + a1 + a0
921
4:  y3 = 8 a3 + 4 a2 + 2 a1 + a0
922
923
924
--- Emacs Calculator Mode ---
925
4:  a1 = y2 - y3 / 6 - y0 / 3 - y1 / 2
926
3:  6 a3 = y3 - y0 + 3 y1 - 3 y2
927
2:  2 a2 = y2 + y0 - 2 y1
928
1:  a0 = y1
929
930
931
/////////////// Results of Calc for Spline 16 ////////////////////////////////////////////////////////////    
932
933
1.  y0 = -c3 + c2 - c1 + c0
934
2.  y1 = c0
935
3.  y1 = a0
936
4.  y2 = a3 + a2 + a1 + a0
937
5.  y2 = b3 + b2 + b1 + b0
938
6.  y3 = 8 b3 + 4 b2 + 2 b1 + b0
939
7.  c1 = a1
940
8.  3 a3 + 2 a2 + a1 = 3 b3 + 2 b2 + b1
941
9.  2 c2 = 2 a2
942
10. 6 a3 + 2 a2 = 6 b3 + 2 b2
943
11. -6 c3 + 2 c2 = 0
944
12. 12 b3 + 2 b2 = 0
945
946
/////////////// Results of Calc for Spline 36 ////////////////////////////////////////////////////////////    
947
948
Equations:
949
950
--- Emacs Calculator Mode ---
951
20: 18 c3 + 2 c2 = 0
952
19: 2 e2 - 12 e3 = 0
953
18: 12 b3 + 2 b2 = 12 c3 + 2 c2
954
17: 6 a3 + 2 a2 = 6 b3 + 2 b2
955
16: 2 d2 = 2 a2
956
15: 2 e2 - 6 e3 = 2 d2 - 6 d3
957
14: 12 b3 + 4 b2 + b1 = 12 c3 + 4 c2 + c1
958
13: 3 a3 + 2 a2 + a1 = 3 b3 + 2 b2 + b1
959
12: d1 = a1
960
11: 3 e3 - 2 e2 + e1 = 3 d3 - 2 d2 + d1
961
10: 27 c3 + 9 c2 + 3 c1 + c0 = y5
962
9:  8 c3 + 4 c2 + 2 c1 + c0 = y4
963
8:  y4 = 8 b3 + 4 b2 + 2 b1 + b0
964
7:  y3 = b3 + b2 + b1 + b0
965
6:  y3 = a3 + a2 + a1 + a0
966
5:  y0 = 4 e2 - 8 e3 - 2 e1 + e0
967
4:  y1 = e2 - e3 - e1 + e0
968
3:  y1 = d2 - d3 - d1 + d0
969
2:  y2 = d0
970
1:  y2 = a0
971
972
--- Emacs Calculator Mode ---
973
4:  11 a3 = 6 y4 - y5 - 13 y3 - 6 y1 + y0 + 13 y2
974
3:  209 a1 
975
      = 7 y5 + 168 y3 - 42 y4 - 3 y2 + 26 y0 
976
          - 156 y1
977
2:  a2 = 12:209 y5 + 288:209 y3 - 72:209 y4 
978
           - 45:209 y0 - 453:209 y2 + 270:209 y1
979
1:  a0 = y2
980
981
982
/////////////// Results of Calc for Spline 64 ////////////////////////////////////////////////////////////    
983
984
Equations:
985
986
1:   y0 = -27 g3 + 9 g2 - 3 g1 + g0
987
2:   y1 = -8 g3 + 4 g2 - 2 g1 + g0
988
3:   y1 = -8 f3 + 4 f2 - 2 f1 + f0
989
4:   y2 = -f3 + f2 - f1 + f0
990
5:   y2 = -e3 + e2 - e1 + e0
991
6:   y3 = e0
992
7:   y3 = a0
993
8:   y4 = a3 + a2 + a1 + a0
994
9:   y4 = b3 + b2 + b1 + b0
995
10:  y5 = 8 b3 + 4 b2 + 2 b1 + b0
996
11:  y5 = 8 c3 + 4 c2 + 2 c1 + c0
997
12:  y6 = 27 c3 + 9 c2 + 3 c1 + c0
998
13:  y6 = 27 d3 + 9 d2 + 3 d1 + d0
999
14:  y7 = 64 d3 + 16 d2 + 4 d1 + d0
1000
15:  12 g3 - 4 g2 + g1 = 12 f3 - 4 f2 + f1
1001
16:  3 f3 - 2 f2 + f1 = 3 e3 - 2 e2 + e1
1002
17:  e1 = a1
1003
18   3 a3 + 2 a2 + a1 = 3 b3 + 2 b2 + b1
1004
19   12 b3 + 4 b2 + b1 = 12 c3 + 4 c2 + c1
1005
20   27 c3 + 6 c2 + c1 = 27 d3 + 6 d2 + d1
1006
21   -12 g3 + 2 g2 = -12 f3 + 2 f2
1007
22   -6 f3 + 2 f2 = -6 e3 + 2 e2 
1008
23   2 e2 = 2 a2
1009
24   6 a3 + 2 a2 = 6 b3 + 2 b2
1010
25   12 b3 + 2 b2 = 12 c3 + 2 c2
1011
26   18 c3 + 2 c2 = 18 d3 + 2 d2
1012
27   -18 g3 + 2 g2 = 0
1013
28   24 d3 + 2 d2 = 0
1014
1015
--- Emacs Calculator Mode ---
1016
4:  41 a3 
1017
      = y7 + 24 y5 - 6 y6 - 49 y4 - 24 y2 - y0 
1018
          + 49 y3 + 6 y1
1019
3:  a2 = 270:2911 y6 + 4050:2911 y4 - 45:2911 y7 
1020
           - 1080:2911 y5 + 168:2911 y0 
1021
           + 4032:2911 y2 - 1008:2911 y1 
1022
           - 6387:2911 y3
1023
2:  2911 a1 
1024
      = 156 y6 + 2340 y4 - 26 y7 - 624 y5 + 582 y1 
1025
          - 3 y3 - 2328 y2 - 97 y0
1026
1:  a0 = y3
1027
1028
1029
1030
#endif
1031
1032