From 0900cc583cc5e51b7c6b378de93cca29e11cb0e9 Mon Sep 17 00:00:00 2001 From: Chris Xiong Date: Fri, 28 Oct 2016 16:35:21 +0800 Subject: Add the GLFW port, still using OpenGL 2.1 though. And a more detailed readme. --- smelt/glfw/CxImage/ximaint.cpp | 1046 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1046 insertions(+) create mode 100644 smelt/glfw/CxImage/ximaint.cpp (limited to 'smelt/glfw/CxImage/ximaint.cpp') diff --git a/smelt/glfw/CxImage/ximaint.cpp b/smelt/glfw/CxImage/ximaint.cpp new file mode 100644 index 0000000..5f93038 --- /dev/null +++ b/smelt/glfw/CxImage/ximaint.cpp @@ -0,0 +1,1046 @@ +// xImaInt.cpp : interpolation functions +/* 02/2004 - Branko Brevensek + * CxImage version 7.0.0 31/Dec/2010 - Davide Pizzolato - www.xdp.it + */ + +#include "ximage.h" +#include "ximath.h" + +#if CXIMAGE_SUPPORT_INTERPOLATION + +//////////////////////////////////////////////////////////////////////////////// +/** + * Recalculates coordinates according to specified overflow method. + * If pixel (x,y) lies within image, nothing changes. + * + * \param x, y - coordinates of pixel + * \param ofMethod - overflow method + * + * \return x, y - new coordinates (pixel (x,y) now lies inside image) + * + * \author ***bd*** 2.2004 + */ +void CxImage::OverflowCoordinates(int32_t &x, int32_t &y, OverflowMethod const ofMethod) +{ + if (IsInside(x,y)) return; //if pixel is within bounds, no change + switch (ofMethod) { + case OM_REPEAT: + //clip coordinates + x=max(x,0); x=min(x, head.biWidth-1); + y=max(y,0); y=min(y, head.biHeight-1); + break; + case OM_WRAP: + //wrap coordinates + x = x % head.biWidth; + y = y % head.biHeight; + if (x<0) x = head.biWidth + x; + if (y<0) y = head.biHeight + y; + break; + case OM_MIRROR: + //mirror pixels near border + if (x<0) x=((-x) % head.biWidth); + else if (x>=head.biWidth) x=head.biWidth-(x % head.biWidth + 1); + if (y<0) y=((-y) % head.biHeight); + else if (y>=head.biHeight) y=head.biHeight-(y % head.biHeight + 1); + break; + default: + return; + }//switch +} + +//////////////////////////////////////////////////////////////////////////////// +/** + * See OverflowCoordinates for integer version + * \author ***bd*** 2.2004 + */ +void CxImage::OverflowCoordinates(float &x, float &y, OverflowMethod const ofMethod) +{ + if (x>=0 && x=0 && y=head.biWidth) x=head.biWidth-((float)fmod(x, (float) head.biWidth) + 1); + if (y<0) y=(float)fmod(-y, (float) head.biHeight); + else if (y>=head.biHeight) y=head.biHeight-((float)fmod(y, (float) head.biHeight) + 1); + break; + default: + return; + }//switch +} + +//////////////////////////////////////////////////////////////////////////////// +/** + * Method return pixel color. Different methods are implemented for out of bounds pixels. + * If an image has alpha channel, alpha value is returned in .RGBReserved. + * + * \param x,y : pixel coordinates + * \param ofMethod : out-of-bounds method: + * - OF_WRAP - wrap over to pixels on other side of the image + * - OF_REPEAT - repeat last pixel on the edge + * - OF_COLOR - return input value of color + * - OF_BACKGROUND - return background color (if not set, return input color) + * - OF_TRANSPARENT - return transparent pixel + * + * \param rplColor : input color (returned for out-of-bound coordinates in OF_COLOR mode and if other mode is not applicable) + * + * \return color : color of pixel + * \author ***bd*** 2.2004 + */ +RGBQUAD CxImage::GetPixelColorWithOverflow(int32_t x, int32_t y, OverflowMethod const ofMethod, RGBQUAD* const rplColor) +{ + RGBQUAD color; //color to return + if ((!IsInside(x,y)) || pDib==NULL) { //is pixel within bouns?: + //pixel is out of bounds or no DIB + if (rplColor!=NULL) + color=*rplColor; + else { + color.rgbRed=color.rgbGreen=color.rgbBlue=255; color.rgbReserved=0; //default replacement colour: white transparent + }//if + if (pDib==NULL) return color; + //pixel is out of bounds: + switch (ofMethod) { + case OM_TRANSPARENT: +#if CXIMAGE_SUPPORT_ALPHA + if (AlphaIsValid()) { + //alpha transparency is supported and image has alpha layer + color.rgbReserved=0; + } else { +#endif //CXIMAGE_SUPPORT_ALPHA + //no alpha transparency + if (GetTransIndex()>=0) { + color=GetTransColor(); //single color transparency enabled (return transparent color) + }//if +#if CXIMAGE_SUPPORT_ALPHA + }//if +#endif //CXIMAGE_SUPPORT_ALPHA + return color; + case OM_BACKGROUND: + //return background color (if it exists, otherwise input value) + if (info.nBkgndIndex >= 0) { + if (head.biBitCount<24) color = GetPaletteColor((uint8_t)info.nBkgndIndex); + else color = info.nBkgndColor; + }//if + return color; + case OM_REPEAT: + case OM_WRAP: + case OM_MIRROR: + OverflowCoordinates(x,y,ofMethod); + break; + default: + //simply return replacement color (OM_COLOR and others) + return color; + }//switch + }//if + //just return specified pixel (it's within bounds) + return BlindGetPixelColor(x,y); +} + +//////////////////////////////////////////////////////////////////////////////// +/** + * This method reconstructs image according to chosen interpolation method and then returns pixel (x,y). + * (x,y) can lie between actual image pixels. If (x,y) lies outside of image, method returns value + * according to overflow method. + * This method is very useful for geometrical image transformations, where destination pixel + * can often assume color value lying between source pixels. + * + * \param (x,y) - coordinates of pixel to return + * GPCI method recreates "analogue" image back from digital data, so x and y + * are float values and color value of point (1.1,1) will generally not be same + * as (1,1). Center of first pixel is at (0,0) and center of pixel right to it is (1,0). + * (0.5,0) is half way between these two pixels. + * \param inMethod - interpolation (reconstruction) method (kernel) to use: + * - IM_NEAREST_NEIGHBOUR - returns colour of nearest lying pixel (causes stairy look of + * processed images) + * - IM_BILINEAR - interpolates colour from four neighbouring pixels (softens image a bit) + * - IM_BICUBIC - interpolates from 16 neighbouring pixels (can produce "halo" artifacts) + * - IM_BICUBIC2 - interpolates from 16 neighbouring pixels (perhaps a bit less halo artifacts + than IM_BICUBIC) + * - IM_BSPLINE - interpolates from 16 neighbouring pixels (softens image, washes colours) + * (As far as I know, image should be prefiltered for this method to give + * good results... some other time :) ) + * This method uses bicubic interpolation kernel from CXImage 5.99a and older + * versions. + * - IM_LANCZOS - interpolates from 12*12 pixels (slow, ringing artifacts) + * + * \param ofMethod - overflow method (see comments at GetPixelColorWithOverflow) + * \param rplColor - pointer to color used for out of borders pixels in OM_COLOR mode + * (and other modes if colour can't calculated in a specified way) + * + * \return interpolated color value (including interpolated alpha value, if image has alpha layer) + * + * \author ***bd*** 2.2004 + */ +RGBQUAD CxImage::GetPixelColorInterpolated( + float x,float y, + InterpolationMethod const inMethod, + OverflowMethod const ofMethod, + RGBQUAD* const rplColor) +{ + //calculate nearest pixel + int32_t xi=(int32_t)(x); if (x<0) xi--; //these replace (incredibly slow) floor (Visual c++ 2003, AMD Athlon) + int32_t yi=(int32_t)(y); if (y<0) yi--; + RGBQUAD color; //calculated colour + + switch (inMethod) { + case IM_NEAREST_NEIGHBOUR: + return GetPixelColorWithOverflow((int32_t)(x+0.5f), (int32_t)(y+0.5f), ofMethod, rplColor); + default: { + //IM_BILINEAR: bilinear interpolation + if (xi<-1 || xi>=head.biWidth || yi<-1 || yi>=head.biHeight) { //all 4 points are outside bounds?: + switch (ofMethod) { + case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND: + //we don't need to interpolate anything with all points outside in this case + return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor); + default: + //recalculate coordinates and use faster method later on + OverflowCoordinates(x,y,ofMethod); + xi=(int32_t)(x); if (x<0) xi--; //x and/or y have changed ... recalculate xi and yi + yi=(int32_t)(y); if (y<0) yi--; + }//switch + }//if + //get four neighbouring pixels + if ((xi+1)=0 && (yi+1)=0 && head.biClrUsed==0) { + //all pixels are inside RGB24 image... optimize reading (and use fixed point arithmetic) + uint16_t wt1=(uint16_t)((x-xi)*256.0f), wt2=(uint16_t)((y-yi)*256.0f); + uint16_t wd=wt1*wt2>>8; + uint16_t wb=wt1-wd; + uint16_t wc=wt2-wd; + uint16_t wa=256-wt1-wc; + uint16_t wrr,wgg,wbb; + uint8_t *pxptr=(uint8_t*)info.pImage+yi*info.dwEffWidth+xi*3; + wbb=wa*(*pxptr++); wgg=wa*(*pxptr++); wrr=wa*(*pxptr++); + wbb+=wb*(*pxptr++); wgg+=wb*(*pxptr++); wrr+=wb*(*pxptr); + pxptr+=(info.dwEffWidth-5); //move to next row + wbb+=wc*(*pxptr++); wgg+=wc*(*pxptr++); wrr+=wc*(*pxptr++); + wbb+=wd*(*pxptr++); wgg+=wd*(*pxptr++); wrr+=wd*(*pxptr); + color.rgbRed=(uint8_t) (wrr>>8); color.rgbGreen=(uint8_t) (wgg>>8); color.rgbBlue=(uint8_t) (wbb>>8); +#if CXIMAGE_SUPPORT_ALPHA + if (pAlpha) { + uint16_t waa; + //image has alpha layer... we have to do the same for alpha data + pxptr=AlphaGetPointer(xi,yi); //pointer to first byte + waa=wa*(*pxptr++); waa+=wb*(*pxptr); //first two pixels + pxptr+=(head.biWidth-1); //move to next row + waa+=wc*(*pxptr++); waa+=wd*(*pxptr); //and second row pixels + color.rgbReserved=(uint8_t) (waa>>8); + } else +#endif + { //Alpha not supported or no alpha at all + color.rgbReserved = 0; + } + return color; + } else { + //default (slower) way to get pixels (not RGB24 or some pixels out of borders) + float t1=x-xi, t2=y-yi; + float d=t1*t2; + float b=t1-d; + float c=t2-d; + float a=1-t1-c; + RGBQUAD rgb11,rgb21,rgb12,rgb22; + rgb11=GetPixelColorWithOverflow(xi, yi, ofMethod, rplColor); + rgb21=GetPixelColorWithOverflow(xi+1, yi, ofMethod, rplColor); + rgb12=GetPixelColorWithOverflow(xi, yi+1, ofMethod, rplColor); + rgb22=GetPixelColorWithOverflow(xi+1, yi+1, ofMethod, rplColor); + //calculate linear interpolation + color.rgbRed=(uint8_t) (a*rgb11.rgbRed+b*rgb21.rgbRed+c*rgb12.rgbRed+d*rgb22.rgbRed); + color.rgbGreen=(uint8_t) (a*rgb11.rgbGreen+b*rgb21.rgbGreen+c*rgb12.rgbGreen+d*rgb22.rgbGreen); + color.rgbBlue=(uint8_t) (a*rgb11.rgbBlue+b*rgb21.rgbBlue+c*rgb12.rgbBlue+d*rgb22.rgbBlue); +#if CXIMAGE_SUPPORT_ALPHA + color.rgbReserved=(uint8_t) (a*rgb11.rgbReserved+b*rgb21.rgbReserved+c*rgb12.rgbReserved+d*rgb22.rgbReserved); +#else + color.rgbReserved = 0; +#endif + return color; + }//if + }//default + case IM_BICUBIC: + case IM_BICUBIC2: + case IM_BSPLINE: + case IM_BOX: + case IM_HERMITE: + case IM_HAMMING: + case IM_SINC: + case IM_BLACKMAN: + case IM_BESSEL: + case IM_GAUSSIAN: + case IM_QUADRATIC: + case IM_MITCHELL: + case IM_CATROM: + case IM_HANNING: + case IM_POWER: + //bicubic interpolation(s) + if (((xi+2)<0) || ((xi-1)>=head.biWidth) || ((yi+2)<0) || ((yi-1)>=head.biHeight)) { //all points are outside bounds?: + switch (ofMethod) { + case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND: + //we don't need to interpolate anything with all points outside in this case + return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor); + break; + default: + //recalculate coordinates and use faster method later on + OverflowCoordinates(x,y,ofMethod); + xi=(int32_t)(x); if (x<0) xi--; //x and/or y have changed ... recalculate xi and yi + yi=(int32_t)(y); if (y<0) yi--; + }//switch + }//if + + //some variables needed from here on + int32_t xii,yii; //x any y integer indexes for loops + float kernel, kernelyc; //kernel cache + float kernelx[12], kernely[4]; //precalculated kernel values + float rr,gg,bb,aa; //accumulated color values + //calculate multiplication factors for all pixels + int32_t i; + switch (inMethod) { + case IM_BICUBIC: + for (i=0; i<4; i++) { + kernelx[i]=KernelCubic((float)(xi+i-1-x)); + kernely[i]=KernelCubic((float)(yi+i-1-y)); + }//for i + break; + case IM_BICUBIC2: + for (i=0; i<4; i++) { + kernelx[i]=KernelGeneralizedCubic((float)(xi+i-1-x), -0.5); + kernely[i]=KernelGeneralizedCubic((float)(yi+i-1-y), -0.5); + }//for i + break; + case IM_BSPLINE: + for (i=0; i<4; i++) { + kernelx[i]=KernelBSpline((float)(xi+i-1-x)); + kernely[i]=KernelBSpline((float)(yi+i-1-y)); + }//for i + break; + case IM_BOX: + for (i=0; i<4; i++) { + kernelx[i]=KernelBox((float)(xi+i-1-x)); + kernely[i]=KernelBox((float)(yi+i-1-y)); + }//for i + break; + case IM_HERMITE: + for (i=0; i<4; i++) { + kernelx[i]=KernelHermite((float)(xi+i-1-x)); + kernely[i]=KernelHermite((float)(yi+i-1-y)); + }//for i + break; + case IM_HAMMING: + for (i=0; i<4; i++) { + kernelx[i]=KernelHamming((float)(xi+i-1-x)); + kernely[i]=KernelHamming((float)(yi+i-1-y)); + }//for i + break; + case IM_SINC: + for (i=0; i<4; i++) { + kernelx[i]=KernelSinc((float)(xi+i-1-x)); + kernely[i]=KernelSinc((float)(yi+i-1-y)); + }//for i + break; + case IM_BLACKMAN: + for (i=0; i<4; i++) { + kernelx[i]=KernelBlackman((float)(xi+i-1-x)); + kernely[i]=KernelBlackman((float)(yi+i-1-y)); + }//for i + break; + case IM_BESSEL: + for (i=0; i<4; i++) { + kernelx[i]=KernelBessel((float)(xi+i-1-x)); + kernely[i]=KernelBessel((float)(yi+i-1-y)); + }//for i + break; + case IM_GAUSSIAN: + for (i=0; i<4; i++) { + kernelx[i]=KernelGaussian((float)(xi+i-1-x)); + kernely[i]=KernelGaussian((float)(yi+i-1-y)); + }//for i + break; + case IM_QUADRATIC: + for (i=0; i<4; i++) { + kernelx[i]=KernelQuadratic((float)(xi+i-1-x)); + kernely[i]=KernelQuadratic((float)(yi+i-1-y)); + }//for i + break; + case IM_MITCHELL: + for (i=0; i<4; i++) { + kernelx[i]=KernelMitchell((float)(xi+i-1-x)); + kernely[i]=KernelMitchell((float)(yi+i-1-y)); + }//for i + break; + case IM_CATROM: + for (i=0; i<4; i++) { + kernelx[i]=KernelCatrom((float)(xi+i-1-x)); + kernely[i]=KernelCatrom((float)(yi+i-1-y)); + }//for i + break; + case IM_HANNING: + for (i=0; i<4; i++) { + kernelx[i]=KernelHanning((float)(xi+i-1-x)); + kernely[i]=KernelHanning((float)(yi+i-1-y)); + }//for i + break; + case IM_POWER: + for (i=0; i<4; i++) { + kernelx[i]=KernelPower((float)(xi+i-1-x)); + kernely[i]=KernelPower((float)(yi+i-1-y)); + }//for i + break; + default:break; + }//switch + rr=gg=bb=aa=0; + if (((xi+2)=1 && ((yi+2)=1) && !IsIndexed()) { + //optimized interpolation (faster pixel reads) for RGB24 images with all pixels inside bounds + for (yii=yi-1; yii255) rr=255; if (rr<0) rr=0; color.rgbRed=(uint8_t) rr; + if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(uint8_t) gg; + if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(uint8_t) bb; +#if CXIMAGE_SUPPORT_ALPHA + if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(uint8_t) aa; +#else + color.rgbReserved = 0; +#endif + return color; + case IM_LANCZOS: + //lanczos window (16*16) sinc interpolation + if (((xi+6)<0) || ((xi-5)>=head.biWidth) || ((yi+6)<0) || ((yi-5)>=head.biHeight)) { + //all points are outside bounds + switch (ofMethod) { + case OM_COLOR: case OM_TRANSPARENT: case OM_BACKGROUND: + //we don't need to interpolate anything with all points outside in this case + return GetPixelColorWithOverflow(-999, -999, ofMethod, rplColor); + break; + default: + //recalculate coordinates and use faster method later on + OverflowCoordinates(x,y,ofMethod); + xi=(int32_t)(x); if (x<0) xi--; //x and/or y have changed ... recalculate xi and yi + yi=(int32_t)(y); if (y<0) yi--; + }//switch + }//if + + for (xii=xi-5; xii=0) && ((yi+6)=0) && !IsIndexed()) { + //optimized interpolation (faster pixel reads) for RGB24 images with all pixels inside bounds + for (yii=yi-5; yii255) rr=255; if (rr<0) rr=0; color.rgbRed=(uint8_t) rr; + if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(uint8_t) gg; + if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(uint8_t) bb; +#if CXIMAGE_SUPPORT_ALPHA + if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(uint8_t) aa; +#else + color.rgbReserved = 0; +#endif + return color; + }//switch +} +//////////////////////////////////////////////////////////////////////////////// +/** + * Helper function for GetAreaColorInterpolated. + * Adds 'surf' portion of image pixel with color 'color' to (rr,gg,bb,aa). + */ +void CxImage::AddAveragingCont(RGBQUAD const &color, float const surf, float &rr, float &gg, float &bb, float &aa) +{ + rr+=color.rgbRed*surf; + gg+=color.rgbGreen*surf; + bb+=color.rgbBlue*surf; +#if CXIMAGE_SUPPORT_ALPHA + aa+=color.rgbReserved*surf; +#endif +} +//////////////////////////////////////////////////////////////////////////////// +/** + * This method is similar to GetPixelColorInterpolated, but this method also properly handles + * subsampling. + * If you need to sample original image with interval of more than 1 pixel (as when shrinking an image), + * you should use this method instead of GetPixelColorInterpolated or aliasing will occur. + * When area width and height are both less than pixel, this method gets pixel color by interpolating + * color of frame center with selected (inMethod) interpolation by calling GetPixelColorInterpolated. + * If width and height are more than 1, method calculates color by averaging color of pixels within area. + * Interpolation method is not used in this case. Pixel color is interpolated by averaging instead. + * If only one of both is more than 1, method uses combination of interpolation and averaging. + * Chosen interpolation method is used, but since it is averaged later on, there is little difference + * between IM_BILINEAR (perhaps best for this case) and better methods. IM_NEAREST_NEIGHBOUR again + * leads to aliasing artifacts. + * This method is a bit slower than GetPixelColorInterpolated and when aliasing is not a problem, you should + * simply use the later. + * + * \param xc, yc - center of (rectangular) area + * \param w, h - width and height of area + * \param inMethod - interpolation method that is used, when interpolation is used (see above) + * \param ofMethod - overflow method used when retrieving individual pixel colors + * \param rplColor - replacement colour to use, in OM_COLOR + * + * \author ***bd*** 2.2004 + */ +RGBQUAD CxImage::GetAreaColorInterpolated( + float const xc, float const yc, float const w, float const h, + InterpolationMethod const inMethod, + OverflowMethod const ofMethod, + RGBQUAD* const rplColor) +{ + RGBQUAD color; //calculated colour + + if (h<=1 && w<=1) { + //both width and height are less than one... we will use interpolation of center point + return GetPixelColorInterpolated(xc, yc, inMethod, ofMethod, rplColor); + } else { + //area is wider and/or taller than one pixel: + CxRect2 area(xc-w/2.0f, yc-h/2.0f, xc+w/2.0f, yc+h/2.0f); //area + int32_t xi1=(int32_t)(area.botLeft.x+0.49999999f); //low x + int32_t yi1=(int32_t)(area.botLeft.y+0.49999999f); //low y + + + int32_t xi2=(int32_t)(area.topRight.x+0.5f); //top x + int32_t yi2=(int32_t)(area.topRight.y+0.5f); //top y (for loops) + + float rr,gg,bb,aa; //red, green, blue and alpha components + rr=gg=bb=aa=0; + int32_t x,y; //loop counters + float s=0; //surface of all pixels + float cps; //surface of current crosssection + if (h>1 && w>1) { + //width and height of area are greater than one pixel, so we can employ "ordinary" averaging + CxRect2 intBL, intTR; //bottom left and top right intersection + intBL=area.CrossSection(CxRect2(((float)xi1)-0.5f, ((float)yi1)-0.5f, ((float)xi1)+0.5f, ((float)yi1)+0.5f)); + intTR=area.CrossSection(CxRect2(((float)xi2)-0.5f, ((float)yi2)-0.5f, ((float)xi2)+0.5f, ((float)yi2)+0.5f)); + float wBL, wTR, hBL, hTR; + wBL=intBL.Width(); //width of bottom left pixel-area intersection + hBL=intBL.Height(); //height of bottom left... + wTR=intTR.Width(); //width of top right... + hTR=intTR.Height(); //height of top right... + + AddAveragingCont(GetPixelColorWithOverflow(xi1,yi1,ofMethod,rplColor), wBL*hBL, rr, gg, bb, aa); //bottom left pixel + AddAveragingCont(GetPixelColorWithOverflow(xi2,yi1,ofMethod,rplColor), wTR*hBL, rr, gg, bb, aa); //bottom right pixel + AddAveragingCont(GetPixelColorWithOverflow(xi1,yi2,ofMethod,rplColor), wBL*hTR, rr, gg, bb, aa); //top left pixel + AddAveragingCont(GetPixelColorWithOverflow(xi2,yi2,ofMethod,rplColor), wTR*hTR, rr, gg, bb, aa); //top right pixel + //bottom and top row + for (x=xi1+1; x255) rr=255; if (rr<0) rr=0; color.rgbRed=(uint8_t) rr; + if (gg>255) gg=255; if (gg<0) gg=0; color.rgbGreen=(uint8_t) gg; + if (bb>255) bb=255; if (bb<0) bb=0; color.rgbBlue=(uint8_t) bb; +#if CXIMAGE_SUPPORT_ALPHA + if (aa>255) aa=255; if (aa<0) aa=0; color.rgbReserved=(uint8_t) aa; +#else + color.rgbReserved = 0; +#endif + }//if + return color; +} + +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBSpline(const float x) +{ + if (x>2.0f) return 0.0f; + // thanks to Kristian Kratzenstein + float a, b, c, d; + float xm1 = x - 1.0f; // Was calculatet anyway cause the "if((x-1.0f) < 0)" + float xp1 = x + 1.0f; + float xp2 = x + 2.0f; + + if ((xp2) <= 0.0f) a = 0.0f; else a = xp2*xp2*xp2; // Only float, not float -> double -> float + if ((xp1) <= 0.0f) b = 0.0f; else b = xp1*xp1*xp1; + if (x <= 0) c = 0.0f; else c = x*x*x; + if ((xm1) <= 0.0f) d = 0.0f; else d = xm1*xm1*xm1; + + return (0.16666666666666666667f * (a - (4.0f * b) + (6.0f * c) - (4.0f * d))); + + /* equivalent + if (x < -2.0) + return(0.0f); + if (x < -1.0) + return((2.0f+x)*(2.0f+x)*(2.0f+x)*0.16666666666666666667f); + if (x < 0.0) + return((4.0f+x*x*(-6.0f-3.0f*x))*0.16666666666666666667f); + if (x < 1.0) + return((4.0f+x*x*(-6.0f+3.0f*x))*0.16666666666666666667f); + if (x < 2.0) + return((2.0f-x)*(2.0f-x)*(2.0f-x)*0.16666666666666666667f); + return(0.0f); + */ +} + +//////////////////////////////////////////////////////////////////////////////// +/** + * Bilinear interpolation kernel: + \verbatim + / + | 1-t , if 0 <= t <= 1 + h(t) = | t+1 , if -1 <= t < 0 + | 0 , otherwise + \ + \endverbatim + * ***bd*** 2.2004 + */ +float CxImage::KernelLinear(const float t) +{ +// if (0<=t && t<=1) return 1-t; +// if (-1<=t && t<0) return 1+t; +// return 0; + + // + if (t < -1.0f) + return 0.0f; + if (t < 0.0f) + return 1.0f+t; + if (t < 1.0f) + return 1.0f-t; + return 0.0f; +} + +//////////////////////////////////////////////////////////////////////////////// +/** + * Bicubic interpolation kernel (a=-1): + \verbatim + / + | 1-2|t|**2+|t|**3 , if |t| < 1 + h(t) = | 4-8|t|+5|t|**2-|t|**3 , if 1<=|t|<2 + | 0 , otherwise + \ + \endverbatim + * ***bd*** 2.2004 + */ +float CxImage::KernelCubic(const float t) +{ + float abs_t = (float)fabs(t); + float abs_t_sq = abs_t * abs_t; + if (abs_t<1) return 1-2*abs_t_sq+abs_t_sq*abs_t; + if (abs_t<2) return 4 - 8*abs_t +5*abs_t_sq - abs_t_sq*abs_t; + return 0; +} + +//////////////////////////////////////////////////////////////////////////////// +/** + * Bicubic kernel (for a=-1 it is the same as BicubicKernel): + \verbatim + / + | (a+2)|t|**3 - (a+3)|t|**2 + 1 , |t| <= 1 + h(t) = | a|t|**3 - 5a|t|**2 + 8a|t| - 4a , 1 < |t| <= 2 + | 0 , otherwise + \ + \endverbatim + * Often used values for a are -1 and -1/2. + */ +float CxImage::KernelGeneralizedCubic(const float t, const float a) +{ + float abs_t = (float)fabs(t); + float abs_t_sq = abs_t * abs_t; + if (abs_t<1) return (a+2)*abs_t_sq*abs_t - (a+3)*abs_t_sq + 1; + if (abs_t<2) return a*abs_t_sq*abs_t - 5*a*abs_t_sq + 8*a*abs_t - 4*a; + return 0; +} + +//////////////////////////////////////////////////////////////////////////////// +/** + * Lanczos windowed sinc interpolation kernel with radius r. + \verbatim + / + h(t) = | sinc(t)*sinc(t/r) , if |t| r) return 0; + if (t==0) return 1; + float pit=PI*t; + float pitd=pit/r; + return (float)((sin(pit)/pit) * (sin(pitd)/pitd)); +} + +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBox(const float x) +{ + if (x < -0.5f) + return 0.0f; + if (x < 0.5f) + return 1.0f; + return 0.0f; +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelHermite(const float x) +{ + if (x < -1.0f) + return 0.0f; + if (x < 0.0f) + return (-2.0f*x-3.0f)*x*x+1.0f; + if (x < 1.0f) + return (2.0f*x-3.0f)*x*x+1.0f; + return 0.0f; +// if (fabs(x)>1) return 0.0f; +// return(0.5f+0.5f*(float)cos(PI*x)); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelHanning(const float x) +{ + if (fabs(x)>1) return 0.0f; + return (0.5f+0.5f*(float)cos(PI*x))*((float)sin(PI*x)/(PI*x)); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelHamming(const float x) +{ + if (x < -1.0f) + return 0.0f; + if (x < 0.0f) + return 0.92f*(-2.0f*x-3.0f)*x*x+1.0f; + if (x < 1.0f) + return 0.92f*(2.0f*x-3.0f)*x*x+1.0f; + return 0.0f; +// if (fabs(x)>1) return 0.0f; +// return(0.54f+0.46f*(float)cos(PI*x)); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelSinc(const float x) +{ + if (x == 0.0) + return(1.0); + return((float)sin(PI*x)/(PI*x)); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBlackman(const float x) +{ + //if (fabs(x)>1) return 0.0f; + return (0.42f+0.5f*(float)cos(PI*x)+0.08f*(float)cos(2.0f*PI*x)); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBessel_J1(const float x) +{ + double p, q; + + register int32_t i; + + static const double + Pone[] = + { + 0.581199354001606143928050809e+21, + -0.6672106568924916298020941484e+20, + 0.2316433580634002297931815435e+19, + -0.3588817569910106050743641413e+17, + 0.2908795263834775409737601689e+15, + -0.1322983480332126453125473247e+13, + 0.3413234182301700539091292655e+10, + -0.4695753530642995859767162166e+7, + 0.270112271089232341485679099e+4 + }, + Qone[] = + { + 0.11623987080032122878585294e+22, + 0.1185770712190320999837113348e+20, + 0.6092061398917521746105196863e+17, + 0.2081661221307607351240184229e+15, + 0.5243710262167649715406728642e+12, + 0.1013863514358673989967045588e+10, + 0.1501793594998585505921097578e+7, + 0.1606931573481487801970916749e+4, + 0.1e+1 + }; + + p = Pone[8]; + q = Qone[8]; + for (i=7; i >= 0; i--) + { + p = p*x*x+Pone[i]; + q = q*x*x+Qone[i]; + } + return (float)(p/q); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBessel_P1(const float x) +{ + double p, q; + + register int32_t i; + + static const double + Pone[] = + { + 0.352246649133679798341724373e+5, + 0.62758845247161281269005675e+5, + 0.313539631109159574238669888e+5, + 0.49854832060594338434500455e+4, + 0.2111529182853962382105718e+3, + 0.12571716929145341558495e+1 + }, + Qone[] = + { + 0.352246649133679798068390431e+5, + 0.626943469593560511888833731e+5, + 0.312404063819041039923015703e+5, + 0.4930396490181088979386097e+4, + 0.2030775189134759322293574e+3, + 0.1e+1 + }; + + p = Pone[5]; + q = Qone[5]; + for (i=4; i >= 0; i--) + { + p = p*(8.0/x)*(8.0/x)+Pone[i]; + q = q*(8.0/x)*(8.0/x)+Qone[i]; + } + return (float)(p/q); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBessel_Q1(const float x) +{ + double p, q; + + register int32_t i; + + static const double + Pone[] = + { + 0.3511751914303552822533318e+3, + 0.7210391804904475039280863e+3, + 0.4259873011654442389886993e+3, + 0.831898957673850827325226e+2, + 0.45681716295512267064405e+1, + 0.3532840052740123642735e-1 + }, + Qone[] = + { + 0.74917374171809127714519505e+4, + 0.154141773392650970499848051e+5, + 0.91522317015169922705904727e+4, + 0.18111867005523513506724158e+4, + 0.1038187585462133728776636e+3, + 0.1e+1 + }; + + p = Pone[5]; + q = Qone[5]; + for (i=4; i >= 0; i--) + { + p = p*(8.0/x)*(8.0/x)+Pone[i]; + q = q*(8.0/x)*(8.0/x)+Qone[i]; + } + return (float)(p/q); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBessel_Order1(float x) +{ + float p, q; + + if (x == 0.0) + return (0.0f); + p = x; + if (x < 0.0) + x=(-x); + if (x < 8.0) + return(p*KernelBessel_J1(x)); + 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)* + (-1.0f/sqrt(2.0f)*(sin(x)+cos(x)))); + if (p < 0.0f) + q = (-q); + return (q); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelBessel(const float x) +{ + if (x == 0.0f) + return(PI/4.0f); + return(KernelBessel_Order1(PI*x)/(2.0f*x)); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelGaussian(const float x) +{ + return (float)(exp(-2.0f*x*x)*0.79788456080287f/*sqrt(2.0f/PI)*/); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelQuadratic(const float x) +{ + if (x < -1.5f) + return(0.0f); + if (x < -0.5f) + return(0.5f*(x+1.5f)*(x+1.5f)); + if (x < 0.5f) + return(0.75f-x*x); + if (x < 1.5f) + return(0.5f*(x-1.5f)*(x-1.5f)); + return(0.0f); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelMitchell(const float x) +{ +#define KM_B (1.0f/3.0f) +#define KM_C (1.0f/3.0f) +#define KM_P0 (( 6.0f - 2.0f * KM_B ) / 6.0f) +#define KM_P2 ((-18.0f + 12.0f * KM_B + 6.0f * KM_C) / 6.0f) +#define KM_P3 (( 12.0f - 9.0f * KM_B - 6.0f * KM_C) / 6.0f) +#define KM_Q0 (( 8.0f * KM_B + 24.0f * KM_C) / 6.0f) +#define KM_Q1 ((-12.0f * KM_B - 48.0f * KM_C) / 6.0f) +#define KM_Q2 (( 6.0f * KM_B + 30.0f * KM_C) / 6.0f) +#define KM_Q3 (( -1.0f * KM_B - 6.0f * KM_C) / 6.0f) + + if (x < -2.0) + return(0.0f); + if (x < -1.0) + return(KM_Q0-x*(KM_Q1-x*(KM_Q2-x*KM_Q3))); + if (x < 0.0f) + return(KM_P0+x*x*(KM_P2-x*KM_P3)); + if (x < 1.0f) + return(KM_P0+x*x*(KM_P2+x*KM_P3)); + if (x < 2.0f) + return(KM_Q0+x*(KM_Q1+x*(KM_Q2+x*KM_Q3))); + return(0.0f); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelCatrom(const float x) +{ + if (x < -2.0) + return(0.0f); + if (x < -1.0) + return(0.5f*(4.0f+x*(8.0f+x*(5.0f+x)))); + if (x < 0.0) + return(0.5f*(2.0f+x*x*(-5.0f-3.0f*x))); + if (x < 1.0) + return(0.5f*(2.0f+x*x*(-5.0f+3.0f*x))); + if (x < 2.0) + return(0.5f*(4.0f+x*(-8.0f+x*(5.0f-x)))); + return(0.0f); +} +//////////////////////////////////////////////////////////////////////////////// +float CxImage::KernelPower(const float x, const float a) +{ + if (fabs(x)>1) return 0.0f; + return (1.0f - (float)fabs(pow(x,a))); +} +//////////////////////////////////////////////////////////////////////////////// + +#endif -- cgit v1.2.3