~ubuntu-branches/ubuntu/gutsy/amsn/gutsy

« back to all changes in this revision

Viewing changes to utils/TkCximage/src/CxImage/ximaint.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Theodore Karkoulis
  • Date: 2006-01-04 15:26:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20060104152602-ipe1yg00rl3nlklv
Tags: 0.95-1
New Upstream Release (closes: #345052, #278575).

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
// xImaInt.cpp : interpolation functions
 
2
/* 02/2004 - Branko Brevensek 
 
3
 * CxImage version 5.99c 17/Oct/2004 - Davide Pizzolato - www.xdp.it
 
4
 */
 
5
 
 
6
#include "ximage.h"
 
7
#include "ximath.h"
 
8
 
 
9
#if CXIMAGE_SUPPORT_INTERPOLATION
 
10
 
 
11
////////////////////////////////////////////////////////////////////////////////
 
12
/**
 
13
 * Recalculates coordinates according to specified overflow method.
 
14
 * If pixel (x,y) lies within image, nothing changes.
 
15
 *
 
16
 *  \param x, y - coordinates of pixel
 
17
 *  \param ofMethod - overflow method
 
18
 * 
 
19
 *  \return x, y - new coordinates (pixel (x,y) now lies inside image)
 
20
 *
 
21
 *  \author ***bd*** 2.2004
 
22
 */
 
23
void CxImage::OverflowCoordinates(long &x, long &y, OverflowMethod const ofMethod)
 
24
{
 
25
  if (IsInside(x,y)) return;  //if pixel is within bounds, no change
 
26
  switch (ofMethod) {
 
27
    case OM_REPEAT:
 
28
      //clip coordinates
 
29
      x=max(x,0); x=min(x, head.biWidth-1);
 
30
      y=max(y,0); y=min(y, head.biHeight-1);
 
31
      break;
 
32
    case OM_WRAP:
 
33
      //wrap coordinates
 
34
      x = x % head.biWidth;
 
35
      y = y % head.biHeight;
 
36
      if (x<0) x = head.biWidth + x;
 
37
      if (y<0) y = head.biHeight + y;
 
38
      break;
 
39
    case OM_MIRROR:
 
40
      //mirror pixels near border
 
41
      if (x<0) x=((-x) % head.biWidth);
 
42
      else if (x>=head.biWidth) x=head.biWidth-(x % head.biWidth + 1);
 
43
      if (y<0) y=((-y) % head.biHeight);
 
44
      else if (y>=head.biHeight) y=head.biHeight-(y % head.biHeight + 1);
 
45
      break;
 
46
    default:
 
47
      return;
 
48
  }//switch
 
49
}
 
50
 
 
51
////////////////////////////////////////////////////////////////////////////////
 
52
/**
 
53
 * See OverflowCoordinates for integer version 
 
54
 * \author ***bd*** 2.2004
 
55
 */
 
56
void CxImage::OverflowCoordinates(float &x, float &y, OverflowMethod const ofMethod)
 
57
{
 
58
  if (x>=0 && x<head.biWidth && y>=0 && y<head.biHeight) return;  //if pixel is within bounds, no change
 
59
  switch (ofMethod) {
 
60
    case OM_REPEAT:
 
61
      //clip coordinates
 
62
      x=max(x,0); x=min(x, head.biWidth-1);
 
63
      y=max(y,0); y=min(y, head.biHeight-1);
 
64
      break;
 
65
    case OM_WRAP:
 
66
      //wrap coordinates
 
67
      x = (float)fmod(x, (float) head.biWidth);
 
68
      y = (float)fmod(y, (float) head.biHeight);
 
69
      if (x<0) x = head.biWidth + x;
 
70
      if (y<0) y = head.biHeight + y;
 
71
      break;
 
72
    case OM_MIRROR:
 
73
      //mirror pixels near border
 
74
      if (x<0) x=(float)fmod(-x, (float) head.biWidth);
 
75
      else if (x>=head.biWidth) x=head.biWidth-((float)fmod(x, (float) head.biWidth) + 1);
 
76
      if (y<0) y=(float)fmod(-y, (float) head.biHeight);
 
77
      else if (y>=head.biHeight) y=head.biHeight-((float)fmod(y, (float) head.biHeight) + 1);
 
78
      break;
 
79
    default:
 
80
      return;
 
81
  }//switch
 
82
}
 
83
 
 
84
////////////////////////////////////////////////////////////////////////////////
 
85
/**
 
86
 * Method return pixel color. Different methods are implemented for out of bounds pixels.
 
87
 * If an image has alpha channel, alpha value is returned in .RGBReserved.
 
88
 *
 
89
 *  \param x,y : pixel coordinates
 
90
 *  \param ofMethod : out-of-bounds method:
 
91
 *    - OF_WRAP - wrap over to pixels on other side of the image
 
92
 *    - OF_REPEAT - repeat last pixel on the edge
 
93
 *    - OF_COLOR - return input value of color
 
94
 *    - OF_BACKGROUND - return background color (if not set, return input color)
 
95
 *    - OF_TRANSPARENT - return transparent pixel
 
96
 *
 
97
 *  \param rplColor : input color (returned for out-of-bound coordinates in OF_COLOR mode and if other mode is not applicable)
 
98
 *
 
99
 * \return color : color of pixel
 
100
 * \author ***bd*** 2.2004
 
101
 */
 
102
RGBQUAD CxImage::GetPixelColorWithOverflow(long x, long y, OverflowMethod const ofMethod, RGBQUAD* const rplColor)
 
103
{
 
104
  RGBQUAD color;          //color to return
 
105
  if ((!IsInside(x,y)) || pDib==NULL) {     //is pixel within bouns?:
 
106
    //pixel is out of bounds or no DIB
 
107
    if (rplColor!=NULL)
 
108
      color=*rplColor;
 
109
    else {
 
110
      color.rgbRed=color.rgbGreen=color.rgbBlue=255; color.rgbReserved=0; //default replacement colour: white transparent
 
111
    }//if
 
112
    if (pDib==NULL) return color;
 
113
    //pixel is out of bounds:
 
114
    switch (ofMethod) {
 
115
      case OM_TRANSPARENT:
 
116
#if CXIMAGE_SUPPORT_ALPHA
 
117
        if (AlphaIsValid()) {
 
118
          //alpha transparency is supported and image has alpha layer
 
119
          color.rgbReserved=0;
 
120
        } else {
 
121
#endif //CXIMAGE_SUPPORT_ALPHA
 
122
          //no alpha transparency
 
123
          if (GetTransIndex()>=0) {
 
124
            color=GetTransColor();    //single color transparency enabled (return transparent color)
 
125
          }//if
 
126
#if CXIMAGE_SUPPORT_ALPHA
 
127
        }//if
 
128
#endif //CXIMAGE_SUPPORT_ALPHA
 
129
        return color;
 
130
      case OM_BACKGROUND:
 
131
                  //return background color (if it exists, otherwise input value)
 
132
                  if (info.nBkgndIndex != -1) {
 
133
                          if (head.biBitCount<24) color = GetPaletteColor((BYTE)info.nBkgndIndex);
 
134
                          else color = info.nBkgndColor;
 
135
                  }//if
 
136
                  return color;
 
137
      case OM_REPEAT:
 
138
      case OM_WRAP:
 
139
      case OM_MIRROR:
 
140
        OverflowCoordinates(x,y,ofMethod);
 
141
        break;
 
142
      default:
 
143
        //simply return replacement color (OM_COLOR and others)
 
144
        return color;
 
145
    }//switch
 
146
  }//if
 
147
  //just return specified pixel (it's within bounds)
 
148
  return BlindGetPixelColor(x,y);
 
149
}
 
150
 
 
151
////////////////////////////////////////////////////////////////////////////////
 
152
/**
 
153
 * This method reconstructs image according to chosen interpolation method and then returns pixel (x,y).
 
154
 * (x,y) can lie between actual image pixels. If (x,y) lies outside of image, method returns value
 
155
 * according to overflow method.
 
156
 * This method is very useful for geometrical image transformations, where destination pixel
 
157
 * can often assume color value lying between source pixels.
 
158
 *
 
159
 *  \param (x,y) - coordinates of pixel to return
 
160
 *           GPCI method recreates "analogue" image back from digital data, so x and y
 
161
 *           are float values and color value of point (1.1,1) will generally not be same
 
162
 *           as (1,1). Center of first pixel is at (0,0) and center of pixel right to it is (1,0).
 
163
 *           (0.5,0) is half way between these two pixels.
 
164
 *  \param inMethod - interpolation (reconstruction) method (kernel) to use:
 
165
 *    - IM_NEAREST_NEIGHBOUR - returns colour of nearest lying pixel (causes stairy look of 
 
166
 *                            processed images)
 
167
 *    - IM_BILINEAR - interpolates colour from four neighbouring pixels (softens image a bit)
 
168
 *    - IM_BICUBIC - interpolates from 16 neighbouring pixels (can produce "halo" artifacts)
 
169
 *    - IM_BICUBIC2 - interpolates from 16 neighbouring pixels (perhaps a bit less halo artifacts 
 
170
                     than IM_BICUBIC)
 
171
 *    - IM_BSPLINE - interpolates from 16 neighbouring pixels (softens image, washes colours)
 
172
 *                  (As far as I know, image should be prefiltered for this method to give 
 
173
 *                   good results... some other time :) )
 
174
 *                  This method uses bicubic interpolation kernel from CXImage 5.99a and older
 
175
 *                  versions.
 
176
 *    - IM_LANCZOS - interpolates from 12*12 pixels (slow, ringing artifacts)
 
177
 *
 
178
 *  \param ofMethod - overflow method (see comments at GetPixelColorWithOverflow)
 
179
 *  \param rplColor - pointer to color used for out of borders pixels in OM_COLOR mode
 
180
 *              (and other modes if colour can't calculated in a specified way)
 
181
 *
 
182
 *  \return interpolated color value (including interpolated alpha value, if image has alpha layer)
 
183
 * 
 
184
 *  \author ***bd*** 2.2004
 
185
 */
 
186
RGBQUAD CxImage::GetPixelColorInterpolated(
 
187
  float x,float y, 
 
188
  InterpolationMethod const inMethod, 
 
189
  OverflowMethod const ofMethod, 
 
190
  RGBQUAD* const rplColor)
 
191
{
 
192
  //calculate nearest pixel
 
193
  int xi=(int)(x); if (x<0) xi--;   //these replace (incredibly slow) floor (Visual c++ 2003, AMD Athlon)
 
194
  int yi=(int)(y); if (y<0) yi--;
 
195
  RGBQUAD color;                    //calculated colour
 
196
 
 
197
  switch (inMethod) {
 
198
    case IM_NEAREST_NEIGHBOUR:
 
199
      return GetPixelColorWithOverflow((long)(x+0.5f), (long)(y+0.5f), ofMethod, rplColor);
 
200
    default: {
 
201
      //bilinear interpolation
 
202
      if (xi<-1 || xi>=head.biWidth || yi<-1 || yi>=head.biHeight) {  //all 4 points are outside bounds?:
 
203
        switch (ofMethod) {
 
204
          case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND:
 
205
            //we don't need to interpolate anything with all points outside in this case
 
206
            return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor);
 
207
          default:
 
208
            //recalculate coordinates and use faster method later on
 
209
            OverflowCoordinates(x,y,ofMethod);
 
210
            xi=(int)(x); if (x<0) xi--;   //x and/or y have changed ... recalculate xi and yi
 
211
            yi=(int)(y); if (y<0) yi--;
 
212
        }//switch
 
213
      }//if
 
214
      //get four neighbouring pixels
 
215
      if ((xi+1)<head.biWidth && xi>=0 && (yi+1)<head.biHeight && yi>=0 && head.biClrUsed==0) {
 
216
        //all pixels are inside RGB24 image... optimize reading (and use fixed point arithmetic)
 
217
        WORD wt1=(WORD)((x-xi)*256.0f), wt2=(WORD)((y-yi)*256.0f);
 
218
        WORD wd=wt1*wt2>>8;
 
219
        WORD wb=wt1-wd;
 
220
        WORD wc=wt2-wd;
 
221
        WORD wa=256-wt1-wc;
 
222
        WORD wrr,wgg,wbb;
 
223
        BYTE *pxptr=(BYTE*)info.pImage+yi*info.dwEffWidth+xi*3;
 
224
        wbb=wa*(*pxptr++); wgg=wa*(*pxptr++); wrr=wa*(*pxptr++);
 
225
        wbb+=wb*(*pxptr++); wgg+=wb*(*pxptr++); wrr+=wb*(*pxptr);
 
226
        pxptr+=(info.dwEffWidth-5); //move to next row
 
227
        wbb+=wc*(*pxptr++); wgg+=wc*(*pxptr++); wrr+=wc*(*pxptr++); 
 
228
        wbb+=wd*(*pxptr++); wgg+=wd*(*pxptr++); wrr+=wd*(*pxptr); 
 
229
        color.rgbRed=(BYTE) (wrr>>8); color.rgbGreen=(BYTE) (wgg>>8); color.rgbBlue=(BYTE) (wbb>>8);
 
230
#if CXIMAGE_SUPPORT_ALPHA
 
231
        if (pAlpha) {
 
232
          WORD waa;
 
233
          //image has alpha layer... we have to do the same for alpha data
 
234
          pxptr=AlphaGetPointer(xi,yi);                           //pointer to first byte
 
235
          waa=wa*(*pxptr++); waa+=wb*(*pxptr);   //first two pixels
 
236
          pxptr+=(head.biWidth-1);                                //move to next row
 
237
          waa+=wc*(*pxptr++); waa+=wd*(*pxptr);   //and second row pixels
 
238
          color.rgbReserved=(BYTE) (waa>>8);
 
239
        } else
 
240
#endif
 
241
                { //Alpha not supported or no alpha at all
 
242
                        color.rgbReserved = 0;
 
243
                }
 
244
        return color;
 
245
      } else {
 
246
        //default (slower) way to get pixels (not RGB24 or some pixels out of borders)
 
247
        float t1=x-xi, t2=y-yi;
 
248
        float d=t1*t2;
 
249
        float b=t1-d;
 
250
        float c=t2-d;
 
251
        float a=1-t1-c;
 
252
        RGBQUAD rgb11,rgb21,rgb12,rgb22;
 
253
        rgb11=GetPixelColorWithOverflow(xi, yi, ofMethod, rplColor);
 
254
        rgb21=GetPixelColorWithOverflow(xi+1, yi, ofMethod, rplColor);
 
255
        rgb12=GetPixelColorWithOverflow(xi, yi+1, ofMethod, rplColor);
 
256
        rgb22=GetPixelColorWithOverflow(xi+1, yi+1, ofMethod, rplColor);
 
257
        //calculate linear interpolation
 
258
        color.rgbRed=(BYTE) (a*rgb11.rgbRed+b*rgb21.rgbRed+c*rgb12.rgbRed+d*rgb22.rgbRed);
 
259
        color.rgbGreen=(BYTE) (a*rgb11.rgbGreen+b*rgb21.rgbGreen+c*rgb12.rgbGreen+d*rgb22.rgbGreen);
 
260
        color.rgbBlue=(BYTE) (a*rgb11.rgbBlue+b*rgb21.rgbBlue+c*rgb12.rgbBlue+d*rgb22.rgbBlue);
 
261
#if CXIMAGE_SUPPORT_ALPHA
 
262
        if (AlphaIsValid())
 
263
                        color.rgbReserved=(BYTE) (a*rgb11.rgbReserved+b*rgb21.rgbReserved+c*rgb12.rgbReserved+d*rgb22.rgbReserved);
 
264
                else
 
265
#endif
 
266
                { //Alpha not supported or no alpha at all
 
267
                        color.rgbReserved = 0;
 
268
                }
 
269
        return color;
 
270
      }//if
 
271
    }//default
 
272
    case IM_BICUBIC: 
 
273
    case IM_BICUBIC2:
 
274
    case IM_BSPLINE:
 
275
        case IM_BOX:
 
276
        case IM_HERMITE:
 
277
        case IM_HAMMING:
 
278
        case IM_SINC:
 
279
        case IM_BLACKMAN:
 
280
        case IM_BESSEL:
 
281
        case IM_GAUSSIAN:
 
282
        case IM_QUADRATIC:
 
283
        case IM_MITCHELL:
 
284
        case IM_CATROM:
 
285
      //bicubic interpolation(s)
 
286
      if (((xi+2)<0) || ((xi-1)>=head.biWidth) || ((yi+2)<0) || ((yi-1)>=head.biHeight)) { //all points are outside bounds?:
 
287
        switch (ofMethod) {
 
288
          case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND:
 
289
            //we don't need to interpolate anything with all points outside in this case
 
290
            return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor);
 
291
            break;
 
292
          default:
 
293
            //recalculate coordinates and use faster method later on
 
294
            OverflowCoordinates(x,y,ofMethod);
 
295
            xi=(int)(x); if (x<0) xi--;   //x and/or y have changed ... recalculate xi and yi
 
296
            yi=(int)(y); if (y<0) yi--;
 
297
        }//switch
 
298
      }//if
 
299
 
 
300
      //some variables needed from here on
 
301
      int xii,yii;                      //x any y integer indexes for loops
 
302
      float kernel, kernelyc;           //kernel cache
 
303
      float kernelx[12], kernely[4];    //precalculated kernel values
 
304
      float rr,gg,bb,aa;                //accumulated color values
 
305
      //calculate multiplication factors for all pixels
 
306
          int i;
 
307
      switch (inMethod) {
 
308
        case IM_BICUBIC:
 
309
          for (i=0; i<4; i++) {
 
310
            kernelx[i]=KernelCubic((float)(xi+i-1-x));
 
311
            kernely[i]=KernelCubic((float)(yi+i-1-y));
 
312
          }//for i
 
313
          break;
 
314
        case IM_BICUBIC2:
 
315
          for (i=0; i<4; i++) {
 
316
            kernelx[i]=KernelGeneralizedCubic((float)(xi+i-1-x), -0.5);
 
317
            kernely[i]=KernelGeneralizedCubic((float)(yi+i-1-y), -0.5);
 
318
          }//for i
 
319
          break;
 
320
        case IM_BSPLINE:
 
321
          for (i=0; i<4; i++) {
 
322
            kernelx[i]=KernelBSpline((float)(xi+i-1-x));
 
323
            kernely[i]=KernelBSpline((float)(yi+i-1-y));
 
324
          }//for i
 
325
          break;
 
326
        case IM_BOX:
 
327
          for (i=0; i<4; i++) {
 
328
            kernelx[i]=KernelBox((float)(xi+i-1-x));
 
329
            kernely[i]=KernelBox((float)(yi+i-1-y));
 
330
          }//for i
 
331
          break;
 
332
        case IM_HERMITE:
 
333
          for (i=0; i<4; i++) {
 
334
            kernelx[i]=KernelHermite((float)(xi+i-1-x));
 
335
            kernely[i]=KernelHermite((float)(yi+i-1-y));
 
336
          }//for i
 
337
          break;
 
338
        case IM_HAMMING:
 
339
          for (i=0; i<4; i++) {
 
340
            kernelx[i]=KernelHamming((float)(xi+i-1-x));
 
341
            kernely[i]=KernelHamming((float)(yi+i-1-y));
 
342
          }//for i
 
343
          break;
 
344
        case IM_SINC:
 
345
          for (i=0; i<4; i++) {
 
346
            kernelx[i]=KernelSinc((float)(xi+i-1-x));
 
347
            kernely[i]=KernelSinc((float)(yi+i-1-y));
 
348
          }//for i
 
349
          break;
 
350
        case IM_BLACKMAN:
 
351
          for (i=0; i<4; i++) {
 
352
            kernelx[i]=KernelBlackman((float)(xi+i-1-x));
 
353
            kernely[i]=KernelBlackman((float)(yi+i-1-y));
 
354
          }//for i
 
355
          break;
 
356
        case IM_BESSEL:
 
357
          for (i=0; i<4; i++) {
 
358
            kernelx[i]=KernelBessel((float)(xi+i-1-x));
 
359
            kernely[i]=KernelBessel((float)(yi+i-1-y));
 
360
          }//for i
 
361
          break;
 
362
        case IM_GAUSSIAN:
 
363
          for (i=0; i<4; i++) {
 
364
            kernelx[i]=KernelGaussian((float)(xi+i-1-x));
 
365
            kernely[i]=KernelGaussian((float)(yi+i-1-y));
 
366
          }//for i
 
367
          break;
 
368
        case IM_QUADRATIC:
 
369
          for (i=0; i<4; i++) {
 
370
            kernelx[i]=KernelQuadratic((float)(xi+i-1-x));
 
371
            kernely[i]=KernelQuadratic((float)(yi+i-1-y));
 
372
          }//for i
 
373
          break;
 
374
        case IM_MITCHELL:
 
375
          for (i=0; i<4; i++) {
 
376
            kernelx[i]=KernelMitchell((float)(xi+i-1-x));
 
377
            kernely[i]=KernelMitchell((float)(yi+i-1-y));
 
378
          }//for i
 
379
          break;
 
380
        case IM_CATROM:
 
381
          for (i=0; i<4; i++) {
 
382
            kernelx[i]=KernelCatrom((float)(xi+i-1-x));
 
383
            kernely[i]=KernelCatrom((float)(yi+i-1-y));
 
384
          }//for i
 
385
          break;
 
386
      }//switch
 
387
      rr=gg=bb=aa=0;
 
388
      if (((xi+2)<head.biWidth) && xi>=1 && ((yi+2)<head.biHeight) && (yi>=1) && !IsIndexed()) {
 
389
        //optimized interpolation (faster pixel reads) for RGB24 images with all pixels inside bounds
 
390
        BYTE *pxptr, *pxptra;
 
391
        for (yii=yi-1; yii<yi+3; yii++) {
 
392
          pxptr=(BYTE *)BlindGetPixelPointer(xi-1, yii);    //calculate pointer to first byte in row
 
393
          kernelyc=kernely[yii-(yi-1)];
 
394
#if CXIMAGE_SUPPORT_ALPHA
 
395
          if (AlphaIsValid()) {
 
396
            //alpha is supported and valid (optimized bicubic int. for image with alpha)
 
397
            pxptra=AlphaGetPointer(xi-1, yii);
 
398
            kernel=kernelyc*kernelx[0];
 
399
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 
400
            kernel=kernelyc*kernelx[1];
 
401
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 
402
            kernel=kernelyc*kernelx[2];
 
403
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 
404
            kernel=kernelyc*kernelx[3];
 
405
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr); aa+=kernel*(*pxptra);
 
406
          } else
 
407
#endif
 
408
          //alpha not supported or valid (optimized bicubic int. for no alpha channel)
 
409
          {
 
410
            kernel=kernelyc*kernelx[0];
 
411
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 
412
            kernel=kernelyc*kernelx[1];
 
413
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 
414
            kernel=kernelyc*kernelx[2];
 
415
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 
416
            kernel=kernelyc*kernelx[3];
 
417
            bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr);
 
418
          }
 
419
        }//yii
 
420
      } else {
 
421
        //slower more flexible interpolation for border pixels and paletted images
 
422
        RGBQUAD rgbs;
 
423
        for (yii=yi-1; yii<yi+3; yii++) {
 
424
          kernelyc=kernely[yii-(yi-1)];
 
425
          for (xii=xi-1; xii<xi+3; xii++) {
 
426
            kernel=kernelyc*kernelx[xii-(xi-1)];
 
427
            rgbs=GetPixelColorWithOverflow(xii, yii, ofMethod, rplColor);
 
428
            rr+=kernel*rgbs.rgbRed;
 
429
            gg+=kernel*rgbs.rgbGreen;
 
430
            bb+=kernel*rgbs.rgbBlue;
 
431
#if CXIMAGE_SUPPORT_ALPHA
 
432
            aa+=kernel*rgbs.rgbReserved;
 
433
#endif
 
434
          }//xii
 
435
        }//yii
 
436
      }//if
 
437
      //for all colors, clip to 0..255 and assign to RGBQUAD
 
438
      if (rr>255) rr=255; if (rr<0) rr=0; color.rgbRed=(BYTE) rr;
 
439
      if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(BYTE) gg;
 
440
      if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(BYTE) bb;
 
441
#if CXIMAGE_SUPPORT_ALPHA
 
442
      if (AlphaIsValid()) {
 
443
        if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(BYTE) aa;
 
444
      } else
 
445
#endif
 
446
                { //Alpha not supported or no alpha at all
 
447
                        color.rgbReserved = 0;
 
448
                }
 
449
      return color;
 
450
    case IM_LANCZOS:
 
451
      //lanczos window (16*16) sinc interpolation
 
452
      if (((xi+6)<0) || ((xi-5)>=head.biWidth) || ((yi+6)<0) || ((yi-5)>=head.biHeight)) {
 
453
        //all points are outside bounds
 
454
        switch (ofMethod) {
 
455
          case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND:
 
456
            //we don't need to interpolate anything with all points outside in this case
 
457
            return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor);
 
458
            break;
 
459
          default:
 
460
            //recalculate coordinates and use faster method later on
 
461
            OverflowCoordinates(x,y,ofMethod);
 
462
            xi=(int)(x); if (x<0) xi--;   //x and/or y have changed ... recalculate xi and yi
 
463
            yi=(int)(y); if (y<0) yi--;
 
464
        }//switch
 
465
      }//if
 
466
 
 
467
      for (xii=xi-5; xii<xi+7; xii++) kernelx[xii-(xi-5)]=KernelLanczosSinc((float)(xii-x), 6.0f);
 
468
      rr=gg=bb=aa=0;
 
469
 
 
470
      if (((xi+6)<head.biWidth) && ((xi-5)>=0) && ((yi+6)<head.biHeight) && ((yi-5)>=0) && !IsIndexed()) {
 
471
        //optimized interpolation (faster pixel reads) for RGB24 images with all pixels inside bounds
 
472
        BYTE *pxptr, *pxptra;
 
473
        for (yii=yi-5; yii<yi+7; yii++) {
 
474
          pxptr=(BYTE *)BlindGetPixelPointer(xi-5, yii);    //calculate pointer to first byte in row
 
475
          kernelyc=KernelLanczosSinc((float)(yii-y),6.0f);
 
476
#if CXIMAGE_SUPPORT_ALPHA
 
477
          if (AlphaIsValid()) {
 
478
            //alpha is supported and valid
 
479
            pxptra=AlphaGetPointer(xi-1, yii);
 
480
            for (xii=0; xii<12; xii++) {
 
481
              kernel=kernelyc*kernelx[xii];
 
482
              bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 
483
            }//for xii
 
484
          } else
 
485
#endif
 
486
          //alpha not supported or valid
 
487
          {
 
488
            for (xii=0; xii<12; xii++) {
 
489
              kernel=kernelyc*kernelx[xii];
 
490
              bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 
491
            }//for xii
 
492
          }
 
493
        }//yii
 
494
      } else {
 
495
        //slower more flexible interpolation for border pixels and paletted images
 
496
        RGBQUAD rgbs;
 
497
        for (yii=yi-5; yii<yi+7; yii++) {
 
498
          kernelyc=KernelLanczosSinc((float)(yii-y),6.0f);
 
499
          for (xii=xi-5; xii<xi+7; xii++) {
 
500
            kernel=kernelyc*kernelx[xii-(xi-5)];
 
501
            rgbs=GetPixelColorWithOverflow(xii, yii, ofMethod, rplColor);
 
502
            rr+=kernel*rgbs.rgbRed;
 
503
            gg+=kernel*rgbs.rgbGreen;
 
504
            bb+=kernel*rgbs.rgbBlue;
 
505
#if CXIMAGE_SUPPORT_ALPHA
 
506
            aa+=kernel*rgbs.rgbReserved;
 
507
#endif
 
508
          }//xii
 
509
        }//yii
 
510
      }//if
 
511
      //for all colors, clip to 0..255 and assign to RGBQUAD
 
512
      if (rr>255) rr=255; if (rr<0) rr=0; color.rgbRed=(BYTE) rr;
 
513
      if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(BYTE) gg;
 
514
      if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(BYTE) bb;
 
515
#if CXIMAGE_SUPPORT_ALPHA
 
516
      if (AlphaIsValid()) {
 
517
        if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(BYTE) aa;   
 
518
      } else
 
519
#endif
 
520
                { //Alpha not supported or no alpha at all
 
521
                        color.rgbReserved = 0;
 
522
                }
 
523
      return color;
 
524
  }//switch
 
525
}
 
526
////////////////////////////////////////////////////////////////////////////////
 
527
/**
 
528
 * Helper function for GetAreaColorInterpolated.
 
529
 * Adds 'surf' portion of image pixel with color 'color' to (rr,gg,bb,aa).
 
530
 */
 
531
void CxImage::AddAveragingCont(RGBQUAD const &color, float const surf, float &rr, float &gg, float &bb, float &aa)
 
532
{
 
533
  rr+=color.rgbRed*surf;
 
534
  gg+=color.rgbGreen*surf;
 
535
  bb+=color.rgbBlue*surf;
 
536
#if CXIMAGE_SUPPORT_ALPHA
 
537
  aa+=color.rgbReserved*surf;
 
538
#endif
 
539
}
 
540
////////////////////////////////////////////////////////////////////////////////
 
541
/**
 
542
 * This method is similar to GetPixelColorInterpolated, but this method also properly handles 
 
543
 * subsampling.
 
544
 * If you need to sample original image with interval of more than 1 pixel (as when shrinking an image), 
 
545
 * you should use this method instead of GetPixelColorInterpolated or aliasing will occur.
 
546
 * When area width and height are both less than pixel, this method gets pixel color by interpolating
 
547
 * color of frame center with selected (inMethod) interpolation by calling GetPixelColorInterpolated. 
 
548
 * If width and height are more than 1, method calculates color by averaging color of pixels within area.
 
549
 * Interpolation method is not used in this case. Pixel color is interpolated by averaging instead.
 
550
 * If only one of both is more than 1, method uses combination of interpolation and averaging.
 
551
 * Chosen interpolation method is used, but since it is averaged later on, there is little difference
 
552
 * between IM_BILINEAR (perhaps best for this case) and better methods. IM_NEAREST_NEIGHBOUR again
 
553
 * leads to aliasing artifacts.
 
554
 * This method is a bit slower than GetPixelColorInterpolated and when aliasing is not a problem, you should
 
555
 * simply use the later. 
 
556
 *
 
557
 * \param  xc, yc - center of (rectangular) area
 
558
 * \param  w, h - width and height of area
 
559
 * \param  inMethod - interpolation method that is used, when interpolation is used (see above)
 
560
 * \param  ofMethod - overflow method used when retrieving individual pixel colors
 
561
 * \param  rplColor - replacement colour to use, in OM_COLOR
 
562
 *
 
563
 * \author ***bd*** 2.2004
 
564
 */
 
565
RGBQUAD CxImage::GetAreaColorInterpolated(
 
566
  float const xc, float const yc, float const w, float const h, 
 
567
  InterpolationMethod const inMethod, 
 
568
  OverflowMethod const ofMethod, 
 
569
  RGBQUAD* const rplColor)
 
570
{
 
571
        RGBQUAD color;      //calculated colour
 
572
        
 
573
        if (h<=1 && w<=1) {
 
574
                //both width and height are less than one... we will use interpolation of center point
 
575
                return GetPixelColorInterpolated(xc, yc, inMethod, ofMethod, rplColor);
 
576
        } else {
 
577
                //area is wider and/or taller than one pixel:
 
578
                CxRect2 area(xc-w/2.0f, yc-h/2.0f, xc+w/2.0f, yc+h/2.0f);   //area
 
579
                int xi1=(int)(area.botLeft.x+0.49999999f);                //low x
 
580
                int yi1=(int)(area.botLeft.y+0.49999999f);                //low y
 
581
                
 
582
                
 
583
                int xi2=(int)(area.topRight.x+0.5f);                      //top x
 
584
                int yi2=(int)(area.topRight.y+0.5f);                      //top y (for loops)
 
585
                
 
586
                float rr,gg,bb,aa;                                        //red, green, blue and alpha components
 
587
                rr=gg=bb=aa=0;
 
588
                int x,y;                                                  //loop counters
 
589
                float s=0;                                                //surface of all pixels
 
590
                float cps;                                                //surface of current crosssection
 
591
                if (h>1 && w>1) {
 
592
                        //width and height of area are greater than one pixel, so we can employ "ordinary" averaging
 
593
                        CxRect2 intBL, intTR;     //bottom left and top right intersection
 
594
                        intBL=area.CrossSection(CxRect2(((float)xi1)-0.5f, ((float)yi1)-0.5f, ((float)xi1)+0.5f, ((float)yi1)+0.5f));
 
595
                        intTR=area.CrossSection(CxRect2(((float)xi2)-0.5f, ((float)yi2)-0.5f, ((float)xi2)+0.5f, ((float)yi2)+0.5f));
 
596
                        float wBL, wTR, hBL, hTR;
 
597
                        wBL=intBL.Width();            //width of bottom left pixel-area intersection
 
598
                        hBL=intBL.Height();           //height of bottom left...
 
599
                        wTR=intTR.Width();            //width of top right...
 
600
                        hTR=intTR.Height();           //height of top right...
 
601
                        
 
602
                        AddAveragingCont(GetPixelColorWithOverflow(xi1,yi1,ofMethod,rplColor), wBL*hBL, rr, gg, bb, aa);    //bottom left pixel
 
603
                        AddAveragingCont(GetPixelColorWithOverflow(xi2,yi1,ofMethod,rplColor), wTR*hBL, rr, gg, bb, aa);    //bottom right pixel
 
604
                        AddAveragingCont(GetPixelColorWithOverflow(xi1,yi2,ofMethod,rplColor), wBL*hTR, rr, gg, bb, aa);    //top left pixel
 
605
                        AddAveragingCont(GetPixelColorWithOverflow(xi2,yi2,ofMethod,rplColor), wTR*hTR, rr, gg, bb, aa);    //top right pixel
 
606
                        //bottom and top row
 
607
                        for (x=xi1+1; x<xi2; x++) {
 
608
                                AddAveragingCont(GetPixelColorWithOverflow(x,yi1,ofMethod,rplColor), hBL, rr, gg, bb, aa);    //bottom row
 
609
                                AddAveragingCont(GetPixelColorWithOverflow(x,yi2,ofMethod,rplColor), hTR, rr, gg, bb, aa);    //top row
 
610
                        }
 
611
                        //leftmost and rightmost column
 
612
                        for (y=yi1+1; y<yi2; y++) {
 
613
                                AddAveragingCont(GetPixelColorWithOverflow(xi1,y,ofMethod,rplColor), wBL, rr, gg, bb, aa);    //left column
 
614
                                AddAveragingCont(GetPixelColorWithOverflow(xi2,y,ofMethod,rplColor), wTR, rr, gg, bb, aa);    //right column
 
615
                        }
 
616
                        for (y=yi1+1; y<yi2; y++) {
 
617
                                for (x=xi1+1; x<xi2; x++) { 
 
618
                                        color=GetPixelColorWithOverflow(x,y,ofMethod,rplColor);
 
619
                                        rr+=color.rgbRed;
 
620
                                        gg+=color.rgbGreen;
 
621
                                        bb+=color.rgbBlue;
 
622
#if CXIMAGE_SUPPORT_ALPHA
 
623
                                        aa+=color.rgbReserved;
 
624
#endif
 
625
                                }//for x
 
626
                        }//for y
 
627
                } else {
 
628
                        //width or height greater than one:
 
629
                        CxRect2 intersect;                                          //intersection with current pixel
 
630
                        CxPoint2 center;
 
631
                        for (y=yi1; y<=yi2; y++) {
 
632
                                for (x=xi1; x<=xi2; x++) {
 
633
                                        intersect=area.CrossSection(CxRect2(((float)x)-0.5f, ((float)y)-0.5f, ((float)x)+0.5f, ((float)y)+0.5f));
 
634
                                        center=intersect.Center();
 
635
                                        color=GetPixelColorInterpolated(center.x, center.y, inMethod, ofMethod, rplColor);
 
636
                                        cps=intersect.Surface();
 
637
                                        rr+=color.rgbRed*cps;
 
638
                                        gg+=color.rgbGreen*cps;
 
639
                                        bb+=color.rgbBlue*cps;
 
640
#if CXIMAGE_SUPPORT_ALPHA
 
641
                                        aa+=color.rgbReserved*cps;
 
642
#endif
 
643
                                }//for x
 
644
                        }//for y      
 
645
                }//if
 
646
                
 
647
                s=area.Surface();
 
648
                rr/=s; gg/=s; bb/=s; aa/=s;
 
649
                if (rr>255) rr=255; if (rr<0) rr=0; color.rgbRed=(BYTE) rr;
 
650
                if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(BYTE) gg;
 
651
                if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(BYTE) bb;
 
652
#if CXIMAGE_SUPPORT_ALPHA
 
653
                if (AlphaIsValid()) {
 
654
                        if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(BYTE) aa;
 
655
                }//if
 
656
#endif
 
657
        }//if
 
658
        return color;
 
659
}
 
660
 
 
661
////////////////////////////////////////////////////////////////////////////////
 
662
float CxImage::KernelBSpline(const float x)
 
663
{
 
664
        if (x>2.0f) return 0.0f;
 
665
        // thanks to Kristian Kratzenstein
 
666
        float a, b, c, d;
 
667
        float xm1 = x - 1.0f; // Was calculatet anyway cause the "if((x-1.0f) < 0)"
 
668
        float xp1 = x + 1.0f;
 
669
        float xp2 = x + 2.0f;
 
670
 
 
671
        if ((xp2) <= 0.0f) a = 0.0f; else a = xp2*xp2*xp2; // Only float, not float -> double -> float
 
672
        if ((xp1) <= 0.0f) b = 0.0f; else b = xp1*xp1*xp1;
 
673
        if (x <= 0) c = 0.0f; else c = x*x*x;  
 
674
        if ((xm1) <= 0.0f) d = 0.0f; else d = xm1*xm1*xm1;
 
675
 
 
676
        return (0.16666666666666666667f * (a - (4.0f * b) + (6.0f * c) - (4.0f * d)));
 
677
 
 
678
        /* equivalent <Vladim�r Kloucek>
 
679
        if (x < -2.0)
 
680
                return(0.0f);
 
681
        if (x < -1.0)
 
682
                return((2.0f+x)*(2.0f+x)*(2.0f+x)*0.16666666666666666667f);
 
683
        if (x < 0.0)
 
684
                return((4.0f+x*x*(-6.0f-3.0f*x))*0.16666666666666666667f);
 
685
        if (x < 1.0)
 
686
                return((4.0f+x*x*(-6.0f+3.0f*x))*0.16666666666666666667f);
 
687
        if (x < 2.0)
 
688
                return((2.0f-x)*(2.0f-x)*(2.0f-x)*0.16666666666666666667f);
 
689
        return(0.0f);
 
690
        */
 
691
}
 
692
 
 
693
////////////////////////////////////////////////////////////////////////////////
 
694
/**
 
695
 * Bilinear interpolation kernel:
 
696
  \verbatim
 
697
          /
 
698
         | 1-t           , if  0 <= t <= 1
 
699
  h(t) = | t+1           , if -1 <= t <  0
 
700
         | 0             , otherwise
 
701
          \
 
702
  \endverbatim
 
703
 * ***bd*** 2.2004
 
704
 */
 
705
float CxImage::KernelLinear(const float t)
 
706
{
 
707
//  if (0<=t && t<=1) return 1-t;
 
708
//  if (-1<=t && t<0) return 1+t;
 
709
//  return 0;
 
710
        
 
711
        //<Vladim�r Kloucek>
 
712
        if (t < -1.0f)
 
713
                return 0.0f;
 
714
        if (t < 0.0f)
 
715
                return 1.0f+t;
 
716
        if (t < 1.0f)
 
717
                return 1.0f-t;
 
718
        return 0.0f;
 
719
}
 
720
 
 
721
////////////////////////////////////////////////////////////////////////////////
 
722
/**
 
723
 * Bicubic interpolation kernel (a=-1):
 
724
  \verbatim
 
725
          /
 
726
         | 1-2|t|**2+|t|**3          , if |t| < 1
 
727
  h(t) = | 4-8|t|+5|t|**2-|t|**3     , if 1<=|t|<2
 
728
         | 0                         , otherwise
 
729
          \
 
730
  \endverbatim
 
731
 * ***bd*** 2.2004
 
732
 */
 
733
float CxImage::KernelCubic(const float t)
 
734
{
 
735
  float abs_t = (float)fabs(t);
 
736
  float abs_t_sq = abs_t * abs_t;
 
737
  if (abs_t<1) return 1-2*abs_t_sq+abs_t_sq*abs_t;
 
738
  if (abs_t<2) return 4 - 8*abs_t +5*abs_t_sq - abs_t_sq*abs_t;
 
739
  return 0;
 
740
}
 
741
 
 
742
////////////////////////////////////////////////////////////////////////////////
 
743
/**
 
744
 * Bicubic kernel (for a=-1 it is the same as BicubicKernel):
 
745
  \verbatim
 
746
          /
 
747
         | (a+2)|t|**3 - (a+3)|t|**2 + 1     , |t| <= 1
 
748
  h(t) = | a|t|**3 - 5a|t|**2 + 8a|t| - 4a   , 1 < |t| <= 2
 
749
         | 0                                 , otherwise
 
750
          \
 
751
  \endverbatim
 
752
 * Often used values for a are -1 and -1/2.
 
753
 */
 
754
float CxImage::KernelGeneralizedCubic(const float t, const float a)
 
755
{
 
756
  float abs_t = (float)fabs(t);
 
757
  float abs_t_sq = abs_t * abs_t;
 
758
  if (abs_t<1) return (a+2)*abs_t_sq*abs_t - (a+3)*abs_t_sq + 1;
 
759
  if (abs_t<2) return a*abs_t_sq*abs_t - 5*a*abs_t_sq + 8*a*abs_t - 4*a;
 
760
  return 0;
 
761
}
 
762
 
 
763
////////////////////////////////////////////////////////////////////////////////
 
764
/**
 
765
 * Lanczos windowed sinc interpolation kernel with radius r.
 
766
  \verbatim
 
767
          /
 
768
  h(t) = | sinc(t)*sinc(t/r)       , if |t|<r
 
769
         | 0                       , otherwise
 
770
          \
 
771
  \endverbatim
 
772
 * ***bd*** 2.2004
 
773
 */
 
774
float CxImage::KernelLanczosSinc(const float t, const float r)
 
775
{
 
776
  if (fabs(t) > r) return 0;
 
777
  if (t==0) return 1;
 
778
  float pit=PI*t;
 
779
  float pitd=pit/r;
 
780
  return (float)((sin(pit)/pit) * (sin(pitd)/pitd));
 
781
}
 
782
 
 
783
////////////////////////////////////////////////////////////////////////////////
 
784
float CxImage::KernelBox(const float x)
 
785
{
 
786
        if (x < -0.5f)
 
787
                return 0.0f;
 
788
        if (x < 0.5f)
 
789
                return 1.0f;
 
790
        return 0.0f;
 
791
}
 
792
////////////////////////////////////////////////////////////////////////////////
 
793
float CxImage::KernelHermite(const float x)
 
794
{
 
795
        if (x < -1.0f)
 
796
                return 0.0f;
 
797
        if (x < 0.0f)
 
798
                return (-2.0f*x-3.0f)*x*x+1.0f;
 
799
        if (x < 1.0f)
 
800
                return (2.0f*x-3.0f)*x*x+1.0f;
 
801
        return 0.0f;
 
802
//      if (fabs(x)>1) return 0.0f;
 
803
//      return(0.5f+0.5f*(float)cos(PI*x));
 
804
}
 
805
////////////////////////////////////////////////////////////////////////////////
 
806
float CxImage::KernelHamming(const float x)
 
807
{
 
808
        if (x < -1.0f)
 
809
                return 0.0f;
 
810
        if (x < 0.0f)
 
811
                return 0.92f*(-2.0f*x-3.0f)*x*x+1.0f;
 
812
        if (x < 1.0f)
 
813
                return 0.92f*(2.0f*x-3.0f)*x*x+1.0f;
 
814
        return 0.0f;
 
815
//      if (fabs(x)>1) return 0.0f;
 
816
//      return(0.54f+0.46f*(float)cos(PI*x));
 
817
}
 
818
////////////////////////////////////////////////////////////////////////////////
 
819
float CxImage::KernelSinc(const float x)
 
820
{
 
821
        if (x == 0.0)
 
822
                return(1.0);
 
823
        return((float)sin(PI*x)/(PI*x));
 
824
}
 
825
////////////////////////////////////////////////////////////////////////////////
 
826
float CxImage::KernelBlackman(const float x)
 
827
{
 
828
        //if (fabs(x)>1) return 0.0f;
 
829
        return (0.42f+0.5f*(float)cos(PI*x)+0.08f*(float)cos(2.0f*PI*x));
 
830
}
 
831
////////////////////////////////////////////////////////////////////////////////
 
832
float CxImage::KernelBessel_J1(const float x)
 
833
{
 
834
        double p, q;
 
835
        
 
836
        register long i;
 
837
        
 
838
        static const double
 
839
        Pone[] =
 
840
        {
 
841
                0.581199354001606143928050809e+21,
 
842
                -0.6672106568924916298020941484e+20,
 
843
                0.2316433580634002297931815435e+19,
 
844
                -0.3588817569910106050743641413e+17,
 
845
                0.2908795263834775409737601689e+15,
 
846
                -0.1322983480332126453125473247e+13,
 
847
                0.3413234182301700539091292655e+10,
 
848
                -0.4695753530642995859767162166e+7,
 
849
                0.270112271089232341485679099e+4
 
850
        },
 
851
        Qone[] =
 
852
        {
 
853
                0.11623987080032122878585294e+22,
 
854
                0.1185770712190320999837113348e+20,
 
855
                0.6092061398917521746105196863e+17,
 
856
                0.2081661221307607351240184229e+15,
 
857
                0.5243710262167649715406728642e+12,
 
858
                0.1013863514358673989967045588e+10,
 
859
                0.1501793594998585505921097578e+7,
 
860
                0.1606931573481487801970916749e+4,
 
861
                0.1e+1
 
862
        };
 
863
                
 
864
        p = Pone[8];
 
865
        q = Qone[8];
 
866
        for (i=7; i >= 0; i--)
 
867
        {
 
868
                p = p*x*x+Pone[i];
 
869
                q = q*x*x+Qone[i];
 
870
        }
 
871
        return (float)(p/q);
 
872
}
 
873
////////////////////////////////////////////////////////////////////////////////
 
874
float CxImage::KernelBessel_P1(const float x)
 
875
{
 
876
        double p, q;
 
877
        
 
878
        register long i;
 
879
        
 
880
        static const double
 
881
        Pone[] =
 
882
        {
 
883
                0.352246649133679798341724373e+5,
 
884
                0.62758845247161281269005675e+5,
 
885
                0.313539631109159574238669888e+5,
 
886
                0.49854832060594338434500455e+4,
 
887
                0.2111529182853962382105718e+3,
 
888
                0.12571716929145341558495e+1
 
889
        },
 
890
        Qone[] =
 
891
        {
 
892
                0.352246649133679798068390431e+5,
 
893
                0.626943469593560511888833731e+5,
 
894
                0.312404063819041039923015703e+5,
 
895
                0.4930396490181088979386097e+4,
 
896
                0.2030775189134759322293574e+3,
 
897
                0.1e+1
 
898
        };
 
899
                
 
900
        p = Pone[5];
 
901
        q = Qone[5];
 
902
        for (i=4; i >= 0; i--)
 
903
        {
 
904
                p = p*(8.0/x)*(8.0/x)+Pone[i];
 
905
                q = q*(8.0/x)*(8.0/x)+Qone[i];
 
906
        }
 
907
        return (float)(p/q);
 
908
}
 
909
////////////////////////////////////////////////////////////////////////////////
 
910
float CxImage::KernelBessel_Q1(const float x)
 
911
{
 
912
        double p, q;
 
913
        
 
914
        register long i;
 
915
        
 
916
        static const double
 
917
        Pone[] =
 
918
        {
 
919
                0.3511751914303552822533318e+3,
 
920
                0.7210391804904475039280863e+3,
 
921
                0.4259873011654442389886993e+3,
 
922
                0.831898957673850827325226e+2,
 
923
                0.45681716295512267064405e+1,
 
924
                0.3532840052740123642735e-1
 
925
        },
 
926
        Qone[] =
 
927
        {
 
928
                0.74917374171809127714519505e+4,
 
929
                0.154141773392650970499848051e+5,
 
930
                0.91522317015169922705904727e+4,
 
931
                0.18111867005523513506724158e+4,
 
932
                0.1038187585462133728776636e+3,
 
933
                0.1e+1
 
934
        };
 
935
                
 
936
        p = Pone[5];
 
937
        q = Qone[5];
 
938
        for (i=4; i >= 0; i--)
 
939
        {
 
940
                p = p*(8.0/x)*(8.0/x)+Pone[i];
 
941
                q = q*(8.0/x)*(8.0/x)+Qone[i];
 
942
        }
 
943
        return (float)(p/q);
 
944
}
 
945
////////////////////////////////////////////////////////////////////////////////
 
946
float CxImage::KernelBessel_Order1(float x)
 
947
{
 
948
        float p, q;
 
949
        
 
950
        if (x == 0.0)
 
951
                return (0.0f);
 
952
        p = x;
 
953
        if (x < 0.0)
 
954
                x=(-x);
 
955
        if (x < 8.0)
 
956
                return(p*KernelBessel_J1(x));
 
957
        q = (float)sqrt(2.0f/(PI*x))*(float)(KernelBessel_P1(x)*(1.0f/sqrt(2.0f)*(sin(x)-cos(x)))-8.0f/x*KernelBessel_Q1(x)*
 
958
                (-1.0f/sqrt(2.0f)*(sin(x)+cos(x))));
 
959
        if (p < 0.0f)
 
960
                q = (-q);
 
961
        return (q);
 
962
}
 
963
////////////////////////////////////////////////////////////////////////////////
 
964
float CxImage::KernelBessel(const float x)
 
965
{
 
966
        if (x == 0.0f)
 
967
                return(PI/4.0f);
 
968
        return(KernelBessel_Order1(PI*x)/(2.0f*x));
 
969
}
 
970
////////////////////////////////////////////////////////////////////////////////
 
971
float CxImage::KernelGaussian(const float x)
 
972
{
 
973
        return (float)(exp(-2.0f*x*x)*0.79788456080287f/*sqrt(2.0f/PI)*/);
 
974
}
 
975
////////////////////////////////////////////////////////////////////////////////
 
976
float CxImage::KernelQuadratic(const float x)
 
977
{
 
978
        if (x < -1.5f)
 
979
                return(0.0f);
 
980
        if (x < -0.5f)
 
981
                return(0.5f*(x+1.5f)*(x+1.5f));
 
982
        if (x < 0.5f)
 
983
                return(0.75f-x*x);
 
984
        if (x < 1.5f)
 
985
                return(0.5f*(x-1.5f)*(x-1.5f));
 
986
        return(0.0f);
 
987
}
 
988
////////////////////////////////////////////////////////////////////////////////
 
989
float CxImage::KernelMitchell(const float x)
 
990
{
 
991
#define KM_B (1.0f/3.0f)
 
992
#define KM_C (1.0f/3.0f)
 
993
#define KM_P0 ((  6.0f - 2.0f * KM_B ) / 6.0f)
 
994
#define KM_P2 ((-18.0f + 12.0f * KM_B + 6.0f * KM_C) / 6.0f)
 
995
#define KM_P3 (( 12.0f - 9.0f  * KM_B - 6.0f * KM_C) / 6.0f)
 
996
#define KM_Q0 ((  8.0f * KM_B + 24.0f * KM_C) / 6.0f)
 
997
#define KM_Q1 ((-12.0f * KM_B - 48.0f * KM_C) / 6.0f)
 
998
#define KM_Q2 ((  6.0f * KM_B + 30.0f * KM_C) / 6.0f)
 
999
#define KM_Q3 (( -1.0f * KM_B -  6.0f * KM_C) / 6.0f)
 
1000
        
 
1001
        if (x < -2.0)
 
1002
                return(0.0f);
 
1003
        if (x < -1.0)
 
1004
                return(KM_Q0-x*(KM_Q1-x*(KM_Q2-x*KM_Q3)));
 
1005
        if (x < 0.0f)
 
1006
                return(KM_P0+x*x*(KM_P2-x*KM_P3));
 
1007
        if (x < 1.0f)
 
1008
                return(KM_P0+x*x*(KM_P2+x*KM_P3));
 
1009
        if (x < 2.0f)
 
1010
                return(KM_Q0+x*(KM_Q1+x*(KM_Q2+x*KM_Q3)));
 
1011
        return(0.0f);
 
1012
}
 
1013
////////////////////////////////////////////////////////////////////////////////
 
1014
float CxImage::KernelCatrom(const float x)
 
1015
{
 
1016
        if (x < -2.0)
 
1017
                return(0.0f);
 
1018
        if (x < -1.0)
 
1019
                return(0.5f*(4.0f+x*(8.0f+x*(5.0f+x))));
 
1020
        if (x < 0.0)
 
1021
                return(0.5f*(2.0f+x*x*(-5.0f-3.0f*x)));
 
1022
        if (x < 1.0)
 
1023
                return(0.5f*(2.0f+x*x*(-5.0f+3.0f*x)));
 
1024
        if (x < 2.0)
 
1025
                return(0.5f*(4.0f+x*(-8.0f+x*(5.0f-x))));
 
1026
        return(0.0f);
 
1027
}
 
1028
////////////////////////////////////////////////////////////////////////////////
 
1029
 
 
1030
#endif