Quantum GIS API Documentation
1.8
|
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 }