Quantum GIS API Documentation  1.8
src/core/qgsrasterprojector.cpp
Go to the documentation of this file.
00001 /***************************************************************************
00002     qgsrasterprojector.cpp - Raster projector
00003      --------------------------------------
00004     Date                 : Jan 16, 2011
00005     Copyright            : (C) 2005 by Radim Blazek
00006     email                : radim dot blazek at gmail dot com
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 "qgslogger.h"
00019 #include "qgsrasterprojector.h"
00020 #include "qgscoordinatetransform.h"
00021 
00022 QgsRasterProjector::QgsRasterProjector(
00023   QgsCoordinateReferenceSystem theSrcCRS,
00024   QgsCoordinateReferenceSystem theDestCRS,
00025   QgsRectangle theDestExtent,
00026   int theDestRows, int theDestCols,
00027   double theMaxSrcXRes, double theMaxSrcYRes,
00028   QgsRectangle theExtent )
00029     : mSrcCRS( theSrcCRS )
00030     , mDestCRS( theDestCRS )
00031     , mCoordinateTransform( theDestCRS, theSrcCRS )
00032     , mDestExtent( theDestExtent )
00033     , mExtent( theExtent )
00034     , mDestRows( theDestRows ), mDestCols( theDestCols )
00035     , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
00036 {
00037   QgsDebugMsg( "Entered" );
00038   QgsDebugMsg( "theDestExtent = " + theDestExtent.toString() );
00039 
00040   mDestXRes = mDestExtent.width() / ( mDestCols );
00041   mDestYRes = mDestExtent.height() / ( mDestRows );
00042 
00043   // Calculate tolerance
00044   // TODO: Think it over better
00045   // Note: we are checking on matrix each even point, that means that the real error
00046   // in that moment is approximately half size
00047   double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes;
00048   mSqrTolerance = myDestRes * myDestRes;
00049 
00050   // Initialize the matrix by corners and middle points
00051   mCPCols = mCPRows = 3;
00052   for ( int i = 0; i < mCPRows; i++ )
00053   {
00054     QList<QgsPoint> myRow;
00055     myRow.append( QgsPoint() );
00056     myRow.append( QgsPoint() );
00057     myRow.append( QgsPoint() );
00058     mCPMatrix.insert( i,  myRow );
00059   }
00060   for ( int i = 0; i < mCPRows; i++ )
00061   {
00062     calcRow( i );
00063   }
00064 
00065   while ( true )
00066   {
00067     bool myColsOK = checkCols();
00068     if ( !myColsOK )
00069     {
00070       insertRows();
00071     }
00072     bool myRowsOK = checkRows();
00073     if ( !myRowsOK )
00074     {
00075       insertCols();
00076     }
00077     if ( myColsOK && myRowsOK )
00078     {
00079       QgsDebugMsg( "CP matrix within tolerance" );
00080       mApproximate = true;
00081       break;
00082     }
00083     // What is the maximum reasonable size of transformatio matrix?
00084     // TODO: consider better when to break - ratio
00085     if ( mCPRows * mCPCols > 0.0625 * mDestRows * mDestCols )
00086     {
00087       QgsDebugMsg( "Too large CP matrix" );
00088       mApproximate = false;
00089       break;
00090     }
00091 
00092   }
00093   QgsDebugMsg( QString( "CPMatrix size: mCPRows = %1 mCPCols = %2" ).arg( mCPRows ).arg( mCPCols ) );
00094   mDestRowsPerMatrixRow = ( float )mDestRows / ( mCPRows - 1 );
00095   mDestColsPerMatrixCol = ( float )mDestCols / ( mCPCols - 1 );
00096 
00097   //QgsDebugMsg( "CPMatrix:\n" + cpToString() );
00098 
00099   // Calculate source dimensions
00100   calcSrcExtent();
00101   calcSrcRowsCols();
00102   mSrcYRes = mSrcExtent.height() / mSrcRows;
00103   mSrcXRes = mSrcExtent.width() / mSrcCols;
00104 
00105   // init helper points
00106   pHelperTop = new QgsPoint[mDestCols];
00107   pHelperBottom = new QgsPoint[mDestCols];
00108   calcHelper( 0, pHelperTop );
00109   calcHelper( 1, pHelperBottom );
00110   mHelperTopRow = 0;
00111 }
00112 
00113 QgsRasterProjector::~QgsRasterProjector()
00114 {
00115   //delete mCoordinateTransform;
00116 }
00117 
00118 void QgsRasterProjector::calcSrcExtent()
00119 {
00120   /* Run around the mCPMatrix and find source extent */
00121   // Attention, source limits are not necessarily on destination edges, e.g.
00122   // for destination EPSG:32661 Polar Stereographic and source EPSG:4326,
00123   // the maximum y may be in the middle of destination extent
00124   // TODO: How to find extent exactly and quickly?
00125   // For now, we runt through all matrix
00126   QgsPoint myPoint = mCPMatrix[0][0];
00127   mSrcExtent = QgsRectangle( myPoint.x(), myPoint.y(), myPoint.x(), myPoint.y() );
00128   for ( int i = 0; i < mCPRows; i++ )
00129   {
00130     for ( int j = 0; j < mCPCols ; j++ )
00131     {
00132       myPoint = mCPMatrix[i][j];
00133       mSrcExtent.combineExtentWith( myPoint.x(), myPoint.y() );
00134     }
00135   }
00136   // Expand a bit to avoid possible approx coords falling out because of representation error?
00137 
00138   // If mMaxSrcXRes, mMaxSrcYRes are defined (fixed src resolution)
00139   // align extent to src resolution to avoid jumping reprojected pixels
00140   // because of shifting resampled grid
00141   // Important especially if we are over mMaxSrcXRes, mMaxSrcYRes limits
00142 
00143   QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
00144 
00145   if ( mMaxSrcXRes > 0 )
00146   {
00147     // with floor/ceil it should work correctly also for mSrcExtent.xMinimum() < mExtent.xMinimum()
00148     double col = floor(( mSrcExtent.xMinimum() - mExtent.xMinimum() ) / mMaxSrcXRes );
00149     double x = mExtent.xMinimum() + col * mMaxSrcXRes;
00150     mSrcExtent.setXMinimum( x );
00151 
00152     col = ceil(( mSrcExtent.xMaximum() - mExtent.xMinimum() ) / mMaxSrcXRes );
00153     x = mExtent.xMinimum() + col * mMaxSrcXRes;
00154     mSrcExtent.setXMaximum( x );
00155   }
00156   if ( mMaxSrcYRes > 0 )
00157   {
00158     double row = floor(( mExtent.yMaximum() - mSrcExtent.yMaximum() ) / mMaxSrcYRes );
00159     double y = mExtent.yMaximum() - row * mMaxSrcYRes;
00160     mSrcExtent.setYMaximum( y );
00161 
00162     row = ceil(( mExtent.yMaximum() - mSrcExtent.yMinimum() ) / mMaxSrcYRes );
00163     y = mExtent.yMaximum() - row * mMaxSrcYRes;
00164     mSrcExtent.setYMinimum( y );
00165   }
00166 
00167 
00168   QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
00169 }
00170 
00171 QString QgsRasterProjector::cpToString()
00172 {
00173   QString myString;
00174   for ( int i = 0; i < mCPRows; i++ )
00175   {
00176     if ( i > 0 )
00177       myString += "\n";
00178     for ( int j = 0; j < mCPCols; j++ )
00179     {
00180       if ( j > 0 )
00181         myString += "  ";
00182       QgsPoint myPoint = mCPMatrix[i][j];
00183       myString += myPoint.toString();
00184     }
00185   }
00186   return myString;
00187 }
00188 
00189 void QgsRasterProjector::calcSrcRowsCols()
00190 {
00191   // Wee need to calculate minimum cell size in the source
00192   // TODO: Think it over better, what is the right source resolution?
00193   //       Taking distances between cell centers projected to source along source
00194   //       axis would result in very high resolution
00195   // TODO: different resolution for rows and cols ?
00196 
00197   // For now, we take cell sizes projected to source but not to source axes
00198   double myDestColsPerMatrixCell = mDestCols / mCPCols;
00199   double myDestRowsPerMatrixCell = mDestRows / mCPRows;
00200   QgsDebugMsg( QString( "myDestColsPerMatrixCell = %1 myDestRowsPerMatrixCell = %2" ).arg( myDestColsPerMatrixCell ).arg( myDestRowsPerMatrixCell ) );
00201 
00202   double myMinSize = DBL_MAX;
00203 
00204   for ( int i = 0; i < mCPRows - 1; i++ )
00205   {
00206     for ( int j = 0; j < mCPCols - 1; j++ )
00207     {
00208       QgsPoint myPointA = mCPMatrix[i][j];
00209       QgsPoint myPointB = mCPMatrix[i][j+1];
00210       QgsPoint myPointC = mCPMatrix[i+1][j];
00211       double mySize = sqrt( myPointA.sqrDist( myPointB ) ) / myDestColsPerMatrixCell;
00212       if ( mySize < myMinSize )
00213         myMinSize = mySize;
00214 
00215       mySize = sqrt( myPointA.sqrDist( myPointC ) ) / myDestRowsPerMatrixCell;
00216       if ( mySize < myMinSize )
00217         myMinSize = mySize;
00218     }
00219   }
00220 
00221   // Make it a bit higher resolution
00222   // TODO: find the best coefficient, attention, increasing resolution for WMS
00223   // is changing WMS content
00224   myMinSize *= 0.75;
00225 
00226   QgsDebugMsg( QString( "mMaxSrcXRes = %1 mMaxSrcYRes = %2" ).arg( mMaxSrcXRes ).arg( mMaxSrcYRes ) );
00227   // mMaxSrcXRes, mMaxSrcYRes may be 0 - no limit (WMS)
00228   double myMinXSize = mMaxSrcXRes > myMinSize ? mMaxSrcXRes : myMinSize;
00229   double myMinYSize = mMaxSrcYRes > myMinSize ? mMaxSrcYRes : myMinSize;
00230   QgsDebugMsg( QString( "myMinXSize = %1 myMinYSize = %2" ).arg( myMinXSize ).arg( myMinYSize ) );
00231   QgsDebugMsg( QString( "mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ) );
00232 
00233   // we have to round to keep alignment set in calcSrcExtent
00234   mSrcRows = ( int ) qRound( mSrcExtent.height() / myMinYSize );
00235   mSrcCols = ( int ) qRound( mSrcExtent.width() / myMinXSize );
00236 
00237   QgsDebugMsg( QString( "mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ) );
00238 }
00239 
00240 
00241 inline void QgsRasterProjector::destPointOnCPMatrix( int theRow, int theCol, double *theX, double *theY )
00242 {
00243   *theX = mDestExtent.xMinimum() + theCol * mDestExtent.width() / ( mCPCols - 1 );
00244   *theY = mDestExtent.yMaximum() - theRow * mDestExtent.height() / ( mCPRows - 1 );
00245 }
00246 
00247 inline int QgsRasterProjector::matrixRow( int theDestRow )
00248 {
00249   return ( int )( floor(( theDestRow + 0.5 ) / mDestRowsPerMatrixRow ) );
00250 }
00251 inline int QgsRasterProjector::matrixCol( int theDestCol )
00252 {
00253   return ( int )( floor(( theDestCol + 0.5 ) / mDestColsPerMatrixCol ) );
00254 }
00255 
00256 QgsPoint QgsRasterProjector::srcPoint( int theDestRow, int theCol )
00257 {
00258   Q_UNUSED( theDestRow );
00259   Q_UNUSED( theCol );
00260   return QgsPoint();
00261 }
00262 
00263 void QgsRasterProjector::calcHelper( int theMatrixRow, QgsPoint *thePoints )
00264 {
00265   // TODO?: should we also precalc dest cell center coordinates for x and y?
00266   for ( int myDestCol = 0; myDestCol < mDestCols; myDestCol++ )
00267   {
00268     double myDestX = mDestExtent.xMinimum() + ( myDestCol + 0.5 ) * mDestXRes;
00269 
00270     int myMatrixCol = matrixCol( myDestCol );
00271 
00272     double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
00273 
00274     destPointOnCPMatrix( theMatrixRow, myMatrixCol, &myDestXMin, &myDestYMin );
00275     destPointOnCPMatrix( theMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
00276 
00277     double xfrac = ( myDestX - myDestXMin ) / ( myDestXMax - myDestXMin );
00278 
00279     QgsPoint &mySrcPoint0 = mCPMatrix[theMatrixRow][myMatrixCol];
00280     QgsPoint &mySrcPoint1 = mCPMatrix[theMatrixRow][myMatrixCol+1];
00281     double s = mySrcPoint0.x() + ( mySrcPoint1.x() - mySrcPoint0.x() ) * xfrac;
00282     double t = mySrcPoint0.y() + ( mySrcPoint1.y() - mySrcPoint0.y() ) * xfrac;
00283 
00284     thePoints[myDestCol].setX( s );
00285     thePoints[myDestCol].setY( t );
00286   }
00287 }
00288 void QgsRasterProjector::nextHelper()
00289 {
00290   QgsPoint *tmp;
00291   tmp = pHelperTop;
00292   pHelperTop = pHelperBottom;
00293   pHelperBottom = tmp;
00294   calcHelper( mHelperTopRow + 2, pHelperBottom );
00295   mHelperTopRow++;
00296 }
00297 
00298 void QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
00299 {
00300   if ( mApproximate )
00301     approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
00302   else preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
00303 }
00304 
00305 void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
00306 {
00307   //QgsDebugMsg( QString( "theDestRow = %1" ).arg(theDestRow) );
00308   //QgsDebugMsg( QString( "theDestRow = %1 mDestExtent.yMaximum() = %2 mDestYRes = %3" ).arg(theDestRow).arg(mDestExtent.yMaximum()).arg(mDestYRes) );
00309 
00310   // Get coordinate of center of destination cell
00311   double x = mDestExtent.xMinimum() + ( theDestCol + 0.5 ) * mDestXRes;
00312   double y = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
00313   double z = 0;
00314 
00315   //QgsDebugMsg( QString( "x = %1 y = %2" ).arg(x).arg(y) );
00316   mCoordinateTransform.transformInPlace( x, y, z );
00317   //QgsDebugMsg( QString( "x = %1 y = %2" ).arg(x).arg(y) );
00318 
00319   // Get source row col
00320   *theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - y ) / mSrcYRes );
00321   *theSrcCol = ( int ) floor(( x - mSrcExtent.xMinimum() ) / mSrcXRes );
00322   //QgsDebugMsg( QString( "mSrcExtent.yMaximum() = %1 mSrcYRes = %2" ).arg(mSrcExtent.yMaximum()).arg(mSrcYRes) );
00323   //QgsDebugMsg( QString( "theSrcRow = %1 theSrcCol = %2" ).arg(*theSrcRow).arg(*theSrcCol) );
00324 
00325   // With epsg 32661 (Polar Stereographic) it was happening that *theSrcCol == mSrcCols
00326   // For now silently correct limits to avoid crashes
00327   // TODO: review
00328   if ( *theSrcRow >= mSrcRows )
00329     *theSrcRow = mSrcRows - 1;
00330   if ( *theSrcRow < 0 )
00331     *theSrcRow = 0;
00332   if ( *theSrcCol >= mSrcCols )
00333     *theSrcCol = mSrcCols - 1;
00334   if ( *theSrcCol < 0 )
00335     *theSrcCol = 0;
00336 
00337   Q_ASSERT( *theSrcRow < mSrcRows );
00338   Q_ASSERT( *theSrcCol < mSrcCols );
00339 }
00340 
00341 void QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
00342 {
00343   int myMatrixRow = matrixRow( theDestRow );
00344   int myMatrixCol = matrixCol( theDestCol );
00345 
00346   if ( myMatrixRow > mHelperTopRow )
00347   {
00348     // TODO: make it more robust (for random, not sequential reading)
00349     nextHelper();
00350   }
00351 
00352   double myDestY = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
00353 
00354   // See the schema in javax.media.jai.WarpGrid doc (but up side down)
00355   // TODO: use some kind of cache of values which can be reused
00356   double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
00357 
00358   destPointOnCPMatrix( myMatrixRow + 1, myMatrixCol, &myDestXMin, &myDestYMin );
00359   destPointOnCPMatrix( myMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
00360 
00361   double yfrac = ( myDestY - myDestYMin ) / ( myDestYMax - myDestYMin );
00362 
00363   QgsPoint &myTop = pHelperTop[theDestCol];
00364   QgsPoint &myBot = pHelperBottom[theDestCol];
00365 
00366   // Warning: this is very SLOW compared to the following code!:
00367   //double mySrcX = myBot.x() + (myTop.x() - myBot.x()) * yfrac;
00368   //double mySrcY = myBot.y() + (myTop.y() - myBot.y()) * yfrac;
00369 
00370   double tx = myTop.x();
00371   double ty = myTop.y();
00372   double bx = myBot.x();
00373   double by = myBot.y();
00374   double mySrcX = bx + ( tx - bx ) * yfrac;
00375   double mySrcY = by + ( ty - by ) * yfrac;
00376 
00377   // TODO: check again cell selection (coor is in the middle)
00378 
00379   *theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - mySrcY ) / mSrcYRes );
00380   *theSrcCol = ( int ) floor(( mySrcX - mSrcExtent.xMinimum() ) / mSrcXRes );
00381 
00382   // For now silently correct limits to avoid crashes
00383   // TODO: review
00384   if ( *theSrcRow >= mSrcRows )
00385     *theSrcRow = mSrcRows - 1;
00386   if ( *theSrcRow < 0 )
00387     *theSrcRow = 0;
00388   if ( *theSrcCol >= mSrcCols )
00389     *theSrcCol = mSrcCols - 1;
00390   if ( *theSrcCol < 0 )
00391     *theSrcCol = 0;
00392   Q_ASSERT( *theSrcRow < mSrcRows );
00393   Q_ASSERT( *theSrcCol < mSrcCols );
00394 }
00395 
00396 void QgsRasterProjector::insertRows()
00397 {
00398   for ( int r = 0; r < mCPRows - 1; r++ )
00399   {
00400     QList<QgsPoint> myRow;
00401     for ( int c = 0; c < mCPCols; c++ )
00402     {
00403       myRow.append( QgsPoint() );
00404     }
00405     QgsDebugMsgLevel( QString( "insert new row at %1" ).arg( 1 + r*2 ), 3 );
00406     mCPMatrix.insert( 1 + r*2,  myRow );
00407   }
00408   mCPRows += mCPRows - 1;
00409   for ( int r = 1; r < mCPRows - 1; r += 2 )
00410   {
00411     calcRow( r );
00412   }
00413 }
00414 
00415 void QgsRasterProjector::insertCols()
00416 {
00417   for ( int r = 0; r < mCPRows; r++ )
00418   {
00419     QList<QgsPoint> myRow;
00420     for ( int c = 0; c < mCPCols - 1; c++ )
00421     {
00422       mCPMatrix[r].insert( 1 + c*2,  QgsPoint() );
00423     }
00424   }
00425   mCPCols += mCPCols - 1;
00426   for ( int c = 1; c < mCPCols - 1; c += 2 )
00427   {
00428     calcCol( c );
00429   }
00430 
00431 }
00432 
00433 void QgsRasterProjector::calcCP( int theRow, int theCol )
00434 {
00435   double myDestX, myDestY;
00436   destPointOnCPMatrix( theRow, theCol, &myDestX, &myDestY );
00437   QgsPoint myDestPoint( myDestX, myDestY );
00438 
00439   mCPMatrix[theRow][theCol] = mCoordinateTransform.transform( myDestPoint );
00440 }
00441 
00442 bool QgsRasterProjector::calcRow( int theRow )
00443 {
00444   QgsDebugMsgLevel( QString( "theRow = %1" ).arg( theRow ), 3 );
00445   for ( int i = 0; i < mCPCols; i++ )
00446   {
00447     calcCP( theRow, i );
00448   }
00449 
00450   return true;
00451 }
00452 
00453 bool QgsRasterProjector::calcCol( int theCol )
00454 {
00455   QgsDebugMsgLevel( QString( "theCol = %1" ).arg( theCol ), 3 );
00456   for ( int i = 0; i < mCPRows; i++ )
00457   {
00458     calcCP( i, theCol );
00459   }
00460 
00461   return true;
00462 }
00463 
00464 bool QgsRasterProjector::checkCols()
00465 {
00466   for ( int c = 0; c < mCPCols; c++ )
00467   {
00468     for ( int r = 1; r < mCPRows - 1; r += 2 )
00469     {
00470       double myDestX, myDestY;
00471       destPointOnCPMatrix( r, c, &myDestX, &myDestY );
00472       QgsPoint myDestPoint( myDestX, myDestY );
00473 
00474       QgsPoint mySrcPoint1 = mCPMatrix[r-1][c];
00475       QgsPoint mySrcPoint2 = mCPMatrix[r][c];
00476       QgsPoint mySrcPoint3 = mCPMatrix[r+1][c];
00477 
00478       QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
00479       QgsPoint myDestApprox = mCoordinateTransform.transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
00480       double mySqrDist = myDestApprox.sqrDist( myDestPoint );
00481       if ( mySqrDist > mSqrTolerance )
00482         return false;
00483     }
00484   }
00485   return true;
00486 }
00487 
00488 bool QgsRasterProjector::checkRows()
00489 {
00490   for ( int r = 0; r < mCPRows; r++ )
00491   {
00492     for ( int c = 1; c < mCPCols - 1; c += 2 )
00493     {
00494       double myDestX, myDestY;
00495       destPointOnCPMatrix( r, c, &myDestX, &myDestY );
00496 
00497       QgsPoint myDestPoint( myDestX, myDestY );
00498       QgsPoint mySrcPoint1 = mCPMatrix[r][c-1];
00499       QgsPoint mySrcPoint2 = mCPMatrix[r][c];
00500       QgsPoint mySrcPoint3 = mCPMatrix[r][c+1];
00501 
00502       QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
00503       QgsPoint myDestApprox = mCoordinateTransform.transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
00504       double mySqrDist = myDestApprox.sqrDist( myDestPoint );
00505       if ( mySqrDist > mSqrTolerance )
00506         return false;
00507     }
00508   }
00509   return true;
00510 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines