utilities/CxImage/ximaint.cpp

   1 // xImaInt.cpp : interpolation functions
   2 /* 02/2004 - Branko Brevensek
   3  * CxImage version 6.0.0 02/Feb/2008 - 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 >= 0) {
 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       //IM_BILINEAR: 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         case IM_HANNING:
 286         case IM_POWER:
 287       //bicubic interpolation(s)
 288       if (((xi+2)<0) || ((xi-1)>=head.biWidth) || ((yi+2)<0) || ((yi-1)>=head.biHeight)) { //all points are outside bounds?:
 289         switch (ofMethod) {
 290           case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND:
 291             //we don't need to interpolate anything with all points outside in this case
 292             return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor);
 293             break;
 294           default:
 295             //recalculate coordinates and use faster method later on
 296             OverflowCoordinates(x,y,ofMethod);
 297             xi=(int)(x); if (x<0) xi--;   //x and/or y have changed ... recalculate xi and yi
 298             yi=(int)(y); if (y<0) yi--;
 299         }//switch
 300       }//if
 301
 302       //some variables needed from here on
 303       int xii,yii;                      //x any y integer indexes for loops
 304       float kernel, kernelyc;           //kernel cache
 305       float kernelx[12], kernely[4];    //precalculated kernel values
 306       float rr,gg,bb,aa;                //accumulated color values
 307       //calculate multiplication factors for all pixels
 308           int i;
 309       switch (inMethod) {
 310         case IM_BICUBIC:
 311           for (i=0; i<4; i++) {
 312             kernelx[i]=KernelCubic((float)(xi+i-1-x));
 313             kernely[i]=KernelCubic((float)(yi+i-1-y));
 314           }//for i
 315           break;
 316         case IM_BICUBIC2:
 317           for (i=0; i<4; i++) {
 318             kernelx[i]=KernelGeneralizedCubic((float)(xi+i-1-x), -0.5);
 319             kernely[i]=KernelGeneralizedCubic((float)(yi+i-1-y), -0.5);
 320           }//for i
 321           break;
 322         case IM_BSPLINE:
 323           for (i=0; i<4; i++) {
 324             kernelx[i]=KernelBSpline((float)(xi+i-1-x));
 325             kernely[i]=KernelBSpline((float)(yi+i-1-y));
 326           }//for i
 327           break;
 328         case IM_BOX:
 329           for (i=0; i<4; i++) {
 330             kernelx[i]=KernelBox((float)(xi+i-1-x));
 331             kernely[i]=KernelBox((float)(yi+i-1-y));
 332           }//for i
 333           break;
 334         case IM_HERMITE:
 335           for (i=0; i<4; i++) {
 336             kernelx[i]=KernelHermite((float)(xi+i-1-x));
 337             kernely[i]=KernelHermite((float)(yi+i-1-y));
 338           }//for i
 339           break;
 340         case IM_HAMMING:
 341           for (i=0; i<4; i++) {
 342             kernelx[i]=KernelHamming((float)(xi+i-1-x));
 343             kernely[i]=KernelHamming((float)(yi+i-1-y));
 344           }//for i
 345           break;
 346         case IM_SINC:
 347           for (i=0; i<4; i++) {
 348             kernelx[i]=KernelSinc((float)(xi+i-1-x));
 349             kernely[i]=KernelSinc((float)(yi+i-1-y));
 350           }//for i
 351           break;
 352         case IM_BLACKMAN:
 353           for (i=0; i<4; i++) {
 354             kernelx[i]=KernelBlackman((float)(xi+i-1-x));
 355             kernely[i]=KernelBlackman((float)(yi+i-1-y));
 356           }//for i
 357           break;
 358         case IM_BESSEL:
 359           for (i=0; i<4; i++) {
 360             kernelx[i]=KernelBessel((float)(xi+i-1-x));
 361             kernely[i]=KernelBessel((float)(yi+i-1-y));
 362           }//for i
 363           break;
 364         case IM_GAUSSIAN:
 365           for (i=0; i<4; i++) {
 366             kernelx[i]=KernelGaussian((float)(xi+i-1-x));
 367             kernely[i]=KernelGaussian((float)(yi+i-1-y));
 368           }//for i
 369           break;
 370         case IM_QUADRATIC:
 371           for (i=0; i<4; i++) {
 372             kernelx[i]=KernelQuadratic((float)(xi+i-1-x));
 373             kernely[i]=KernelQuadratic((float)(yi+i-1-y));
 374           }//for i
 375           break;
 376         case IM_MITCHELL:
 377           for (i=0; i<4; i++) {
 378             kernelx[i]=KernelMitchell((float)(xi+i-1-x));
 379             kernely[i]=KernelMitchell((float)(yi+i-1-y));
 380           }//for i
 381           break;
 382         case IM_CATROM:
 383           for (i=0; i<4; i++) {
 384             kernelx[i]=KernelCatrom((float)(xi+i-1-x));
 385             kernely[i]=KernelCatrom((float)(yi+i-1-y));
 386           }//for i
 387           break;
 388         case IM_HANNING:
 389           for (i=0; i<4; i++) {
 390             kernelx[i]=KernelHanning((float)(xi+i-1-x));
 391             kernely[i]=KernelHanning((float)(yi+i-1-y));
 392           }//for i
 393           break;
 394         case IM_POWER:
 395           for (i=0; i<4; i++) {
 396             kernelx[i]=KernelPower((float)(xi+i-1-x));
 397             kernely[i]=KernelPower((float)(yi+i-1-y));
 398           }//for i
 399           break;
 400         default:
 401           i=0;
 402       }//switch
 403       rr=gg=bb=aa=0;
 404       if (((xi+2)<head.biWidth) && xi>=1 && ((yi+2)<head.biHeight) && (yi>=1) && !IsIndexed()) {
 405         //optimized interpolation (faster pixel reads) for RGB24 images with all pixels inside bounds
 406         BYTE *pxptr, *pxptra;
 407         for (yii=yi-1; yii<yi+3; yii++) {
 408           pxptr=(BYTE *)BlindGetPixelPointer(xi-1, yii);    //calculate pointer to first byte in row
 409           kernelyc=kernely[yii-(yi-1)];
 410 #if CXIMAGE_SUPPORT_ALPHA
 411           if (AlphaIsValid()) {
 412             //alpha is supported and valid (optimized bicubic int. for image with alpha)
 413             pxptra=AlphaGetPointer(xi-1, yii);
 414             kernel=kernelyc*kernelx[0];
 415             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 416             kernel=kernelyc*kernelx[1];
 417             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 418             kernel=kernelyc*kernelx[2];
 419             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 420             kernel=kernelyc*kernelx[3];
 421             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr); aa+=kernel*(*pxptra);
 422           } else
 423 #endif
 424           //alpha not supported or valid (optimized bicubic int. for no alpha channel)
 425           {
 426             kernel=kernelyc*kernelx[0];
 427             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 428             kernel=kernelyc*kernelx[1];
 429             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 430             kernel=kernelyc*kernelx[2];
 431             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 432             kernel=kernelyc*kernelx[3];
 433             bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr);
 434           }
 435         }//yii
 436       } else {
 437         //slower more flexible interpolation for border pixels and paletted images
 438         RGBQUAD rgbs;
 439         for (yii=yi-1; yii<yi+3; yii++) {
 440           kernelyc=kernely[yii-(yi-1)];
 441           for (xii=xi-1; xii<xi+3; xii++) {
 442             kernel=kernelyc*kernelx[xii-(xi-1)];
 443             rgbs=GetPixelColorWithOverflow(xii, yii, ofMethod, rplColor);
 444             rr+=kernel*rgbs.rgbRed;
 445             gg+=kernel*rgbs.rgbGreen;
 446             bb+=kernel*rgbs.rgbBlue;
 447 #if CXIMAGE_SUPPORT_ALPHA
 448             aa+=kernel*rgbs.rgbReserved;
 449 #endif
 450           }//xii
 451         }//yii
 452       }//if
 453       //for all colors, clip to 0..255 and assign to RGBQUAD
 454       if (rr>255) rr=255; if (rr<0) rr=0; color.rgbRed=(BYTE) rr;
 455       if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(BYTE) gg;
 456       if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(BYTE) bb;
 457 #if CXIMAGE_SUPPORT_ALPHA
 458       if (AlphaIsValid()) {
 459         if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(BYTE) aa;
 460       } else
 461 #endif
 462                 { //Alpha not supported or no alpha at all
 463                         color.rgbReserved = 0;
 464                 }
 465       return color;
 466     case IM_LANCZOS:
 467       //lanczos window (16*16) sinc interpolation
 468       if (((xi+6)<0) || ((xi-5)>=head.biWidth) || ((yi+6)<0) || ((yi-5)>=head.biHeight)) {
 469         //all points are outside bounds
 470         switch (ofMethod) {
 471           case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND:
 472             //we don't need to interpolate anything with all points outside in this case
 473             return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor);
 474             break;
 475           default:
 476             //recalculate coordinates and use faster method later on
 477             OverflowCoordinates(x,y,ofMethod);
 478             xi=(int)(x); if (x<0) xi--;   //x and/or y have changed ... recalculate xi and yi
 479             yi=(int)(y); if (y<0) yi--;
 480         }//switch
 481       }//if
 482
 483       for (xii=xi-5; xii<xi+7; xii++) kernelx[xii-(xi-5)]=KernelLanczosSinc((float)(xii-x), 6.0f);
 484       rr=gg=bb=aa=0;
 485
 486       if (((xi+6)<head.biWidth) && ((xi-5)>=0) && ((yi+6)<head.biHeight) && ((yi-5)>=0) && !IsIndexed()) {
 487         //optimized interpolation (faster pixel reads) for RGB24 images with all pixels inside bounds
 488         BYTE *pxptr, *pxptra;
 489         for (yii=yi-5; yii<yi+7; yii++) {
 490           pxptr=(BYTE *)BlindGetPixelPointer(xi-5, yii);    //calculate pointer to first byte in row
 491           kernelyc=KernelLanczosSinc((float)(yii-y),6.0f);
 492 #if CXIMAGE_SUPPORT_ALPHA
 493           if (AlphaIsValid()) {
 494             //alpha is supported and valid
 495             pxptra=AlphaGetPointer(xi-1, yii);
 496             for (xii=0; xii<12; xii++) {
 497               kernel=kernelyc*kernelx[xii];
 498               bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++); aa+=kernel*(*pxptra++);
 499             }//for xii
 500           } else
 501 #endif
 502           //alpha not supported or valid
 503           {
 504             for (xii=0; xii<12; xii++) {
 505               kernel=kernelyc*kernelx[xii];
 506               bb+=kernel*(*pxptr++); gg+=kernel*(*pxptr++); rr+=kernel*(*pxptr++);
 507             }//for xii
 508           }
 509         }//yii
 510       } else {
 511         //slower more flexible interpolation for border pixels and paletted images
 512         RGBQUAD rgbs;
 513         for (yii=yi-5; yii<yi+7; yii++) {
 514           kernelyc=KernelLanczosSinc((float)(yii-y),6.0f);
 515           for (xii=xi-5; xii<xi+7; xii++) {
 516             kernel=kernelyc*kernelx[xii-(xi-5)];
 517             rgbs=GetPixelColorWithOverflow(xii, yii, ofMethod, rplColor);
 518             rr+=kernel*rgbs.rgbRed;
 519             gg+=kernel*rgbs.rgbGreen;
 520             bb+=kernel*rgbs.rgbBlue;
 521 #if CXIMAGE_SUPPORT_ALPHA
 522             aa+=kernel*rgbs.rgbReserved;
 523 #endif
 524           }//xii
 525         }//yii
 526       }//if
 527       //for all colors, clip to 0..255 and assign to RGBQUAD
 528       if (rr>255) rr=255; if (rr<0) rr=0; color.rgbRed=(BYTE) rr;
 529       if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(BYTE) gg;
 530       if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(BYTE) bb;
 531 #if CXIMAGE_SUPPORT_ALPHA
 532       if (AlphaIsValid()) {
 533         if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(BYTE) aa;
 534       } else
 535 #endif
 536                 { //Alpha not supported or no alpha at all
 537                         color.rgbReserved = 0;
 538                 }
 539       return color;
 540   }//switch
 541 }
 542 ////////////////////////////////////////////////////////////////////////////////
 543 /**
 544  * Helper function for GetAreaColorInterpolated.
 545  * Adds 'surf' portion of image pixel with color 'color' to (rr,gg,bb,aa).
 546  */
 547 void CxImage::AddAveragingCont(RGBQUAD const &color, float const surf, float &rr, float &gg, float &bb, float &aa)
 548 {
 549   rr+=color.rgbRed*surf;
 550   gg+=color.rgbGreen*surf;
 551   bb+=color.rgbBlue*surf;
 552 #if CXIMAGE_SUPPORT_ALPHA
 553   aa+=color.rgbReserved*surf;
 554 #endif
 555 }
 556 ////////////////////////////////////////////////////////////////////////////////
 557 /**
 558  * This method is similar to GetPixelColorInterpolated, but this method also properly handles
 559  * subsampling.
 560  * If you need to sample original image with interval of more than 1 pixel (as when shrinking an image),
 561  * you should use this method instead of GetPixelColorInterpolated or aliasing will occur.
 562  * When area width and height are both less than pixel, this method gets pixel color by interpolating
 563  * color of frame center with selected (inMethod) interpolation by calling GetPixelColorInterpolated.
 564  * If width and height are more than 1, method calculates color by averaging color of pixels within area.
 565  * Interpolation method is not used in this case. Pixel color is interpolated by averaging instead.
 566  * If only one of both is more than 1, method uses combination of interpolation and averaging.
 567  * Chosen interpolation method is used, but since it is averaged later on, there is little difference
 568  * between IM_BILINEAR (perhaps best for this case) and better methods. IM_NEAREST_NEIGHBOUR again
 569  * leads to aliasing artifacts.
 570  * This method is a bit slower than GetPixelColorInterpolated and when aliasing is not a problem, you should
 571  * simply use the later.
 572  *
 573  * \param  xc, yc - center of (rectangular) area
 574  * \param  w, h - width and height of area
 575  * \param  inMethod - interpolation method that is used, when interpolation is used (see above)
 576  * \param  ofMethod - overflow method used when retrieving individual pixel colors
 577  * \param  rplColor - replacement colour to use, in OM_COLOR
 578  *
 579  * \author ***bd*** 2.2004
 580  */
 581 RGBQUAD CxImage::GetAreaColorInterpolated(
 582   float const xc, float const yc, float const w, float const h,
 583   InterpolationMethod const inMethod,
 584   OverflowMethod const ofMethod,
 585   RGBQUAD* const rplColor)
 586 {
 587         RGBQUAD color;      //calculated colour
 588
 589         if (h<=1 && w<=1) {
 590                 //both width and height are less than one... we will use interpolation of center point
 591                 return GetPixelColorInterpolated(xc, yc, inMethod, ofMethod, rplColor);
 592         } else {
 593                 //area is wider and/or taller than one pixel:
 594                 CxRect2 area(xc-w/2.0f, yc-h/2.0f, xc+w/2.0f, yc+h/2.0f);   //area
 595                 int xi1=(int)(area.botLeft.x+0.49999999f);                //low x
 596                 int yi1=(int)(area.botLeft.y+0.49999999f);                //low y
 597
 598
 599                 int xi2=(int)(area.topRight.x+0.5f);                      //top x
 600                 int yi2=(int)(area.topRight.y+0.5f);                      //top y (for loops)
 601
 602                 float rr,gg,bb,aa;                                        //red, green, blue and alpha components
 603                 rr=gg=bb=aa=0;
 604                 int x,y;                                                  //loop counters
 605                 float s=0;                                                //surface of all pixels
 606                 float cps;                                                //surface of current crosssection
 607                 if (h>1 && w>1) {
 608                         //width and height of area are greater than one pixel, so we can employ "ordinary" averaging
 609                         CxRect2 intBL, intTR;     //bottom left and top right intersection
 610                         intBL=area.CrossSection(CxRect2(((float)xi1)-0.5f, ((float)yi1)-0.5f, ((float)xi1)+0.5f, ((float)yi1)+0.5f));
 611                         intTR=area.CrossSection(CxRect2(((float)xi2)-0.5f, ((float)yi2)-0.5f, ((float)xi2)+0.5f, ((float)yi2)+0.5f));
 612                         float wBL, wTR, hBL, hTR;
 613                         wBL=intBL.Width();            //width of bottom left pixel-area intersection
 614                         hBL=intBL.Height();           //height of bottom left...
 615                         wTR=intTR.Width();            //width of top right...
 616                         hTR=intTR.Height();           //height of top right...
 617
 618                         AddAveragingCont(GetPixelColorWithOverflow(xi1,yi1,ofMethod,rplColor), wBL*hBL, rr, gg, bb, aa);    //bottom left pixel
 619                         AddAveragingCont(GetPixelColorWithOverflow(xi2,yi1,ofMethod,rplColor), wTR*hBL, rr, gg, bb, aa);    //bottom right pixel
 620                         AddAveragingCont(GetPixelColorWithOverflow(xi1,yi2,ofMethod,rplColor), wBL*hTR, rr, gg, bb, aa);    //top left pixel
 621                         AddAveragingCont(GetPixelColorWithOverflow(xi2,yi2,ofMethod,rplColor), wTR*hTR, rr, gg, bb, aa);    //top right pixel
 622                         //bottom and top row
 623                         for (x=xi1+1; x<xi2; x++) {
 624                                 AddAveragingCont(GetPixelColorWithOverflow(x,yi1,ofMethod,rplColor), hBL, rr, gg, bb, aa);    //bottom row
 625                                 AddAveragingCont(GetPixelColorWithOverflow(x,yi2,ofMethod,rplColor), hTR, rr, gg, bb, aa);    //top row
 626                         }
 627                         //leftmost and rightmost column
 628                         for (y=yi1+1; y<yi2; y++) {
 629                                 AddAveragingCont(GetPixelColorWithOverflow(xi1,y,ofMethod,rplColor), wBL, rr, gg, bb, aa);    //left column
 630                                 AddAveragingCont(GetPixelColorWithOverflow(xi2,y,ofMethod,rplColor), wTR, rr, gg, bb, aa);    //right column
 631                         }
 632                         for (y=yi1+1; y<yi2; y++) {
 633                                 for (x=xi1+1; x<xi2; x++) {
 634                                         color=GetPixelColorWithOverflow(x,y,ofMethod,rplColor);
 635                                         rr+=color.rgbRed;
 636                                         gg+=color.rgbGreen;
 637                                         bb+=color.rgbBlue;
 638 #if CXIMAGE_SUPPORT_ALPHA
 639                                         aa+=color.rgbReserved;
 640 #endif
 641                                 }//for x
 642                         }//for y
 643                 } else {
 644                         //width or height greater than one:
 645                         CxRect2 intersect;                                          //intersection with current pixel
 646                         CxPoint2 center;
 647                         for (y=yi1; y<=yi2; y++) {
 648                                 for (x=xi1; x<=xi2; x++) {
 649                                         intersect=area.CrossSection(CxRect2(((float)x)-0.5f, ((float)y)-0.5f, ((float)x)+0.5f, ((float)y)+0.5f));
 650                                         center=intersect.Center();
 651                                         color=GetPixelColorInterpolated(center.x, center.y, inMethod, ofMethod, rplColor);
 652                                         cps=intersect.Surface();
 653                                         rr+=color.rgbRed*cps;
 654                                         gg+=color.rgbGreen*cps;
 655                                         bb+=color.rgbBlue*cps;
 656 #if CXIMAGE_SUPPORT_ALPHA
 657                                         aa+=color.rgbReserved*cps;
 658 #endif
 659                                 }//for x
 660                         }//for y
 661                 }//if
 662
 663                 s=area.Surface();
 664                 rr/=s; gg/=s; bb/=s; aa/=s;
 665                 if (rr>255) rr=255; if (rr<0) rr=0; color.rgbRed=(BYTE) rr;
 666                 if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(BYTE) gg;
 667                 if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(BYTE) bb;
 668 #if CXIMAGE_SUPPORT_ALPHA
 669                 if (AlphaIsValid()) {
 670                         if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(BYTE) aa;
 671                 }//if
 672 #endif
 673         }//if
 674         return color;
 675 }
 676
 677 ////////////////////////////////////////////////////////////////////////////////
 678 float CxImage::KernelBSpline(const float x)
 679 {
 680         if (x>2.0f) return 0.0f;
 681         // thanks to Kristian Kratzenstein
 682         float a, b, c, d;
 683         float xm1 = x - 1.0f; // Was calculatet anyway cause the "if((x-1.0f) < 0)"
 684         float xp1 = x + 1.0f;
 685         float xp2 = x + 2.0f;
 686
 687         if ((xp2) <= 0.0f) a = 0.0f; else a = xp2*xp2*xp2; // Only float, not float -> double -> float
 688         if ((xp1) <= 0.0f) b = 0.0f; else b = xp1*xp1*xp1;
 689         if (x <= 0) c = 0.0f; else c = x*x*x;
 690         if ((xm1) <= 0.0f) d = 0.0f; else d = xm1*xm1*xm1;
 691
 692         return (0.16666666666666666667f * (a - (4.0f * b) + (6.0f * c) - (4.0f * d)));
 693
 694         /* equivalent <Vladimir Kloucek>
 695         if (x < -2.0)
 696                 return(0.0f);
 697         if (x < -1.0)
 698                 return((2.0f+x)*(2.0f+x)*(2.0f+x)*0.16666666666666666667f);
 699         if (x < 0.0)
 700                 return((4.0f+x*x*(-6.0f-3.0f*x))*0.16666666666666666667f);
 701         if (x < 1.0)
 702                 return((4.0f+x*x*(-6.0f+3.0f*x))*0.16666666666666666667f);
 703         if (x < 2.0)
 704                 return((2.0f-x)*(2.0f-x)*(2.0f-x)*0.16666666666666666667f);
 705         return(0.0f);
 706         */
 707 }
 708
 709 ////////////////////////////////////////////////////////////////////////////////
 710 /**
 711  * Bilinear interpolation kernel:
 712   \verbatim
 713           /
 714          | 1-t           , if  0 <= t <= 1
 715   h(t) = | t+1           , if -1 <= t <  0
 716          | 0             , otherwise
 717           \
 718   \endverbatim
 719  * ***bd*** 2.2004
 720  */
 721 float CxImage::KernelLinear(const float t)
 722 {
 723 //  if (0<=t && t<=1) return 1-t;
 724 //  if (-1<=t && t<0) return 1+t;
 725 //  return 0;
 726
 727         //<Vladimir Kloucek>
 728         if (t < -1.0f)
 729                 return 0.0f;
 730         if (t < 0.0f)
 731                 return 1.0f+t;
 732         if (t < 1.0f)
 733                 return 1.0f-t;
 734         return 0.0f;
 735 }
 736
 737 ////////////////////////////////////////////////////////////////////////////////
 738 /**
 739  * Bicubic interpolation kernel (a=-1):
 740   \verbatim
 741           /
 742          | 1-2|t|**2+|t|**3          , if |t| < 1
 743   h(t) = | 4-8|t|+5|t|**2-|t|**3     , if 1<=|t|<2
 744          | 0                         , otherwise
 745           \
 746   \endverbatim
 747  * ***bd*** 2.2004
 748  */
 749 float CxImage::KernelCubic(const float t)
 750 {
 751   float abs_t = (float)fabs(t);
 752   float abs_t_sq = abs_t * abs_t;
 753   if (abs_t<1) return 1-2*abs_t_sq+abs_t_sq*abs_t;
 754   if (abs_t<2) return 4 - 8*abs_t +5*abs_t_sq - abs_t_sq*abs_t;
 755   return 0;
 756 }
 757
 758 ////////////////////////////////////////////////////////////////////////////////
 759 /**
 760  * Bicubic kernel (for a=-1 it is the same as BicubicKernel):
 761   \verbatim
 762           /
 763          | (a+2)|t|**3 - (a+3)|t|**2 + 1     , |t| <= 1
 764   h(t) = | a|t|**3 - 5a|t|**2 + 8a|t| - 4a   , 1 < |t| <= 2
 765          | 0                                 , otherwise
 766           \
 767   \endverbatim
 768  * Often used values for a are -1 and -1/2.
 769  */
 770 float CxImage::KernelGeneralizedCubic(const float t, const float a)
 771 {
 772   float abs_t = (float)fabs(t);
 773   float abs_t_sq = abs_t * abs_t;
 774   if (abs_t<1) return (a+2)*abs_t_sq*abs_t - (a+3)*abs_t_sq + 1;
 775   if (abs_t<2) return a*abs_t_sq*abs_t - 5*a*abs_t_sq + 8*a*abs_t - 4*a;
 776   return 0;
 777 }
 778
 779 ////////////////////////////////////////////////////////////////////////////////
 780 /**
 781  * Lanczos windowed sinc interpolation kernel with radius r.
 782   \verbatim
 783           /
 784   h(t) = | sinc(t)*sinc(t/r)       , if |t|<r
 785          | 0                       , otherwise
 786           \
 787   \endverbatim
 788  * ***bd*** 2.2004
 789  */
 790 float CxImage::KernelLanczosSinc(const float t, const float r)
 791 {
 792   if (fabs(t) > r) return 0;
 793   if (t==0) return 1;
 794   float pit=PI*t;
 795   float pitd=pit/r;
 796   return (float)((sin(pit)/pit) * (sin(pitd)/pitd));
 797 }
 798
 799 ////////////////////////////////////////////////////////////////////////////////
 800 float CxImage::KernelBox(const float x)
 801 {
 802         if (x < -0.5f)
 803                 return 0.0f;
 804         if (x < 0.5f)
 805                 return 1.0f;
 806         return 0.0f;
 807 }
 808 ////////////////////////////////////////////////////////////////////////////////
 809 float CxImage::KernelHermite(const float x)
 810 {
 811         if (x < -1.0f)
 812                 return 0.0f;
 813         if (x < 0.0f)
 814                 return (-2.0f*x-3.0f)*x*x+1.0f;
 815         if (x < 1.0f)
 816                 return (2.0f*x-3.0f)*x*x+1.0f;
 817         return 0.0f;
 818 //      if (fabs(x)>1) return 0.0f;
 819 //      return(0.5f+0.5f*(float)cos(PI*x));
 820 }
 821 ////////////////////////////////////////////////////////////////////////////////
 822 float CxImage::KernelHanning(const float x)
 823 {
 824         if (fabs(x)>1) return 0.0f;
 825         return (0.5f+0.5f*(float)cos(PI*x))*((float)sin(PI*x)/(PI*x));
 826 }
 827 ////////////////////////////////////////////////////////////////////////////////
 828 float CxImage::KernelHamming(const float x)
 829 {
 830         if (x < -1.0f)
 831                 return 0.0f;
 832         if (x < 0.0f)
 833                 return 0.92f*(-2.0f*x-3.0f)*x*x+1.0f;
 834         if (x < 1.0f)
 835                 return 0.92f*(2.0f*x-3.0f)*x*x+1.0f;
 836         return 0.0f;
 837 //      if (fabs(x)>1) return 0.0f;
 838 //      return(0.54f+0.46f*(float)cos(PI*x));
 839 }
 840 ////////////////////////////////////////////////////////////////////////////////
 841 float CxImage::KernelSinc(const float x)
 842 {
 843         if (x == 0.0)
 844                 return(1.0);
 845         return((float)sin(PI*x)/(PI*x));
 846 }
 847 ////////////////////////////////////////////////////////////////////////////////
 848 float CxImage::KernelBlackman(const float x)
 849 {
 850         //if (fabs(x)>1) return 0.0f;
 851         return (0.42f+0.5f*(float)cos(PI*x)+0.08f*(float)cos(2.0f*PI*x));
 852 }
 853 ////////////////////////////////////////////////////////////////////////////////
 854 float CxImage::KernelBessel_J1(const float x)
 855 {
 856         double p, q;
 857
 858         long i;
 859
 860         static const double
 861         Pone[] =
 862         {
 863                 0.581199354001606143928050809e+21,
 864                 -0.6672106568924916298020941484e+20,
 865                 0.2316433580634002297931815435e+19,
 866                 -0.3588817569910106050743641413e+17,
 867                 0.2908795263834775409737601689e+15,
 868                 -0.1322983480332126453125473247e+13,
 869                 0.3413234182301700539091292655e+10,
 870                 -0.4695753530642995859767162166e+7,
 871                 0.270112271089232341485679099e+4
 872         },
 873         Qone[] =
 874         {
 875                 0.11623987080032122878585294e+22,
 876                 0.1185770712190320999837113348e+20,
 877                 0.6092061398917521746105196863e+17,
 878                 0.2081661221307607351240184229e+15,
 879                 0.5243710262167649715406728642e+12,
 880                 0.1013863514358673989967045588e+10,
 881                 0.1501793594998585505921097578e+7,
 882                 0.1606931573481487801970916749e+4,
 883                 0.1e+1
 884         };
 885
 886         p = Pone[8];
 887         q = Qone[8];
 888         for (i=7; i >= 0; i--)
 889         {
 890                 p = p*x*x+Pone[i];
 891                 q = q*x*x+Qone[i];
 892         }
 893         return (float)(p/q);
 894 }
 895 ////////////////////////////////////////////////////////////////////////////////
 896 float CxImage::KernelBessel_P1(const float x)
 897 {
 898         double p, q;
 899
 900         long i;
 901
 902         static const double
 903         Pone[] =
 904         {
 905                 0.352246649133679798341724373e+5,
 906                 0.62758845247161281269005675e+5,
 907                 0.313539631109159574238669888e+5,
 908                 0.49854832060594338434500455e+4,
 909                 0.2111529182853962382105718e+3,
 910                 0.12571716929145341558495e+1
 911         },
 912         Qone[] =
 913         {
 914                 0.352246649133679798068390431e+5,
 915                 0.626943469593560511888833731e+5,
 916                 0.312404063819041039923015703e+5,
 917                 0.4930396490181088979386097e+4,
 918                 0.2030775189134759322293574e+3,
 919                 0.1e+1
 920         };
 921
 922         p = Pone[5];
 923         q = Qone[5];
 924         for (i=4; i >= 0; i--)
 925         {
 926                 p = p*(8.0/x)*(8.0/x)+Pone[i];
 927                 q = q*(8.0/x)*(8.0/x)+Qone[i];
 928         }
 929         return (float)(p/q);
 930 }
 931 ////////////////////////////////////////////////////////////////////////////////
 932 float CxImage::KernelBessel_Q1(const float x)
 933 {
 934         double p, q;
 935
 936         long i;
 937
 938         static const double
 939         Pone[] =
 940         {
 941                 0.3511751914303552822533318e+3,
 942                 0.7210391804904475039280863e+3,
 943                 0.4259873011654442389886993e+3,
 944                 0.831898957673850827325226e+2,
 945                 0.45681716295512267064405e+1,
 946                 0.3532840052740123642735e-1
 947         },
 948         Qone[] =
 949         {
 950                 0.74917374171809127714519505e+4,
 951                 0.154141773392650970499848051e+5,
 952                 0.91522317015169922705904727e+4,
 953                 0.18111867005523513506724158e+4,
 954                 0.1038187585462133728776636e+3,
 955                 0.1e+1
 956         };
 957
 958         p = Pone[5];
 959         q = Qone[5];
 960         for (i=4; i >= 0; i--)
 961         {
 962                 p = p*(8.0/x)*(8.0/x)+Pone[i];
 963                 q = q*(8.0/x)*(8.0/x)+Qone[i];
 964         }
 965         return (float)(p/q);
 966 }
 967 ////////////////////////////////////////////////////////////////////////////////
 968 float CxImage::KernelBessel_Order1(float x)
 969 {
 970         float p, q;
 971
 972         if (x == 0.0)
 973                 return (0.0f);
 974         p = x;
 975         if (x < 0.0)
 976                 x=(-x);
 977         if (x < 8.0)
 978                 return(p*KernelBessel_J1(x));
 979         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)*
 980                 (-1.0f/sqrt(2.0f)*(sin(x)+cos(x))));
 981         if (p < 0.0f)
 982                 q = (-q);
 983         return (q);
 984 }
 985 ////////////////////////////////////////////////////////////////////////////////
 986 float CxImage::KernelBessel(const float x)
 987 {
 988         if (x == 0.0f)
 989                 return(PI/4.0f);
 990         return(KernelBessel_Order1(PI*x)/(2.0f*x));
 991 }
 992 ////////////////////////////////////////////////////////////////////////////////
 993 float CxImage::KernelGaussian(const float x)
 994 {
 995         return (float)(exp(-2.0f*x*x)*0.79788456080287f/*sqrt(2.0f/PI)*/);
 996 }
 997 ////////////////////////////////////////////////////////////////////////////////
 998 float CxImage::KernelQuadratic(const float x)
 999 {
1000         if (x < -1.5f)
1001                 return(0.0f);
1002         if (x < -0.5f)
1003                 return(0.5f*(x+1.5f)*(x+1.5f));
1004         if (x < 0.5f)
1005                 return(0.75f-x*x);
1006         if (x < 1.5f)
1007                 return(0.5f*(x-1.5f)*(x-1.5f));
1008         return(0.0f);
1009 }
1010 ////////////////////////////////////////////////////////////////////////////////
1011 float CxImage::KernelMitchell(const float x)
1012 {
1013 #define KM_B (1.0f/3.0f)
1014 #define KM_C (1.0f/3.0f)
1015 #define KM_P0 ((  6.0f - 2.0f * KM_B ) / 6.0f)
1016 #define KM_P2 ((-18.0f + 12.0f * KM_B + 6.0f * KM_C) / 6.0f)
1017 #define KM_P3 (( 12.0f - 9.0f  * KM_B - 6.0f * KM_C) / 6.0f)
1018 #define KM_Q0 ((  8.0f * KM_B + 24.0f * KM_C) / 6.0f)
1019 #define KM_Q1 ((-12.0f * KM_B - 48.0f * KM_C) / 6.0f)
1020 #define KM_Q2 ((  6.0f * KM_B + 30.0f * KM_C) / 6.0f)
1021 #define KM_Q3 (( -1.0f * KM_B -  6.0f * KM_C) / 6.0f)
1022
1023         if (x < -2.0)
1024                 return(0.0f);
1025         if (x < -1.0)
1026                 return(KM_Q0-x*(KM_Q1-x*(KM_Q2-x*KM_Q3)));
1027         if (x < 0.0f)
1028                 return(KM_P0+x*x*(KM_P2-x*KM_P3));
1029         if (x < 1.0f)
1030                 return(KM_P0+x*x*(KM_P2+x*KM_P3));
1031         if (x < 2.0f)
1032                 return(KM_Q0+x*(KM_Q1+x*(KM_Q2+x*KM_Q3)));
1033         return(0.0f);
1034 }
1035 ////////////////////////////////////////////////////////////////////////////////
1036 float CxImage::KernelCatrom(const float x)
1037 {
1038         if (x < -2.0)
1039                 return(0.0f);
1040         if (x < -1.0)
1041                 return(0.5f*(4.0f+x*(8.0f+x*(5.0f+x))));
1042         if (x < 0.0)
1043                 return(0.5f*(2.0f+x*x*(-5.0f-3.0f*x)));
1044         if (x < 1.0)
1045                 return(0.5f*(2.0f+x*x*(-5.0f+3.0f*x)));
1046         if (x < 2.0)
1047                 return(0.5f*(4.0f+x*(-8.0f+x*(5.0f-x))));
1048         return(0.0f);
1049 }
1050 ////////////////////////////////////////////////////////////////////////////////
1051 float CxImage::KernelPower(const float x, const float a)
1052 {
1053         if (fabs(x)>1) return 0.0f;
1054         return (1.0f - (float)fabs(pow(x,a)));
1055 }
1056 ////////////////////////////////////////////////////////////////////////////////
1057
1058 #endif