utilities/CxImage/ximaint.cpp

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