api/qgscubicrasterresampler_8cpp_source.html

00001 /***************************************************************************
00002                          qgscubicrasterresampler.cpp
00003                          ----------------------------
00004     begin                : December 2011
00005     copyright            : (C) 2011 by Marco Hugentobler
00006     email                : marco at sourcepole dot ch
00007  ***************************************************************************/
00008
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017
00018 #include "qgscubicrasterresampler.h"
00019 #include <QImage>
00020 #include <cmath>
00021
00022 QgsCubicRasterResampler::QgsCubicRasterResampler()
00023 {
00024 }
00025
00026 QgsCubicRasterResampler::~QgsCubicRasterResampler()
00027 {
00028 }
00029
00030 QgsRasterResampler * QgsCubicRasterResampler::clone() const
00031 {
00032   return new QgsCubicRasterResampler();
00033 }
00034
00035 void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage )
00036 {
00037   int nCols = srcImage.width();
00038   int nRows = srcImage.height();
00039
00040   int pos = 0;
00041   QRgb px;
00042   int* redMatrix = new int[ nCols * nRows ];
00043   int* greenMatrix = new int[ nCols * nRows ];
00044   int* blueMatrix = new int[ nCols * nRows ];
00045   int* alphaMatrix = new int[ nCols * nRows ];
00046
00047   for ( int i = 0; i < nRows; ++i )
00048   {
00049     for ( int j = 0; j < nCols; ++j )
00050     {
00051       px = srcImage.pixel( j, i );
00052       redMatrix[pos] = qRed( px );
00053       greenMatrix[pos] = qGreen( px );
00054       blueMatrix[pos] = qBlue( px );
00055       alphaMatrix[pos] = qAlpha( px );
00056       ++pos;
00057     }
00058   }
00059
00060   //derivative x
00061   double* xDerivativeMatrixRed = new double[ nCols * nRows ];
00062   xDerivativeMatrix( nCols, nRows, xDerivativeMatrixRed, redMatrix );
00063   double* xDerivativeMatrixGreen = new double[ nCols * nRows ];
00064   xDerivativeMatrix( nCols, nRows, xDerivativeMatrixGreen, greenMatrix );
00065   double* xDerivativeMatrixBlue = new double[ nCols * nRows ];
00066   xDerivativeMatrix( nCols, nRows, xDerivativeMatrixBlue, blueMatrix );
00067   double* xDerivativeMatrixAlpha = new double[ nCols * nRows ];
00068   xDerivativeMatrix( nCols, nRows, xDerivativeMatrixAlpha, alphaMatrix );
00069
00070   //derivative y
00071   double* yDerivativeMatrixRed = new double[ nCols * nRows ];
00072   yDerivativeMatrix( nCols, nRows, yDerivativeMatrixRed, redMatrix );
00073   double* yDerivativeMatrixGreen = new double[ nCols * nRows ];
00074   yDerivativeMatrix( nCols, nRows, yDerivativeMatrixGreen, greenMatrix );
00075   double* yDerivativeMatrixBlue = new double[ nCols * nRows ];
00076   yDerivativeMatrix( nCols, nRows, yDerivativeMatrixBlue, blueMatrix );
00077   double* yDerivativeMatrixAlpha = new double[ nCols * nRows ];
00078   yDerivativeMatrix( nCols, nRows, yDerivativeMatrixAlpha, alphaMatrix );
00079
00080   //compute output
00081   double nSrcPerDstX = ( double ) srcImage.width() / ( double ) dstImage.width();
00082   double nSrcPerDstY = ( double ) srcImage.height() / ( double ) dstImage.height();
00083
00084   double currentSrcRow = nSrcPerDstY / 2.0 - 0.5;
00085   double currentSrcCol;
00086   int currentSrcColInt;
00087   int currentSrcRowInt;
00088   int lastSrcColInt = -1;
00089   int lastSrcRowInt = -1;
00090
00091   int r, g, b, a;
00092   //bernstein polynomials
00093   double bp0u, bp1u, bp2u, bp3u, bp0v, bp1v, bp2v, bp3v;
00094   double u, v;
00095
00096   for ( int i = 0; i < dstImage.height(); ++i )
00097   {
00098     currentSrcRowInt = floor( currentSrcRow );
00099     v = currentSrcRow - currentSrcRowInt;
00100
00101     currentSrcCol = nSrcPerDstX / 2.0 - 0.5;
00102     for ( int j = 0; j < dstImage.width(); ++j )
00103     {
00104       currentSrcColInt = floor( currentSrcCol );
00105       u = currentSrcCol - currentSrcColInt;
00106
00107       //handle eight edge-cases
00108       if ( currentSrcRowInt < 0 || currentSrcRowInt >= ( srcImage.height() - 1 ) || currentSrcColInt < 0 || currentSrcColInt >= ( srcImage.width() - 1 ) )
00109       {
00110         QRgb px1, px2;
00111         //pixels at the border of the source image needs to be handled in a special way
00112         if ( currentSrcRowInt < 0 && currentSrcColInt < 0 )
00113         {
00114           dstImage.setPixel( j, i, srcImage.pixel( 0, 0 ) );
00115         }
00116         else if ( currentSrcRowInt < 0 && currentSrcColInt >= ( srcImage.width() - 1 ) )
00117         {
00118           dstImage.setPixel( j, i, srcImage.pixel( srcImage.width() - 1, 0 ) );
00119         }
00120         else if ( currentSrcRowInt >= ( srcImage.height() - 1 ) && currentSrcColInt >= ( srcImage.width() - 1 ) )
00121         {
00122           dstImage.setPixel( j, i, srcImage.pixel( srcImage.width() - 1, srcImage.height() - 1 ) );
00123         }
00124         else if ( currentSrcRowInt >= ( srcImage.height() - 1 ) && currentSrcColInt < 0 )
00125         {
00126           dstImage.setPixel( j, i, srcImage.pixel( 0, srcImage.height() - 1 ) );
00127         }
00128         else if ( currentSrcRowInt < 0 )
00129         {
00130           px1 = srcImage.pixel( currentSrcColInt, 0 );
00131           px2 = srcImage.pixel( currentSrcColInt + 1, 0 );
00132           dstImage.setPixel( j, i, curveInterpolation( px1, px2, u, xDerivativeMatrixRed[ currentSrcColInt ], xDerivativeMatrixGreen[ currentSrcColInt ],
00133                              xDerivativeMatrixBlue[ currentSrcColInt ], xDerivativeMatrixAlpha[ currentSrcColInt ], xDerivativeMatrixRed[ currentSrcColInt + 1 ], xDerivativeMatrixGreen[ currentSrcColInt + 1 ],
00134                              xDerivativeMatrixBlue[ currentSrcColInt + 1 ], xDerivativeMatrixAlpha[ currentSrcColInt + 1 ] ) );
00135         }
00136         else if ( currentSrcRowInt >= ( srcImage.height() - 1 ) )
00137         {
00138           int idx = ( srcImage.height() - 1 ) * srcImage.width() + currentSrcColInt;
00139           px1 = srcImage.pixel( currentSrcColInt, srcImage.height() - 1 );
00140           px2 = srcImage.pixel( currentSrcColInt + 1, srcImage.height() - 1 );
00141           dstImage.setPixel( j, i, curveInterpolation( px1, px2, u, xDerivativeMatrixRed[ idx ], xDerivativeMatrixGreen[ idx ], xDerivativeMatrixBlue[idx],
00142                              xDerivativeMatrixAlpha[idx], xDerivativeMatrixRed[ idx + 1 ], xDerivativeMatrixGreen[ idx + 1 ], xDerivativeMatrixBlue[idx + 1],
00143                              xDerivativeMatrixAlpha[idx + 1] ) );
00144         }
00145         else if ( currentSrcColInt < 0 )
00146         {
00147           int idx1 = currentSrcRowInt * srcImage.width();
00148           int idx2 = idx1 + srcImage.width();
00149           px1 = srcImage.pixel( 0, currentSrcRowInt );
00150           px2 = srcImage.pixel( 0, currentSrcRowInt + 1 );
00151           dstImage.setPixel( j, i, curveInterpolation( px1, px2, v, yDerivativeMatrixRed[ idx1 ], yDerivativeMatrixGreen[ idx1 ], yDerivativeMatrixBlue[ idx1],
00152                              yDerivativeMatrixAlpha[ idx1], yDerivativeMatrixRed[ idx2 ], yDerivativeMatrixGreen[ idx2 ], yDerivativeMatrixBlue[ idx2],
00153                              yDerivativeMatrixAlpha[ idx2] ) );
00154         }
00155         else if ( currentSrcColInt >= ( srcImage.width() - 1 ) )
00156         {
00157           int idx1 = currentSrcRowInt * srcImage.width() + srcImage.width() - 1;
00158           int idx2 = idx1 + srcImage.width();
00159           px1 = srcImage.pixel( srcImage.width() - 1, currentSrcRowInt );
00160           px2 = srcImage.pixel( srcImage.width() - 1, currentSrcRowInt + 1 );
00161           dstImage.setPixel( j, i, curveInterpolation( px1, px2, v, yDerivativeMatrixRed[ idx1 ], yDerivativeMatrixGreen[ idx1 ], yDerivativeMatrixBlue[ idx1],
00162                              yDerivativeMatrixAlpha[ idx1], yDerivativeMatrixRed[ idx2 ], yDerivativeMatrixGreen[ idx2 ], yDerivativeMatrixBlue[ idx2],
00163                              yDerivativeMatrixAlpha[ idx2] ) );
00164         }
00165         currentSrcCol += nSrcPerDstX;
00166         continue;
00167       }
00168
00169       //first update the control points if necessary
00170       if ( currentSrcColInt != lastSrcColInt || currentSrcRowInt != lastSrcRowInt )
00171       {
00172         calculateControlPoints( nCols, nRows, currentSrcRowInt, currentSrcColInt, redMatrix, greenMatrix, blueMatrix, alphaMatrix,
00173                                 xDerivativeMatrixRed, xDerivativeMatrixGreen, xDerivativeMatrixBlue, xDerivativeMatrixAlpha,
00174                                 yDerivativeMatrixRed, yDerivativeMatrixGreen, yDerivativeMatrixBlue, yDerivativeMatrixAlpha );
00175       }
00176
00177       //bernstein polynomials
00178       bp0u = calcBernsteinPoly( 3, 0, u ); bp1u = calcBernsteinPoly( 3, 1, u );
00179       bp2u = calcBernsteinPoly( 3, 2, u ); bp3u = calcBernsteinPoly( 3, 3, u );
00180       bp0v = calcBernsteinPoly( 3, 0, v ); bp1v = calcBernsteinPoly( 3, 1, v );
00181       bp2v = calcBernsteinPoly( 3, 2, v ); bp3v = calcBernsteinPoly( 3, 3, v );
00182
00183       //then calculate value based on bernstein form of Bezier patch
00184       //todo: move into function
00185       r = bp0u * bp0v * cRed00 +
00186           bp1u * bp0v * cRed10 +
00187           bp2u * bp0v * cRed20 +
00188           bp3u * bp0v * cRed30 +
00189           bp0u * bp1v * cRed01 +
00190           bp1u * bp1v * cRed11 +
00191           bp2u * bp1v * cRed21 +
00192           bp3u * bp1v * cRed31 +
00193           bp0u * bp2v * cRed02 +
00194           bp1u * bp2v * cRed12 +
00195           bp2u * bp2v * cRed22 +
00196           bp3u * bp2v * cRed32 +
00197           bp0u * bp3v * cRed03 +
00198           bp1u * bp3v * cRed13 +
00199           bp2u * bp3v * cRed23 +
00200           bp3u * bp3v * cRed33;
00201
00202       g = bp0u * bp0v * cGreen00 +
00203           bp1u * bp0v * cGreen10 +
00204           bp2u * bp0v * cGreen20 +
00205           bp3u * bp0v * cGreen30 +
00206           bp0u * bp1v * cGreen01 +
00207           bp1u * bp1v * cGreen11 +
00208           bp2u * bp1v * cGreen21 +
00209           bp3u * bp1v * cGreen31 +
00210           bp0u * bp2v * cGreen02 +
00211           bp1u * bp2v * cGreen12 +
00212           bp2u * bp2v * cGreen22 +
00213           bp3u * bp2v * cGreen32 +
00214           bp0u * bp3v * cGreen03 +
00215           bp1u * bp3v * cGreen13 +
00216           bp2u * bp3v * cGreen23 +
00217           bp3u * bp3v * cGreen33;
00218
00219       b = bp0u * bp0v * cBlue00 +
00220           bp1u * bp0v * cBlue10 +
00221           bp2u * bp0v * cBlue20 +
00222           bp3u * bp0v * cBlue30 +
00223           bp0u * bp1v * cBlue01 +
00224           bp1u * bp1v * cBlue11 +
00225           bp2u * bp1v * cBlue21 +
00226           bp3u * bp1v * cBlue31 +
00227           bp0u * bp2v * cBlue02 +
00228           bp1u * bp2v * cBlue12 +
00229           bp2u * bp2v * cBlue22 +
00230           bp3u * bp2v * cBlue32 +
00231           bp0u * bp3v * cBlue03 +
00232           bp1u * bp3v * cBlue13 +
00233           bp2u * bp3v * cBlue23 +
00234           bp3u * bp3v * cBlue33;
00235
00236       a = bp0u * bp0v * cAlpha00 +
00237           bp1u * bp0v * cAlpha10 +
00238           bp2u * bp0v * cAlpha20 +
00239           bp3u * bp0v * cAlpha30 +
00240           bp0u * bp1v * cAlpha01 +
00241           bp1u * bp1v * cAlpha11 +
00242           bp2u * bp1v * cAlpha21 +
00243           bp3u * bp1v * cAlpha31 +
00244           bp0u * bp2v * cAlpha02 +
00245           bp1u * bp2v * cAlpha12 +
00246           bp2u * bp2v * cAlpha22 +
00247           bp3u * bp2v * cAlpha32 +
00248           bp0u * bp3v * cAlpha03 +
00249           bp1u * bp3v * cAlpha13 +
00250           bp2u * bp3v * cAlpha23 +
00251           bp3u * bp3v * cAlpha33;
00252
00253       dstImage.setPixel( j, i, qRgba( r, g, b, a ) );
00254       lastSrcColInt = currentSrcColInt;
00255       currentSrcCol += nSrcPerDstX;
00256     }
00257     lastSrcRowInt = currentSrcRowInt;
00258     currentSrcRow += nSrcPerDstY;
00259   }
00260
00261
00262   //cleanup memory
00263   delete[] redMatrix;
00264   delete[] greenMatrix;
00265   delete[] blueMatrix;
00266   delete[] xDerivativeMatrixRed;
00267   delete[] xDerivativeMatrixGreen;
00268   delete[] xDerivativeMatrixBlue;
00269   delete[] yDerivativeMatrixRed;
00270   delete[] yDerivativeMatrixGreen;
00271   delete[] yDerivativeMatrixBlue;
00272 }
00273
00274 void QgsCubicRasterResampler::xDerivativeMatrix( int nCols, int nRows, double* matrix, const int* colorMatrix )
00275 {
00276   double val = 0;
00277   int index = 0;
00278
00279   for ( int i = 0; i < nRows; ++i )
00280   {
00281     for ( int j = 0; j < nCols; ++j )
00282     {
00283       if ( j < 1 )
00284       {
00285         val = colorMatrix[index + 1] - colorMatrix[index];
00286       }
00287       else if ( j == ( nCols - 1 ) )
00288       {
00289         val = colorMatrix[index] - colorMatrix[ index - 1 ];
00290       }
00291       else
00292       {
00293         val = ( colorMatrix[index + 1] - colorMatrix[index - 1] ) / 2.0;
00294       }
00295       matrix[index] = val;
00296       ++index;
00297     }
00298   }
00299 }
00300
00301 void QgsCubicRasterResampler::yDerivativeMatrix( int nCols, int nRows, double* matrix, const int* colorMatrix )
00302 {
00303   double val = 0;
00304   int index = 0;
00305
00306   for ( int i = 0; i < nRows; ++i )
00307   {
00308     for ( int j = 0; j < nCols; ++j )
00309     {
00310       if ( i == 0 )
00311       {
00312         val = colorMatrix[ index + nCols ] - colorMatrix[ index ];
00313       }
00314       else if ( i == ( nRows - 1 ) )
00315       {
00316         val = colorMatrix[ index ] - colorMatrix[ index - nCols ];
00317       }
00318       else
00319       {
00320         val = ( colorMatrix[ index + nCols ] - colorMatrix[ index - nCols ] ) / 2.0;
00321       }
00322       matrix[index] = val;
00323       ++index;
00324     }
00325   }
00326 }
00327
00328 void QgsCubicRasterResampler::calculateControlPoints( int nCols, int nRows, int currentRow, int currentCol, int* redMatrix, int* greenMatrix, int* blueMatrix,
00329     int* alphaMatrix, double* xDerivativeMatrixRed, double* xDerivativeMatrixGreen, double* xDerivativeMatrixBlue, double* xDerivativeMatrixAlpha,
00330     double* yDerivativeMatrixRed, double* yDerivativeMatrixGreen, double* yDerivativeMatrixBlue, double* yDerivativeMatrixAlpha )
00331 {
00332   Q_UNUSED( nRows );
00333   int idx00 = currentRow * nCols + currentCol;
00334   int idx10 = idx00 + 1;
00335   int idx01 = idx00 + nCols;
00336   int idx11 = idx01 + 1;
00337
00338   //corner points
00339   cRed00 = redMatrix[idx00]; cGreen00 = greenMatrix[idx00]; cBlue00 = blueMatrix[idx00]; cAlpha00 = alphaMatrix[idx00];
00340   cRed30 = redMatrix[idx10]; cGreen30 = greenMatrix[idx10]; cBlue30 = blueMatrix[idx10]; cAlpha30 = alphaMatrix[idx10];
00341   cRed03 = redMatrix[idx01]; cGreen03 = greenMatrix[idx01]; cBlue03 = blueMatrix[idx01]; cAlpha03 = alphaMatrix[idx01];
00342   cRed33 = redMatrix[idx11]; cGreen33 = greenMatrix[idx11]; cBlue33 = blueMatrix[idx11]; cAlpha33 = alphaMatrix[idx11];
00343
00344   //control points near c00
00345   cRed10 = cRed00 + 0.333 * xDerivativeMatrixRed[idx00]; cGreen10 = cGreen00 + 0.333 * xDerivativeMatrixGreen[idx00];
00346   cBlue10 = cBlue00 + 0.333 * xDerivativeMatrixBlue[idx00];cAlpha10 = cAlpha00 + 0.333 * xDerivativeMatrixAlpha[idx00];
00347   cRed01 = cRed00 + 0.333 * yDerivativeMatrixRed[idx00]; cGreen01 = cGreen00 + 0.333 * yDerivativeMatrixGreen[idx00];
00348   cBlue01 = cBlue00 + 0.333 * yDerivativeMatrixBlue[idx00];cAlpha01 = cAlpha00 + 0.333 * yDerivativeMatrixAlpha[idx00];
00349   cRed11 = cRed10 + 0.333 * yDerivativeMatrixRed[idx00]; cGreen11 = cGreen10 + 0.333 * yDerivativeMatrixGreen[idx00];
00350   cBlue11 = cBlue10 + 0.333 * yDerivativeMatrixBlue[idx00];cAlpha11 = cAlpha10 + 0.333 * yDerivativeMatrixAlpha[idx00];
00351
00352   //control points near c30
00353   cRed20 = cRed30 - 0.333 * xDerivativeMatrixRed[idx10]; cGreen20 = cGreen30 - 0.333 * xDerivativeMatrixGreen[idx10];
00354   cBlue20 = cBlue30 - 0.333 * xDerivativeMatrixBlue[idx10]; cAlpha20 = cAlpha30 - 0.333 * xDerivativeMatrixAlpha[idx10];
00355   cRed31 = cRed30 + 0.333 * yDerivativeMatrixRed[idx10]; cGreen31 = cGreen30 + 0.333 * yDerivativeMatrixGreen[idx10];
00356   cBlue31 = cBlue30 + 0.333 * yDerivativeMatrixBlue[idx10]; cAlpha31 = cAlpha30 + 0.333 * yDerivativeMatrixAlpha[idx10];
00357   cRed21 = cRed20 + 0.333 * yDerivativeMatrixRed[idx10]; cGreen21 = cGreen20 + 0.333 * yDerivativeMatrixGreen[idx10];
00358   cBlue21 = cBlue20 + 0.333 * yDerivativeMatrixBlue[idx10]; cAlpha21 = cAlpha20 + 0.333 * yDerivativeMatrixAlpha[idx10];
00359
00360   //control points near c03
00361   cRed13 = cRed03 + 0.333 * xDerivativeMatrixRed[idx01]; cGreen13 = cGreen03 + 0.333 * xDerivativeMatrixGreen[idx01];
00362   cBlue13 = cBlue03 + 0.333 * xDerivativeMatrixBlue[idx01]; cAlpha13 = cAlpha03 + 0.333 * xDerivativeMatrixAlpha[idx01];
00363   cRed02 = cRed03 - 0.333 * yDerivativeMatrixRed[idx01]; cGreen02 = cGreen03 - 0.333 * yDerivativeMatrixGreen[idx01];
00364   cBlue02 = cBlue03 - 0.333 * yDerivativeMatrixBlue[idx01]; cAlpha02 = cAlpha03 - 0.333 * yDerivativeMatrixAlpha[idx01];
00365   cRed12 = cRed02 + 0.333 * xDerivativeMatrixRed[idx01]; cGreen12 = cGreen02 + 0.333 * xDerivativeMatrixGreen[idx01];
00366   cBlue12 = cBlue02 + 0.333 * xDerivativeMatrixBlue[idx01]; cAlpha12 = cAlpha02 + 0.333 * xDerivativeMatrixAlpha[idx01];
00367
00368   //control points near c33
00369   cRed23 = cRed33 - 0.333 * xDerivativeMatrixRed[idx11]; cGreen23 = cGreen33 - 0.333 * xDerivativeMatrixGreen[idx11];
00370   cBlue23 = cBlue33 - 0.333 * xDerivativeMatrixBlue[idx11]; cAlpha23 = cAlpha33 - 0.333 * xDerivativeMatrixAlpha[idx11];
00371   cRed32 = cRed33 - 0.333 * yDerivativeMatrixRed[idx11]; cGreen32 = cGreen33 - 0.333 * yDerivativeMatrixGreen[idx11];
00372   cBlue32 = cBlue33 - 0.333 * yDerivativeMatrixBlue[idx11]; cAlpha32 = cAlpha33 - 0.333 * yDerivativeMatrixAlpha[idx11];
00373   cRed22 = cRed32 - 0.333 * xDerivativeMatrixRed[idx11]; cGreen22 = cGreen32 - 0.333 * xDerivativeMatrixGreen[idx11];
00374   cBlue22 = cBlue32 - 0.333 * xDerivativeMatrixBlue[idx11]; cAlpha22 = cAlpha32 - 0.333 * xDerivativeMatrixAlpha[idx11];
00375 }
00376
00377 QRgb QgsCubicRasterResampler::curveInterpolation( QRgb pt1, QRgb pt2, double t, double d1red, double d1green, double d1blue, double d1alpha,
00378     double d2red, double d2green, double d2blue, double d2alpha )
00379 {
00380   //control points
00381   double p0r = qRed( pt1 ); double p1r = p0r + 0.333 * d1red; double p3r = qRed( pt2 ); double p2r = p3r - 0.333 * d2red;
00382   double p0g = qGreen( pt1 ); double p1g = p0g + 0.333 * d1green; double p3g = qGreen( pt2 ); double p2g = p3g - 0.333 * d2green;
00383   double p0b = qBlue( pt1 ); double p1b = p0b + 0.333 * d1blue; double p3b = qBlue( pt2 ); double p2b = p3b - 0.333 * d2blue;
00384   double p0a = qAlpha( pt1 ); double p1a = p0a + 0.333 * d1alpha; double p3a = qAlpha( pt2 ); double p2a = p3a - 0.333 * d2alpha;
00385
00386   //bernstein polynomials
00387   double bp0 = calcBernsteinPoly( 3, 0, t );
00388   double bp1 = calcBernsteinPoly( 3, 1, t );
00389   double bp2 = calcBernsteinPoly( 3, 2, t );
00390   double bp3 = calcBernsteinPoly( 3, 3, t );
00391
00392   int red = bp0 * p0r + bp1 * p1r + bp2 * p2r + bp3 * p3r;
00393   int green = bp0 * p0g + bp1 * p1g + bp2 * p2g + bp3 * p3g;
00394   int blue = bp0 * p0b + bp1 * p1b + bp2 * p2b + bp3 * p3b;
00395   int alpha = bp0 * p0a + bp1 * p1a + bp2 * p2a + bp3 * p3a;
00396
00397   return qRgba( red, green, blue, alpha );
00398 }
00399
00400 double QgsCubicRasterResampler::calcBernsteinPoly( int n, int i, double t )
00401 {
00402   if ( i < 0 )
00403   {
00404     return 0;
00405   }
00406
00407   return lower( n, i )*power( t, i )*power(( 1 - t ), ( n - i ) );
00408 }
00409
00410 int QgsCubicRasterResampler::lower( int n, int i )
00411 {
00412   if ( i >= 0 && i <= n )
00413   {
00414     return faculty( n ) / ( faculty( i )*faculty( n - i ) );
00415   }
00416   else
00417   {
00418     return 0;
00419   }
00420 }
00421
00422 double QgsCubicRasterResampler::power( double a, int b )
00423 {
00424   if ( b == 0 )
00425   {
00426     return 1;
00427   }
00428   double tmp = a;
00429   for ( int i = 2; i <= qAbs(( double )b ); i++ )
00430   {
00431
00432     a *= tmp;
00433   }
00434   if ( b > 0 )
00435   {
00436     return a;
00437   }
00438   else
00439   {
00440     return ( 1.0 / a );
00441   }
00442 }
00443
00444 int QgsCubicRasterResampler::faculty( int n )
00445 {
00446   if ( n < 0 )//Is faculty also defined for negative integers?
00447   {
00448     return 0;
00449   }
00450   int i;
00451   int result = n;
00452
00453   if ( n == 0 || n == 1 )
00454     {return 1;}//faculty of 0 is 1!
00455
00456   for ( i = n - 1; i >= 2; i-- )
00457   {
00458     result *= i;
00459   }
00460   return result;
00461 }