Use bezier interpolation for cubic resampling. Still buggy

qgis · Dec 21, 2011 · 668b73d · 668b73d
1 parent 3d0b85a
commit 668b73d
Show file tree

Hide file tree

Showing 6 changed files with 243 additions and 99 deletions.
diff --git a/src/core/raster/qgsbilinearrasterresampler.cpp b/src/core/raster/qgsbilinearrasterresampler.cpp
@@ -26,7 +26,7 @@ QgsBilinearRasterResampler::~QgsBilinearRasterResampler()
 {
 }
 
-void QgsBilinearRasterResampler::resample( const QImage& srcImage, QImage& dstImage ) const
+void QgsBilinearRasterResampler::resample( const QImage& srcImage, QImage& dstImage )
 {
   double nSrcPerDstX = ( double ) srcImage.width() / ( double ) dstImage.width();
   double nSrcPerDstY = ( double ) srcImage.height() / ( double ) dstImage.height();

diff --git a/src/core/raster/qgsbilinearrasterresampler.h b/src/core/raster/qgsbilinearrasterresampler.h
@@ -27,7 +27,7 @@ class QgsBilinearRasterResampler: public QgsRasterResampler
     QgsBilinearRasterResampler();
     ~QgsBilinearRasterResampler();
 
-    void resample( const QImage& srcImage, QImage& dstImage ) const;
+    void resample( const QImage& srcImage, QImage& dstImage );
 
   private:
     QRgb resampleColorValue( double u, double v, QRgb col1, QRgb col2, QRgb col3, QRgb col4 ) const;

diff --git a/src/core/raster/qgscubicrasterresampler.cpp b/src/core/raster/qgscubicrasterresampler.cpp
@@ -26,7 +26,7 @@ QgsCubicRasterResampler::~QgsCubicRasterResampler()
 {
 }
 
-void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage ) const
+void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage )
 {
   int nCols = srcImage.width();
   int nRows = srcImage.height();
@@ -37,9 +37,9 @@ void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage
   int* greenMatrix = new int[ nCols * nRows ];
   int* blueMatrix = new int[ nCols * nRows ];
 
-  for( int i = 0; i < nRows; ++i )
+  for ( int i = 0; i < nRows; ++i )
   {
-    for( int j = 0; j < nCols; ++j )
+    for ( int j = 0; j < nCols; ++j )
     {
       px = srcImage.pixel( j, i );
       redMatrix[pos] = qRed( px );
@@ -59,19 +59,11 @@ void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage
 
   //derivative y
   double* yDerivativeMatrixRed = new double[ nCols * nRows ];
-  yDerivativeMatrix<int>( nCols, nRows, yDerivativeMatrixRed, redMatrix );
+  yDerivativeMatrix( nCols, nRows, yDerivativeMatrixRed, redMatrix );
   double* yDerivativeMatrixGreen = new double[ nCols * nRows ];
-  yDerivativeMatrix<int>( nCols, nRows, yDerivativeMatrixGreen, greenMatrix );
+  yDerivativeMatrix( nCols, nRows, yDerivativeMatrixGreen, greenMatrix );
   double* yDerivativeMatrixBlue = new double[ nCols * nRows ];
-  yDerivativeMatrix<int>( nCols, nRows, yDerivativeMatrixBlue, blueMatrix );
-
-  //derivative xy
-  double* xyDerivativeMatrixRed = new double[ nCols * nRows ];
-  yDerivativeMatrix<double>( nCols, nRows, xyDerivativeMatrixRed, xDerivativeMatrixRed );
-  double* xyDerivativeMatrixGreen = new double[ nCols * nRows ];
-  yDerivativeMatrix<double>( nCols, nRows, xyDerivativeMatrixGreen, xDerivativeMatrixGreen );
-  double* xyDerivativeMatrixBlue = new double[ nCols * nRows ];
-  yDerivativeMatrix<double>( nCols, nRows, xyDerivativeMatrixBlue, xDerivativeMatrixBlue );
+  yDerivativeMatrix( nCols, nRows, yDerivativeMatrixBlue, blueMatrix );
 
   //compute output
   double nSrcPerDstX = ( double ) srcImage.width() / ( double ) dstImage.width();
@@ -81,47 +73,106 @@ void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage
   double currentSrcCol;
   int currentSrcColInt;
   int currentSrcRowInt;
+  int lastSrcColInt = -1;
+  int lastSrcRowInt = -1;
 
   int r, g, b;
-  int index;
+  //bernstein polynomials
+  double bp0u, bp1u, bp2u, bp3u, bp0v, bp1v, bp2v, bp3v;
 
   for ( int i = 0; i < dstImage.height(); ++i )
   {
-    currentSrcRowInt = (int) currentSrcRow;
+    currentSrcRowInt = ( int ) currentSrcRow;
 
     currentSrcCol = nSrcPerDstY / 2.0 - 0.5;
-    for( int j = 0; j < dstImage.width(); ++j )
+    for ( int j = 0; j < dstImage.width(); ++j )
     {
+      double u = currentSrcCol - ( int )currentSrcCol;
+      double v = currentSrcRow - ( int )currentSrcRow;
+
       //out of bounds check
-      currentSrcColInt = (int)currentSrcCol;
-      if( currentSrcColInt < 0 || currentSrcColInt > srcImage.width() - 2
-        || currentSrcRowInt < 0 || currentSrcRowInt > srcImage.height() - 2 )
+      currentSrcColInt = ( int )currentSrcCol;
+      if ( currentSrcColInt < 0 || currentSrcColInt > srcImage.width() - 2
+           || currentSrcRowInt < 0 || currentSrcRowInt > srcImage.height() - 2 )
       {
+        lastSrcColInt = currentSrcColInt;
         currentSrcCol += nSrcPerDstX;
-        ++index;
         dstImage.setPixel( j, i, qRgb( 0, 0, 255 ) );
         continue;
       }
 
+      //first update the control points if necessary
+      if ( currentSrcColInt != lastSrcColInt || currentSrcRowInt != lastSrcRowInt )
+      {
+        calculateControlPoints( nCols, nRows, currentSrcRowInt, currentSrcColInt, redMatrix, greenMatrix, blueMatrix,
+                                xDerivativeMatrixRed, xDerivativeMatrixGreen, xDerivativeMatrixBlue,
+                                yDerivativeMatrixRed, yDerivativeMatrixGreen, yDerivativeMatrixBlue );
+      }
+
+      //bernstein polynomials
+      bp0u = calcBernsteinPoly( 3, 0, u ); bp1u = calcBernsteinPoly( 3, 1, u );
+      bp2u = calcBernsteinPoly( 3, 2, u ); bp3u = calcBernsteinPoly( 3, 3, u );
+      bp0v = calcBernsteinPoly( 3, 0, v ); bp1v = calcBernsteinPoly( 3, 1, v );
+      bp2v = calcBernsteinPoly( 3, 2, v ); bp3v = calcBernsteinPoly( 3, 3, v );
 
-      //interpolation
-      index = currentSrcRowInt * nCols + currentSrcColInt;
-      r = cubicInterpolation( redMatrix[index], redMatrix[index + 1], redMatrix[index + nCols + 1 ], redMatrix[ index + nCols ],
-                              xDerivativeMatrixRed[index], xDerivativeMatrixRed[index + 1], xDerivativeMatrixRed[index + nCols + 1 ], xDerivativeMatrixRed[ index + nCols ],
-                              yDerivativeMatrixRed[index], yDerivativeMatrixRed[index + 1], yDerivativeMatrixRed[index + nCols + 1 ], yDerivativeMatrixRed[ index + nCols ],
-                              xyDerivativeMatrixRed[index], xyDerivativeMatrixRed[index + 1], xyDerivativeMatrixRed[index + nCols + 1 ], xyDerivativeMatrixRed[ index + nCols ] );
-      g = cubicInterpolation( greenMatrix[index], greenMatrix[index + 1], greenMatrix[index + nCols + 1 ], greenMatrix[ index + nCols ],
-                              xDerivativeMatrixGreen[index], xDerivativeMatrixGreen[index + 1], xDerivativeMatrixGreen[index + nCols + 1 ], xDerivativeMatrixGreen[ index + nCols ],
-                              yDerivativeMatrixGreen[index], yDerivativeMatrixGreen[index + 1], yDerivativeMatrixGreen[index + nCols + 1 ], yDerivativeMatrixGreen[ index + nCols ],
-                              xyDerivativeMatrixGreen[index], xyDerivativeMatrixGreen[index + 1], xyDerivativeMatrixGreen[index + nCols + 1 ], xyDerivativeMatrixGreen[ index + nCols ] );
-      b = cubicInterpolation( blueMatrix[index], blueMatrix[index + 1], blueMatrix[index + nCols + 1 ], blueMatrix[ index + nCols ],
-                              xDerivativeMatrixBlue[index], xDerivativeMatrixBlue[index + 1], xDerivativeMatrixBlue[index + nCols + 1 ], xDerivativeMatrixBlue[ index + nCols ],
-                              yDerivativeMatrixBlue[index], yDerivativeMatrixBlue[index + 1], yDerivativeMatrixBlue[index + nCols + 1 ], yDerivativeMatrixBlue[ index + nCols ],
-                              xyDerivativeMatrixBlue[index], xyDerivativeMatrixBlue[index + 1], xyDerivativeMatrixBlue[index + nCols + 1 ], xyDerivativeMatrixBlue[ index + nCols ] );
-
-      dstImage.setPixel( j, i, qRgb( 255, 0, 0 ) );
+      //then calculate value based on bernstein form of Bezier patch
+      //todo: move into function
+      r = bp0u * bp0v * cRed00 +
+          bp1u * bp0v * cRed10 +
+          bp2u * bp0v * cRed20 +
+          bp3u * bp0v * cRed30 +
+          bp0u * bp1v * cRed01 +
+          bp1u * bp1v * cRed11 +
+          bp2u * bp1v * cRed21 +
+          bp3u * bp1v * cRed31 +
+          bp0u * bp2v * cRed02 +
+          bp1u * bp2v * cRed12 +
+          bp2u * bp2v * cRed22 +
+          bp3u * bp2v * cRed32 +
+          bp0u * bp3v * cRed03 +
+          bp1u * bp3v * cRed13 +
+          bp2u * bp3v * cRed23 +
+          bp3u * bp3v * cRed33;
+
+      g = bp0u * bp0v * cGreen00 +
+          bp1u * bp0v * cGreen10 +
+          bp2u * bp0v * cGreen20 +
+          bp3u * bp0v * cGreen30 +
+          bp0u * bp1v * cGreen01 +
+          bp1u * bp1v * cGreen11 +
+          bp2u * bp1v * cGreen21 +
+          bp3u * bp1v * cGreen31 +
+          bp0u * bp2v * cGreen02 +
+          bp1u * bp2v * cGreen12 +
+          bp2u * bp2v * cGreen22 +
+          bp3u * bp2v * cGreen32 +
+          bp0u * bp3v * cGreen03 +
+          bp1u * bp3v * cGreen13 +
+          bp2u * bp3v * cGreen23 +
+          bp3u * bp3v * cGreen33;
+
+      b = bp0u * bp0v * cBlue00 +
+          bp1u * bp0v * cBlue10 +
+          bp2u * bp0v * cBlue20 +
+          bp3u * bp0v * cBlue30 +
+          bp0u * bp1v * cBlue01 +
+          bp1u * bp1v * cBlue11 +
+          bp2u * bp1v * cBlue21 +
+          bp3u * bp1v * cBlue31 +
+          bp0u * bp2v * cBlue02 +
+          bp1u * bp2v * cBlue12 +
+          bp2u * bp2v * cBlue22 +
+          bp3u * bp2v * cBlue32 +
+          bp0u * bp3v * cBlue03 +
+          bp1u * bp3v * cBlue13 +
+          bp2u * bp3v * cBlue23 +
+          bp3u * bp3v * cBlue33;
+
+      dstImage.setPixel( j, i, qRgb( r, g, b ) );
+      lastSrcColInt = currentSrcColInt;
       currentSrcCol += nSrcPerDstX;
     }
+    lastSrcRowInt = currentSrcRowInt;
     currentSrcRow += nSrcPerDstY;
   }
 
@@ -136,25 +187,22 @@ void QgsCubicRasterResampler::resample( const QImage& srcImage, QImage& dstImage
   delete[] yDerivativeMatrixRed;
   delete[] yDerivativeMatrixGreen;
   delete[] yDerivativeMatrixBlue;
-  delete[] xyDerivativeMatrixRed;
-  delete[] xyDerivativeMatrixGreen;
-  delete[] xyDerivativeMatrixBlue;
 }
 
 void QgsCubicRasterResampler::xDerivativeMatrix( int nCols, int nRows, double* matrix, const int* colorMatrix )
 {
   double val;
   int index = 0;
 
-  for( int i = 0; i < nRows; ++i )
+  for ( int i = 0; i < nRows; ++i )
   {
-    for( int j = 0; j < nCols; ++j )
+    for ( int j = 0; j < nCols; ++j )
     {
-      if( j < 1 )
+      if ( j < 1 )
       {
         val = colorMatrix[index + 1] - colorMatrix[index];
       }
-      else if( j == ( nCols - 1 ) )
+      else if ( j == ( nCols - 1 ) )
       {
         val = colorMatrix[index] - colorMatrix[ index - 1 ];
       }
@@ -168,25 +216,128 @@ void QgsCubicRasterResampler::xDerivativeMatrix( int nCols, int nRows, double* m
   }
 }
 
-int QgsCubicRasterResampler::cubicInterpolation( double p1, double p2, double p3, double p4, double p1x, double p2x, double p3x, double p4x,
-                                   double p1y, double p2y, double p3y, double p4y, double p1xy, double p2xy, double p3xy, double p4xy )
+void QgsCubicRasterResampler::yDerivativeMatrix( int nCols, int nRows, double* matrix, const int* colorMatrix )
+{
+  double val;
+  int index = 0;
+
+  for ( int i = 0; i < nRows; ++i )
+  {
+    for ( int j = 0; j < nCols; ++j )
+    {
+      if ( i == 0 )
+      {
+        val = colorMatrix[ index + nRows ] - colorMatrix[ index ];
+      }
+      else if ( i == nRows - 1 )
+      {
+        val = colorMatrix[ index ] - colorMatrix[ index - nRows ];
+      }
+      else
+      {
+        val = ( colorMatrix[ index + nRows ] - colorMatrix[ index - nRows ] ) / 2.0;
+      }
+      matrix[index] = val;
+      ++index;
+    }
+  }
+}
+
+void QgsCubicRasterResampler::calculateControlPoints( int nCols, int nRows, int currentRow, int currentCol, int* redMatrix, int* greenMatrix, int* blueMatrix,
+    double* xDerivativeMatrixRed, double* xDerivativeMatrixGreen, double* xDerivativeMatrixBlue,
+    double* yDerivativeMatrixRed, double* yDerivativeMatrixGreen, double* yDerivativeMatrixBlue )
+{
+  int idx00 = currentRow * nCols + currentCol;
+  int idx10 = idx00 + 1;
+  int idx01 = idx00 + nCols;
+  int idx11 = idx01 + 1;
+
+  //corner points
+  cRed00 = redMatrix[idx00]; cGreen00 = greenMatrix[idx00]; cBlue00 = blueMatrix[idx00];
+  cRed30 = redMatrix[idx10]; cGreen30 = greenMatrix[idx10]; cBlue30 = blueMatrix[idx10];
+  cRed03 = redMatrix[idx01]; cGreen03 = greenMatrix[idx01]; cBlue03 = blueMatrix[idx01];
+  cRed33 = redMatrix[idx11]; cGreen33 = greenMatrix[idx11]; cBlue33 = blueMatrix[idx11];
+
+  //control points near c00
+  cRed10 = cRed00 + 0.3 * xDerivativeMatrixRed[idx00]; cGreen10 = cGreen00 + 0.3 * xDerivativeMatrixGreen[idx00]; cBlue10 = cBlue00 + 0.3 * xDerivativeMatrixBlue[idx00];
+  cRed01 = cRed00 + 0.3 * yDerivativeMatrixRed[idx00]; cGreen01 = cGreen00 + 0.3 * yDerivativeMatrixGreen[idx00]; cBlue01 = cBlue00 + 0.3 * yDerivativeMatrixBlue[idx00];
+  cRed11 = cRed10 + 0.3 * yDerivativeMatrixRed[idx00]; cGreen11 = cGreen10 + 0.3 * yDerivativeMatrixGreen[idx00]; cBlue11 = cBlue10 + 0.3 * yDerivativeMatrixBlue[idx00];
+
+  //control points near c30
+  cRed20 = cRed30 - 0.3 * xDerivativeMatrixRed[idx10]; cGreen20 = cGreen30 - 0.3 * xDerivativeMatrixGreen[idx10]; cBlue20 = cBlue30 - 0.3 * xDerivativeMatrixBlue[idx10];
+  cRed31 = cRed30 + 0.3 * yDerivativeMatrixRed[idx10]; cGreen31 = cGreen30 + 0.3 * yDerivativeMatrixGreen[idx10]; cBlue31 = cBlue30 + 0.3 * yDerivativeMatrixBlue[idx10];
+  cRed21 = cRed20 + 0.3 * yDerivativeMatrixRed[idx10]; cGreen21 = cGreen20 + 0.3 * yDerivativeMatrixGreen[idx10]; cBlue21 = cBlue20 + 0.3 * yDerivativeMatrixBlue[idx10];
+
+  //control points near c03
+  cRed13 = cRed03 + 0.3 * xDerivativeMatrixRed[idx01]; cGreen13 = cGreen03 + 0.3 * xDerivativeMatrixGreen[idx01]; cBlue13 = cBlue03 + 0.3 * xDerivativeMatrixBlue[idx01];
+  cRed02 = cRed03 - 0.3 * yDerivativeMatrixRed[idx01]; cGreen02 = cGreen03 - 0.3 * yDerivativeMatrixGreen[idx01]; cBlue02 = cBlue03 - 0.3 * yDerivativeMatrixBlue[idx01];
+  cRed12 = cRed02 + 0.3 * xDerivativeMatrixRed[idx01]; cGreen12 = cGreen02 + 0.3 * xDerivativeMatrixGreen[idx01]; cBlue12 = cBlue02 + 0.3 * xDerivativeMatrixBlue[idx01];
+
+  //control points near c33
+  cRed23 = cRed33 - 0.3 * xDerivativeMatrixRed[idx11]; cGreen23 = cGreen33 - 0.3 * xDerivativeMatrixGreen[idx11]; cBlue23 = cBlue33 - 0.3 * xDerivativeMatrixBlue[idx11];
+  cRed32 = cRed33 - 0.3 * yDerivativeMatrixRed[idx11]; cGreen32 = cGreen33 - 0.3 * yDerivativeMatrixGreen[idx11]; cBlue32 = cBlue33 - 0.3 * yDerivativeMatrixBlue[idx11];
+  cRed22 = cRed32 - 0.3 * xDerivativeMatrixRed[idx11]; cGreen22 = cGreen32 - 0.3 * xDerivativeMatrixGreen[idx11]; cBlue22 = cBlue32 - 0.3 * xDerivativeMatrixBlue[idx11];
+}
+
+double QgsCubicRasterResampler::calcBernsteinPoly( int n, int i, double t )
+{
+  if ( i < 0 )
+  {
+    return 0;
+  }
+
+  return lower( n, i )*power( t, i )*power(( 1 - t ), ( n - i ) );
+}
+
+int QgsCubicRasterResampler::lower( int n, int i )
 {
-#if 0
-  //calculate 16 coefficients
-  double a00 = p1 - 3 * p4 + 2* p3 - 3 * p1y + 9 * p4y - 6 * p3y + 2 * p1xy - 6 * p4xy+ 4 * p3xy;
-  double a10 = 3 * p4 - 2 * p3 - 9* p4y + 6 * p3y + 6 * p4xy - 4 * p3xy;
-  double a20 = 3 * p1y - 9 * p4y + 6 * p3y - 2 * p1xy + 6 * p4xy - 4 * p3xy;
-  double a30 = 9 * p4y - 6 * p3y - 6 * p4xy + 4 * p3xy;
-  double a01 = p2 - 2 * p4 + p3 - 3 * p2y + 6 * p4y -3 * p3y + 2 * p2xy - 4 * p4xy + 2 * p3xy;
-  double a11 = - p4 + p3 + 3 * p4y + -3 * p3y - 2 * p4xy;
-  double a21 = 3 * p2y - 6 * p4y + 3 * p3y -2 * p2xy + 4 * p4xy - 2 * p3xy;
-  double a31 = - 3 * p4y + 3 * p3y + 2 * p4xy -2 * p3xy;
-  double a02 = p1x - 3 * p4x + 2 * p3x - 2 * p1y + 6 * p4y - 4 * p3y + p1xy - 3 * p4xy + 2 * p3xy;
-  double a12 = 3 * p4x - 2 * p3x - 6 * p4y + 4 * p3y + 3 * p4xy - 2 * p3xy;
-  double a22 = - p1y + 3 * p4y - 2 * p3y + p1xy - 3 * p4xy + 2 * p3xy;
-  double a32 = -3 * p4y + 2 * p3y + 3 * p4xy - 2 * p3xy;
-  double a03 = p2x - 2 * p4x + p3x - 2 * p4y + 4 * p4y - 2 * p3y + p2xy - 2 * p4xy + p3xy;
-  double a13 = - p4x + p3x + 2 * p4y - 2 * p3y - p4xy + p3xy;
-#endif //0
-  return 0;
+  if ( i >= 0 && i <= n )
+  {
+    return faculty( n ) / ( faculty( i )*faculty( n - i ) );
+  }
+  else
+  {
+    return 0;
+  }
+}
+
+double QgsCubicRasterResampler::power( double a, int b )
+{
+  if ( b == 0 )
+  {
+    return 1;
+  }
+  double tmp = a;
+  for ( int i = 2; i <= qAbs(( double )b ); i++ )
+  {
+
+    a *= tmp;
+  }
+  if ( b > 0 )
+  {
+    return a;
+  }
+  else
+  {
+    return ( 1.0 / a );
+  }
+}
+
+int QgsCubicRasterResampler::faculty( int n )
+{
+  if ( n < 0 )//Is faculty also defined for negative integers?
+  {
+    return 0;
+  }
+  int i;
+  int result = n;
+
+  if ( n == 0 || n == 1 )
+    {return 1;}//faculty of 0 is 1!
+
+  for ( i = n - 1; i >= 2; i-- )
+  {
+    result *= i;
+  }
+  return result;
 }