QGIS API Documentation  2.2.0-Valmiera
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgscubicrasterresampler.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgscubicrasterresampler.cpp
3  ----------------------------
4  begin : December 2011
5  copyright : (C) 2011 by Marco Hugentobler
6  email : marco at sourcepole dot ch
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgscubicrasterresampler.h"
19 #include <QImage>
20 #include <cmath>
21 
22 QgsCubicRasterResampler::QgsCubicRasterResampler()
23 {
24 }
25 
26 QgsCubicRasterResampler::~QgsCubicRasterResampler()
27 {
28 }
29 
30 QgsRasterResampler * QgsCubicRasterResampler::clone() const
31 {
32  return new QgsCubicRasterResampler();
33 }
34 
35 void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage )
36 {
37  int nCols = srcImage.width();
38  int nRows = srcImage.height();
39 
40  int pos = 0;
41  QRgb px;
42  int* redMatrix = new int[ nCols * nRows ];
43  int* greenMatrix = new int[ nCols * nRows ];
44  int* blueMatrix = new int[ nCols * nRows ];
45  int* alphaMatrix = new int[ nCols * nRows ];
46 
47  for ( int i = 0; i < nRows; ++i )
48  {
49  for ( int j = 0; j < nCols; ++j )
50  {
51  px = srcImage.pixel( j, i );
52  redMatrix[pos] = qRed( px );
53  greenMatrix[pos] = qGreen( px );
54  blueMatrix[pos] = qBlue( px );
55  alphaMatrix[pos] = qAlpha( px );
56  ++pos;
57  }
58  }
59 
60  //derivative x
61  double* xDerivativeMatrixRed = new double[ nCols * nRows ];
62  xDerivativeMatrix( nCols, nRows, xDerivativeMatrixRed, redMatrix );
63  double* xDerivativeMatrixGreen = new double[ nCols * nRows ];
64  xDerivativeMatrix( nCols, nRows, xDerivativeMatrixGreen, greenMatrix );
65  double* xDerivativeMatrixBlue = new double[ nCols * nRows ];
66  xDerivativeMatrix( nCols, nRows, xDerivativeMatrixBlue, blueMatrix );
67  double* xDerivativeMatrixAlpha = new double[ nCols * nRows ];
68  xDerivativeMatrix( nCols, nRows, xDerivativeMatrixAlpha, alphaMatrix );
69 
70  //derivative y
71  double* yDerivativeMatrixRed = new double[ nCols * nRows ];
72  yDerivativeMatrix( nCols, nRows, yDerivativeMatrixRed, redMatrix );
73  double* yDerivativeMatrixGreen = new double[ nCols * nRows ];
74  yDerivativeMatrix( nCols, nRows, yDerivativeMatrixGreen, greenMatrix );
75  double* yDerivativeMatrixBlue = new double[ nCols * nRows ];
76  yDerivativeMatrix( nCols, nRows, yDerivativeMatrixBlue, blueMatrix );
77  double* yDerivativeMatrixAlpha = new double[ nCols * nRows ];
78  yDerivativeMatrix( nCols, nRows, yDerivativeMatrixAlpha, alphaMatrix );
79 
80  //compute output
81  double nSrcPerDstX = ( double ) srcImage.width() / ( double ) dstImage.width();
82  double nSrcPerDstY = ( double ) srcImage.height() / ( double ) dstImage.height();
83 
84  double currentSrcRow = nSrcPerDstY / 2.0 - 0.5;
85  double currentSrcCol;
86  int currentSrcColInt;
87  int currentSrcRowInt;
88  int lastSrcColInt = -1;
89  int lastSrcRowInt = -1;
90 
91  int r, g, b, a;
92  //bernstein polynomials
93  double bp0u, bp1u, bp2u, bp3u, bp0v, bp1v, bp2v, bp3v;
94  double u, v;
95 
96  for ( int i = 0; i < dstImage.height(); ++i )
97  {
98  currentSrcRowInt = floor( currentSrcRow );
99  v = currentSrcRow - currentSrcRowInt;
100 
101  currentSrcCol = nSrcPerDstX / 2.0 - 0.5;
102  for ( int j = 0; j < dstImage.width(); ++j )
103  {
104  currentSrcColInt = floor( currentSrcCol );
105  u = currentSrcCol - currentSrcColInt;
106 
107  //handle eight edge-cases
108  if ( currentSrcRowInt < 0 || currentSrcRowInt >= ( srcImage.height() - 1 ) || currentSrcColInt < 0 || currentSrcColInt >= ( srcImage.width() - 1 ) )
109  {
110  QRgb px1, px2;
111  //pixels at the border of the source image needs to be handled in a special way
112  if ( currentSrcRowInt < 0 && currentSrcColInt < 0 )
113  {
114  dstImage.setPixel( j, i, srcImage.pixel( 0, 0 ) );
115  }
116  else if ( currentSrcRowInt < 0 && currentSrcColInt >= ( srcImage.width() - 1 ) )
117  {
118  dstImage.setPixel( j, i, srcImage.pixel( srcImage.width() - 1, 0 ) );
119  }
120  else if ( currentSrcRowInt >= ( srcImage.height() - 1 ) && currentSrcColInt >= ( srcImage.width() - 1 ) )
121  {
122  dstImage.setPixel( j, i, srcImage.pixel( srcImage.width() - 1, srcImage.height() - 1 ) );
123  }
124  else if ( currentSrcRowInt >= ( srcImage.height() - 1 ) && currentSrcColInt < 0 )
125  {
126  dstImage.setPixel( j, i, srcImage.pixel( 0, srcImage.height() - 1 ) );
127  }
128  else if ( currentSrcRowInt < 0 )
129  {
130  px1 = srcImage.pixel( currentSrcColInt, 0 );
131  px2 = srcImage.pixel( currentSrcColInt + 1, 0 );
132  dstImage.setPixel( j, i, curveInterpolation( px1, px2, u, xDerivativeMatrixRed[ currentSrcColInt ], xDerivativeMatrixGreen[ currentSrcColInt ],
133  xDerivativeMatrixBlue[ currentSrcColInt ], xDerivativeMatrixAlpha[ currentSrcColInt ], xDerivativeMatrixRed[ currentSrcColInt + 1 ], xDerivativeMatrixGreen[ currentSrcColInt + 1 ],
134  xDerivativeMatrixBlue[ currentSrcColInt + 1 ], xDerivativeMatrixAlpha[ currentSrcColInt + 1 ] ) );
135  }
136  else if ( currentSrcRowInt >= ( srcImage.height() - 1 ) )
137  {
138  int idx = ( srcImage.height() - 1 ) * srcImage.width() + currentSrcColInt;
139  px1 = srcImage.pixel( currentSrcColInt, srcImage.height() - 1 );
140  px2 = srcImage.pixel( currentSrcColInt + 1, srcImage.height() - 1 );
141  dstImage.setPixel( j, i, curveInterpolation( px1, px2, u, xDerivativeMatrixRed[ idx ], xDerivativeMatrixGreen[ idx ], xDerivativeMatrixBlue[idx],
142  xDerivativeMatrixAlpha[idx], xDerivativeMatrixRed[ idx + 1 ], xDerivativeMatrixGreen[ idx + 1 ], xDerivativeMatrixBlue[idx + 1],
143  xDerivativeMatrixAlpha[idx + 1] ) );
144  }
145  else if ( currentSrcColInt < 0 )
146  {
147  int idx1 = currentSrcRowInt * srcImage.width();
148  int idx2 = idx1 + srcImage.width();
149  px1 = srcImage.pixel( 0, currentSrcRowInt );
150  px2 = srcImage.pixel( 0, currentSrcRowInt + 1 );
151  dstImage.setPixel( j, i, curveInterpolation( px1, px2, v, yDerivativeMatrixRed[ idx1 ], yDerivativeMatrixGreen[ idx1 ], yDerivativeMatrixBlue[ idx1],
152  yDerivativeMatrixAlpha[ idx1], yDerivativeMatrixRed[ idx2 ], yDerivativeMatrixGreen[ idx2 ], yDerivativeMatrixBlue[ idx2],
153  yDerivativeMatrixAlpha[ idx2] ) );
154  }
155  else if ( currentSrcColInt >= ( srcImage.width() - 1 ) )
156  {
157  int idx1 = currentSrcRowInt * srcImage.width() + srcImage.width() - 1;
158  int idx2 = idx1 + srcImage.width();
159  px1 = srcImage.pixel( srcImage.width() - 1, currentSrcRowInt );
160  px2 = srcImage.pixel( srcImage.width() - 1, currentSrcRowInt + 1 );
161  dstImage.setPixel( j, i, curveInterpolation( px1, px2, v, yDerivativeMatrixRed[ idx1 ], yDerivativeMatrixGreen[ idx1 ], yDerivativeMatrixBlue[ idx1],
162  yDerivativeMatrixAlpha[ idx1], yDerivativeMatrixRed[ idx2 ], yDerivativeMatrixGreen[ idx2 ], yDerivativeMatrixBlue[ idx2],
163  yDerivativeMatrixAlpha[ idx2] ) );
164  }
165  currentSrcCol += nSrcPerDstX;
166  continue;
167  }
168 
169  //first update the control points if necessary
170  if ( currentSrcColInt != lastSrcColInt || currentSrcRowInt != lastSrcRowInt )
171  {
172  calculateControlPoints( nCols, nRows, currentSrcRowInt, currentSrcColInt, redMatrix, greenMatrix, blueMatrix, alphaMatrix,
173  xDerivativeMatrixRed, xDerivativeMatrixGreen, xDerivativeMatrixBlue, xDerivativeMatrixAlpha,
174  yDerivativeMatrixRed, yDerivativeMatrixGreen, yDerivativeMatrixBlue, yDerivativeMatrixAlpha );
175  }
176 
177  //bernstein polynomials
178  bp0u = calcBernsteinPoly( 3, 0, u ); bp1u = calcBernsteinPoly( 3, 1, u );
179  bp2u = calcBernsteinPoly( 3, 2, u ); bp3u = calcBernsteinPoly( 3, 3, u );
180  bp0v = calcBernsteinPoly( 3, 0, v ); bp1v = calcBernsteinPoly( 3, 1, v );
181  bp2v = calcBernsteinPoly( 3, 2, v ); bp3v = calcBernsteinPoly( 3, 3, v );
182 
183  //then calculate value based on bernstein form of Bezier patch
184  //todo: move into function
185  r = bp0u * bp0v * cRed00 +
186  bp1u * bp0v * cRed10 +
187  bp2u * bp0v * cRed20 +
188  bp3u * bp0v * cRed30 +
189  bp0u * bp1v * cRed01 +
190  bp1u * bp1v * cRed11 +
191  bp2u * bp1v * cRed21 +
192  bp3u * bp1v * cRed31 +
193  bp0u * bp2v * cRed02 +
194  bp1u * bp2v * cRed12 +
195  bp2u * bp2v * cRed22 +
196  bp3u * bp2v * cRed32 +
197  bp0u * bp3v * cRed03 +
198  bp1u * bp3v * cRed13 +
199  bp2u * bp3v * cRed23 +
200  bp3u * bp3v * cRed33;
201 
202  g = bp0u * bp0v * cGreen00 +
203  bp1u * bp0v * cGreen10 +
204  bp2u * bp0v * cGreen20 +
205  bp3u * bp0v * cGreen30 +
206  bp0u * bp1v * cGreen01 +
207  bp1u * bp1v * cGreen11 +
208  bp2u * bp1v * cGreen21 +
209  bp3u * bp1v * cGreen31 +
210  bp0u * bp2v * cGreen02 +
211  bp1u * bp2v * cGreen12 +
212  bp2u * bp2v * cGreen22 +
213  bp3u * bp2v * cGreen32 +
214  bp0u * bp3v * cGreen03 +
215  bp1u * bp3v * cGreen13 +
216  bp2u * bp3v * cGreen23 +
217  bp3u * bp3v * cGreen33;
218 
219  b = bp0u * bp0v * cBlue00 +
220  bp1u * bp0v * cBlue10 +
221  bp2u * bp0v * cBlue20 +
222  bp3u * bp0v * cBlue30 +
223  bp0u * bp1v * cBlue01 +
224  bp1u * bp1v * cBlue11 +
225  bp2u * bp1v * cBlue21 +
226  bp3u * bp1v * cBlue31 +
227  bp0u * bp2v * cBlue02 +
228  bp1u * bp2v * cBlue12 +
229  bp2u * bp2v * cBlue22 +
230  bp3u * bp2v * cBlue32 +
231  bp0u * bp3v * cBlue03 +
232  bp1u * bp3v * cBlue13 +
233  bp2u * bp3v * cBlue23 +
234  bp3u * bp3v * cBlue33;
235 
236  a = bp0u * bp0v * cAlpha00 +
237  bp1u * bp0v * cAlpha10 +
238  bp2u * bp0v * cAlpha20 +
239  bp3u * bp0v * cAlpha30 +
240  bp0u * bp1v * cAlpha01 +
241  bp1u * bp1v * cAlpha11 +
242  bp2u * bp1v * cAlpha21 +
243  bp3u * bp1v * cAlpha31 +
244  bp0u * bp2v * cAlpha02 +
245  bp1u * bp2v * cAlpha12 +
246  bp2u * bp2v * cAlpha22 +
247  bp3u * bp2v * cAlpha32 +
248  bp0u * bp3v * cAlpha03 +
249  bp1u * bp3v * cAlpha13 +
250  bp2u * bp3v * cAlpha23 +
251  bp3u * bp3v * cAlpha33;
252 
253  dstImage.setPixel( j, i, qRgba( r, g, b, a ) );
254  lastSrcColInt = currentSrcColInt;
255  currentSrcCol += nSrcPerDstX;
256  }
257  lastSrcRowInt = currentSrcRowInt;
258  currentSrcRow += nSrcPerDstY;
259  }
260 
261 
262  //cleanup memory
263  delete[] redMatrix;
264  delete[] greenMatrix;
265  delete[] blueMatrix;
266  delete[] xDerivativeMatrixRed;
267  delete[] xDerivativeMatrixGreen;
268  delete[] xDerivativeMatrixBlue;
269  delete[] yDerivativeMatrixRed;
270  delete[] yDerivativeMatrixGreen;
271  delete[] yDerivativeMatrixBlue;
272 }
273 
274 void QgsCubicRasterResampler::xDerivativeMatrix( int nCols, int nRows, double* matrix, const int* colorMatrix )
275 {
276  double val = 0;
277  int index = 0;
278 
279  for ( int i = 0; i < nRows; ++i )
280  {
281  for ( int j = 0; j < nCols; ++j )
282  {
283  if ( j < 1 )
284  {
285  val = colorMatrix[index + 1] - colorMatrix[index];
286  }
287  else if ( j == ( nCols - 1 ) )
288  {
289  val = colorMatrix[index] - colorMatrix[ index - 1 ];
290  }
291  else
292  {
293  val = ( colorMatrix[index + 1] - colorMatrix[index - 1] ) / 2.0;
294  }
295  matrix[index] = val;
296  ++index;
297  }
298  }
299 }
300 
301 void QgsCubicRasterResampler::yDerivativeMatrix( int nCols, int nRows, double* matrix, const int* colorMatrix )
302 {
303  double val = 0;
304  int index = 0;
305 
306  for ( int i = 0; i < nRows; ++i )
307  {
308  for ( int j = 0; j < nCols; ++j )
309  {
310  if ( i == 0 )
311  {
312  val = colorMatrix[ index + nCols ] - colorMatrix[ index ];
313  }
314  else if ( i == ( nRows - 1 ) )
315  {
316  val = colorMatrix[ index ] - colorMatrix[ index - nCols ];
317  }
318  else
319  {
320  val = ( colorMatrix[ index + nCols ] - colorMatrix[ index - nCols ] ) / 2.0;
321  }
322  matrix[index] = val;
323  ++index;
324  }
325  }
326 }
327 
328 void QgsCubicRasterResampler::calculateControlPoints( int nCols, int nRows, int currentRow, int currentCol, int* redMatrix, int* greenMatrix, int* blueMatrix,
329  int* alphaMatrix, double* xDerivativeMatrixRed, double* xDerivativeMatrixGreen, double* xDerivativeMatrixBlue, double* xDerivativeMatrixAlpha,
330  double* yDerivativeMatrixRed, double* yDerivativeMatrixGreen, double* yDerivativeMatrixBlue, double* yDerivativeMatrixAlpha )
331 {
332  Q_UNUSED( nRows );
333  int idx00 = currentRow * nCols + currentCol;
334  int idx10 = idx00 + 1;
335  int idx01 = idx00 + nCols;
336  int idx11 = idx01 + 1;
337 
338  //corner points
339  cRed00 = redMatrix[idx00]; cGreen00 = greenMatrix[idx00]; cBlue00 = blueMatrix[idx00]; cAlpha00 = alphaMatrix[idx00];
340  cRed30 = redMatrix[idx10]; cGreen30 = greenMatrix[idx10]; cBlue30 = blueMatrix[idx10]; cAlpha30 = alphaMatrix[idx10];
341  cRed03 = redMatrix[idx01]; cGreen03 = greenMatrix[idx01]; cBlue03 = blueMatrix[idx01]; cAlpha03 = alphaMatrix[idx01];
342  cRed33 = redMatrix[idx11]; cGreen33 = greenMatrix[idx11]; cBlue33 = blueMatrix[idx11]; cAlpha33 = alphaMatrix[idx11];
343 
344  //control points near c00
345  cRed10 = cRed00 + 0.333 * xDerivativeMatrixRed[idx00]; cGreen10 = cGreen00 + 0.333 * xDerivativeMatrixGreen[idx00];
346  cBlue10 = cBlue00 + 0.333 * xDerivativeMatrixBlue[idx00];cAlpha10 = cAlpha00 + 0.333 * xDerivativeMatrixAlpha[idx00];
347  cRed01 = cRed00 + 0.333 * yDerivativeMatrixRed[idx00]; cGreen01 = cGreen00 + 0.333 * yDerivativeMatrixGreen[idx00];
348  cBlue01 = cBlue00 + 0.333 * yDerivativeMatrixBlue[idx00];cAlpha01 = cAlpha00 + 0.333 * yDerivativeMatrixAlpha[idx00];
349  cRed11 = cRed10 + 0.333 * yDerivativeMatrixRed[idx00]; cGreen11 = cGreen10 + 0.333 * yDerivativeMatrixGreen[idx00];
350  cBlue11 = cBlue10 + 0.333 * yDerivativeMatrixBlue[idx00];cAlpha11 = cAlpha10 + 0.333 * yDerivativeMatrixAlpha[idx00];
351 
352  //control points near c30
353  cRed20 = cRed30 - 0.333 * xDerivativeMatrixRed[idx10]; cGreen20 = cGreen30 - 0.333 * xDerivativeMatrixGreen[idx10];
354  cBlue20 = cBlue30 - 0.333 * xDerivativeMatrixBlue[idx10]; cAlpha20 = cAlpha30 - 0.333 * xDerivativeMatrixAlpha[idx10];
355  cRed31 = cRed30 + 0.333 * yDerivativeMatrixRed[idx10]; cGreen31 = cGreen30 + 0.333 * yDerivativeMatrixGreen[idx10];
356  cBlue31 = cBlue30 + 0.333 * yDerivativeMatrixBlue[idx10]; cAlpha31 = cAlpha30 + 0.333 * yDerivativeMatrixAlpha[idx10];
357  cRed21 = cRed20 + 0.333 * yDerivativeMatrixRed[idx10]; cGreen21 = cGreen20 + 0.333 * yDerivativeMatrixGreen[idx10];
358  cBlue21 = cBlue20 + 0.333 * yDerivativeMatrixBlue[idx10]; cAlpha21 = cAlpha20 + 0.333 * yDerivativeMatrixAlpha[idx10];
359 
360  //control points near c03
361  cRed13 = cRed03 + 0.333 * xDerivativeMatrixRed[idx01]; cGreen13 = cGreen03 + 0.333 * xDerivativeMatrixGreen[idx01];
362  cBlue13 = cBlue03 + 0.333 * xDerivativeMatrixBlue[idx01]; cAlpha13 = cAlpha03 + 0.333 * xDerivativeMatrixAlpha[idx01];
363  cRed02 = cRed03 - 0.333 * yDerivativeMatrixRed[idx01]; cGreen02 = cGreen03 - 0.333 * yDerivativeMatrixGreen[idx01];
364  cBlue02 = cBlue03 - 0.333 * yDerivativeMatrixBlue[idx01]; cAlpha02 = cAlpha03 - 0.333 * yDerivativeMatrixAlpha[idx01];
365  cRed12 = cRed02 + 0.333 * xDerivativeMatrixRed[idx01]; cGreen12 = cGreen02 + 0.333 * xDerivativeMatrixGreen[idx01];
366  cBlue12 = cBlue02 + 0.333 * xDerivativeMatrixBlue[idx01]; cAlpha12 = cAlpha02 + 0.333 * xDerivativeMatrixAlpha[idx01];
367 
368  //control points near c33
369  cRed23 = cRed33 - 0.333 * xDerivativeMatrixRed[idx11]; cGreen23 = cGreen33 - 0.333 * xDerivativeMatrixGreen[idx11];
370  cBlue23 = cBlue33 - 0.333 * xDerivativeMatrixBlue[idx11]; cAlpha23 = cAlpha33 - 0.333 * xDerivativeMatrixAlpha[idx11];
371  cRed32 = cRed33 - 0.333 * yDerivativeMatrixRed[idx11]; cGreen32 = cGreen33 - 0.333 * yDerivativeMatrixGreen[idx11];
372  cBlue32 = cBlue33 - 0.333 * yDerivativeMatrixBlue[idx11]; cAlpha32 = cAlpha33 - 0.333 * yDerivativeMatrixAlpha[idx11];
373  cRed22 = cRed32 - 0.333 * xDerivativeMatrixRed[idx11]; cGreen22 = cGreen32 - 0.333 * xDerivativeMatrixGreen[idx11];
374  cBlue22 = cBlue32 - 0.333 * xDerivativeMatrixBlue[idx11]; cAlpha22 = cAlpha32 - 0.333 * xDerivativeMatrixAlpha[idx11];
375 }
376 
377 QRgb QgsCubicRasterResampler::curveInterpolation( QRgb pt1, QRgb pt2, double t, double d1red, double d1green, double d1blue, double d1alpha,
378  double d2red, double d2green, double d2blue, double d2alpha )
379 {
380  //control points
381  double p0r = qRed( pt1 ); double p1r = p0r + 0.333 * d1red; double p3r = qRed( pt2 ); double p2r = p3r - 0.333 * d2red;
382  double p0g = qGreen( pt1 ); double p1g = p0g + 0.333 * d1green; double p3g = qGreen( pt2 ); double p2g = p3g - 0.333 * d2green;
383  double p0b = qBlue( pt1 ); double p1b = p0b + 0.333 * d1blue; double p3b = qBlue( pt2 ); double p2b = p3b - 0.333 * d2blue;
384  double p0a = qAlpha( pt1 ); double p1a = p0a + 0.333 * d1alpha; double p3a = qAlpha( pt2 ); double p2a = p3a - 0.333 * d2alpha;
385 
386  //bernstein polynomials
387  double bp0 = calcBernsteinPoly( 3, 0, t );
388  double bp1 = calcBernsteinPoly( 3, 1, t );
389  double bp2 = calcBernsteinPoly( 3, 2, t );
390  double bp3 = calcBernsteinPoly( 3, 3, t );
391 
392  int red = bp0 * p0r + bp1 * p1r + bp2 * p2r + bp3 * p3r;
393  int green = bp0 * p0g + bp1 * p1g + bp2 * p2g + bp3 * p3g;
394  int blue = bp0 * p0b + bp1 * p1b + bp2 * p2b + bp3 * p3b;
395  int alpha = bp0 * p0a + bp1 * p1a + bp2 * p2a + bp3 * p3a;
396 
397  return qRgba( red, green, blue, alpha );
398 }
399 
400 double QgsCubicRasterResampler::calcBernsteinPoly( int n, int i, double t )
401 {
402  if ( i < 0 )
403  {
404  return 0;
405  }
406 
407  return lower( n, i )*power( t, i )*power(( 1 - t ), ( n - i ) );
408 }
409 
410 int QgsCubicRasterResampler::lower( int n, int i )
411 {
412  if ( i >= 0 && i <= n )
413  {
414  return faculty( n ) / ( faculty( i )*faculty( n - i ) );
415  }
416  else
417  {
418  return 0;
419  }
420 }
421 
422 double QgsCubicRasterResampler::power( double a, int b )
423 {
424  if ( b == 0 )
425  {
426  return 1;
427  }
428  double tmp = a;
429  for ( int i = 2; i <= qAbs(( double )b ); i++ )
430  {
431 
432  a *= tmp;
433  }
434  if ( b > 0 )
435  {
436  return a;
437  }
438  else
439  {
440  return ( 1.0 / a );
441  }
442 }
443 
444 int QgsCubicRasterResampler::faculty( int n )
445 {
446  if ( n < 0 )//Is faculty also defined for negative integers?
447  {
448  return 0;
449  }
450  int i;
451  int result = n;
452 
453  if ( n == 0 || n == 1 )
454  {return 1;}//faculty of 0 is 1!
455 
456  for ( i = n - 1; i >= 2; i-- )
457  {
458  result *= i;
459  }
460  return result;
461 }
Generated on Sun Feb 23 2014 14:27:19 for QGIS API Documentation by   doxygen 1.8.1.2