Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
qgsrasterprojector.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrasterprojector.cpp - Raster projector
3  --------------------------------------
4  Date : Jan 16, 2011
5  Copyright : (C) 2005 by Radim Blazek
6  email : radim dot blazek at gmail dot com
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 /* $Id: qgsrasterprojector.cpp 15005 2011-01-08 16:35:21Z rblazek $ */
18 
19 #include <cassert>
20 
21 #include "qgslogger.h"
22 #include "qgsrasterprojector.h"
23 #include "qgscoordinatetransform.h"
24 
28  QgsRectangle theDestExtent,
29  int theDestRows, int theDestCols,
30  double theMaxSrcXRes, double theMaxSrcYRes,
31  QgsRectangle theExtent )
32  : mSrcCRS( theSrcCRS )
33  , mDestCRS( theDestCRS )
34  , mCoordinateTransform( theDestCRS, theSrcCRS )
35  , mDestExtent( theDestExtent )
36  , mExtent( theExtent )
37  , mDestRows( theDestRows ), mDestCols( theDestCols )
38  , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
39 {
40  QgsDebugMsg( "Entered" );
41  QgsDebugMsg( "theDestExtent = " + theDestExtent.toString() );
42 
45 
46  // Calculate tolerance
47  // TODO: Think it over better
48  // Note: we are checking on matrix each even point, that means taht the real error
49  // in that moment is approximately half size
50  double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes;
51  mSqrTolerance = myDestRes * myDestRes;
52 
53  // Initialize the matrix by corners and middle points
54  mCPCols = mCPRows = 3;
55  for ( int i = 0; i < mCPRows; i++ )
56  {
57  QList<QgsPoint> myRow;
58  myRow.append( QgsPoint() );
59  myRow.append( QgsPoint() );
60  myRow.append( QgsPoint() );
61  mCPMatrix.insert( i, myRow );
62  }
63  for ( int i = 0; i < mCPRows; i++ )
64  {
65  calcRow( i );
66  }
67 
68  while ( true )
69  {
70  bool myColsOK = checkCols();
71  if ( !myColsOK )
72  {
73  insertRows();
74  }
75  bool myRowsOK = checkRows();
76  if ( !myRowsOK )
77  {
78  insertCols();
79  }
80  if ( myColsOK && myRowsOK )
81  {
82  QgsDebugMsg( "CP matrix within tolerance" );
83  mApproximate = true;
84  break;
85  }
86  // What is the maximum reasonable size of transformatio matrix?
87  // TODO: consider better when to break - ratio
88  if ( mCPRows * mCPCols > 0.0625 * mDestRows * mDestCols )
89  {
90  QgsDebugMsg( "Too large CP matrix" );
91  mApproximate = false;
92  break;
93  }
94 
95  }
96  QgsDebugMsg( QString( "CPMatrix size: mCPRows = %1 mCPCols = %2" ).arg( mCPRows ).arg( mCPCols ) );
97  mDestRowsPerMatrixRow = ( float )mDestRows / ( mCPRows - 1 );
98  mDestColsPerMatrixCol = ( float )mDestCols / ( mCPCols - 1 );
99 
100  //QgsDebugMsg( "CPMatrix:\n" + cpToString() );
101 
102  // Calculate source dimensions
103  calcSrcExtent();
104  calcSrcRowsCols();
107 
108  // init helper points
111  calcHelper( 0, pHelperTop );
113  mHelperTopRow = 0;
114 }
115 
117 {
118  //delete mCoordinateTransform;
119 }
120 
122 {
123  /* Run around the mCPMatrix and find source extent */
124  // Attention, source limits are not necessarily on destination edges, e.g.
125  // for destination EPSG:32661 Polar Stereographic and source EPSG:4326,
126  // the maximum y may be in the middle of destination extent
127  // TODO: How to find extent exactly and quickly?
128  // For now, we runt through all matrix
129  QgsPoint myPoint = mCPMatrix[0][0];
130  mSrcExtent = QgsRectangle( myPoint.x(), myPoint.y(), myPoint.x(), myPoint.y() );
131  for ( int i = 0; i < mCPRows; i++ )
132  {
133  for ( int j = 0; j < mCPCols ; j++ )
134  {
135  myPoint = mCPMatrix[i][j];
136  mSrcExtent.combineExtentWith( myPoint.x(), myPoint.y() );
137  }
138  }
139  // Expand a bit to avoid possible approx coords falling out because of representation error?
140 
141  // If mMaxSrcXRes, mMaxSrcYRes are defined (fixed src resolution)
142  // align extent to src resolution to avoid jumping reprojected pixels
143  // because of shifting resampled grid
144  // Important especially if we are over mMaxSrcXRes, mMaxSrcYRes limits
145 
146  QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
147 
148  if ( mMaxSrcXRes > 0 )
149  {
150  // with floor/ceil it should work correctly also for mSrcExtent.xMinimum() < mExtent.xMinimum()
151  double col = floor(( mSrcExtent.xMinimum() - mExtent.xMinimum() ) / mMaxSrcXRes );
152  double x = mExtent.xMinimum() + col * mMaxSrcXRes;
153  mSrcExtent.setXMinimum( x );
154 
155  col = ceil(( mSrcExtent.xMaximum() - mExtent.xMinimum() ) / mMaxSrcXRes );
156  x = mExtent.xMinimum() + col * mMaxSrcXRes;
157  mSrcExtent.setXMaximum( x );
158  }
159  if ( mMaxSrcYRes > 0 )
160  {
161  double row = floor(( mExtent.yMaximum() - mSrcExtent.yMaximum() ) / mMaxSrcYRes );
162  double y = mExtent.yMaximum() - row * mMaxSrcYRes;
163  mSrcExtent.setYMaximum( y );
164 
165  row = ceil(( mExtent.yMaximum() - mSrcExtent.yMinimum() ) / mMaxSrcYRes );
166  y = mExtent.yMaximum() - row * mMaxSrcYRes;
167  mSrcExtent.setYMinimum( y );
168  }
169 
170 
171  QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
172 }
173 
175 {
176  QString myString;
177  for ( int i = 0; i < mCPRows; i++ )
178  {
179  if ( i > 0 ) myString += "\n";
180  for ( int j = 0; j < mCPCols; j++ )
181  {
182  if ( j > 0 ) myString += " ";
183  QgsPoint myPoint = mCPMatrix[i][j];
184  myString += myPoint.toString();
185  }
186  }
187  return myString;
188 }
189 
191 {
192  // Wee need to calculate minimum cell size in the source
193  // TODO: Think it over better, what is the right source resolution?
194  // Taking distances between cell centers projected to source along source
195  // axis would result in very high resolution
196  // TODO: different resolution for rows and cols ?
197 
198  // For now, we take cell sizes projected to source but not to source axes
199  double myDestColsPerMatrixCell = mDestCols / mCPCols;
200  double myDestRowsPerMatrixCell = mDestRows / mCPRows;
201  QgsDebugMsg( QString( "myDestColsPerMatrixCell = %1 myDestRowsPerMatrixCell = %2" ).arg( myDestColsPerMatrixCell ).arg( myDestRowsPerMatrixCell ) );
202 
203  double myMinSize = DBL_MAX;
204 
205  for ( int i = 0; i < mCPRows - 1; i++ )
206  {
207  for ( int j = 0; j < mCPCols - 1; j++ )
208  {
209  QgsPoint myPointA = mCPMatrix[i][j];
210  QgsPoint myPointB = mCPMatrix[i][j+1];
211  QgsPoint myPointC = mCPMatrix[i+1][j];
212  double mySize = sqrt( myPointA.sqrDist( myPointB ) ) / myDestColsPerMatrixCell;
213  if ( mySize < myMinSize ) { myMinSize = mySize; }
214 
215  mySize = sqrt( myPointA.sqrDist( myPointC ) ) / myDestRowsPerMatrixCell;
216  if ( mySize < myMinSize ) { myMinSize = mySize; }
217  }
218  }
219 
220  // Make it a bit higher resolution
221  // TODO: find the best coefficient, attention, increasing resolution for WMS
222  // is changing WMS content
223  myMinSize *= 0.75;
224 
225  QgsDebugMsg( QString( "mMaxSrcXRes = %1 mMaxSrcYRes = %2" ).arg( mMaxSrcXRes ).arg( mMaxSrcYRes ) );
226  // mMaxSrcXRes, mMaxSrcYRes may be 0 - no limit (WMS)
227  double myMinXSize = mMaxSrcXRes > myMinSize ? mMaxSrcXRes : myMinSize;
228  double myMinYSize = mMaxSrcYRes > myMinSize ? mMaxSrcYRes : myMinSize;
229  QgsDebugMsg( QString( "myMinXSize = %1 myMinYSize = %2" ).arg( myMinXSize ).arg( myMinYSize ) );
230  QgsDebugMsg( QString( "mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ) );
231 
232  // we have to round to keep alignment set in calcSrcExtent
233  mSrcRows = ( int ) qRound( mSrcExtent.height() / myMinYSize );
234  mSrcCols = ( int ) qRound( mSrcExtent.width() / myMinXSize );
235 
236  QgsDebugMsg( QString( "mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ) );
237 }
238 
239 
240 inline void QgsRasterProjector::destPointOnCPMatrix( int theRow, int theCol, double *theX, double *theY )
241 {
242  *theX = mDestExtent.xMinimum() + theCol * mDestExtent.width() / ( mCPCols - 1 );
243  *theY = mDestExtent.yMaximum() - theRow * mDestExtent.height() / ( mCPRows - 1 );
244 }
245 
246 inline int QgsRasterProjector::matrixRow( int theDestRow )
247 {
248  return ( int )( floor(( theDestRow + 0.5 ) / mDestRowsPerMatrixRow ) );
249 }
250 inline int QgsRasterProjector::matrixCol( int theDestCol )
251 {
252  return ( int )( floor(( theDestCol + 0.5 ) / mDestColsPerMatrixCol ) );
253 }
254 
255 QgsPoint QgsRasterProjector::srcPoint( int theDestRow, int theCol )
256 {
257  return QgsPoint();
258 }
259 
260 void QgsRasterProjector::calcHelper( int theMatrixRow, QgsPoint *thePoints )
261 {
262  // TODO?: should we also precalc dest cell center coordinates for x and y?
263  for ( int myDestCol = 0; myDestCol < mDestCols; myDestCol++ )
264  {
265  double myDestX = mDestExtent.xMinimum() + ( myDestCol + 0.5 ) * mDestXRes;
266 
267  int myMatrixCol = matrixCol( myDestCol );
268 
269  double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
270 
271  destPointOnCPMatrix( theMatrixRow, myMatrixCol, &myDestXMin, &myDestYMin );
272  destPointOnCPMatrix( theMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
273 
274  double xfrac = ( myDestX - myDestXMin ) / ( myDestXMax - myDestXMin );
275 
276  QgsPoint &mySrcPoint0 = mCPMatrix[theMatrixRow][myMatrixCol];
277  QgsPoint &mySrcPoint1 = mCPMatrix[theMatrixRow][myMatrixCol+1];
278  double s = mySrcPoint0.x() + ( mySrcPoint1.x() - mySrcPoint0.x() ) * xfrac;
279  double t = mySrcPoint0.y() + ( mySrcPoint1.y() - mySrcPoint0.y() ) * xfrac;
280 
281  thePoints[myDestCol].setX( s );
282  thePoints[myDestCol].setY( t );
283  }
284 }
286 {
287  QgsPoint *tmp;
288  tmp = pHelperTop;
290  pHelperBottom = tmp;
292  mHelperTopRow++;
293 }
294 
295 void QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
296 {
297  if ( mApproximate ) approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
298  else preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
299 }
300 
301 void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
302 {
303  //QgsDebugMsg( QString( "theDestRow = %1" ).arg(theDestRow) );
304  //QgsDebugMsg( QString( "theDestRow = %1 mDestExtent.yMaximum() = %2 mDestYRes = %3" ).arg(theDestRow).arg(mDestExtent.yMaximum()).arg(mDestYRes) );
305 
306  // Get coordinate of center of destination cell
307  double x = mDestExtent.xMinimum() + ( theDestCol + 0.5 ) * mDestXRes;
308  double y = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
309  double z = 0;
310 
311  //QgsDebugMsg( QString( "x = %1 y = %2" ).arg(x).arg(y) );
313  //QgsDebugMsg( QString( "x = %1 y = %2" ).arg(x).arg(y) );
314 
315  // Get source row col
316  *theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - y ) / mSrcYRes );
317  *theSrcCol = ( int ) floor(( x - mSrcExtent.xMinimum() ) / mSrcXRes );
318  //QgsDebugMsg( QString( "mSrcExtent.yMaximum() = %1 mSrcYRes = %2" ).arg(mSrcExtent.yMaximum()).arg(mSrcYRes) );
319  //QgsDebugMsg( QString( "theSrcRow = %1 theSrcCol = %2" ).arg(*theSrcRow).arg(*theSrcCol) );
320 
321  // With epsg 32661 (Polar Stereographic) it was happening that *theSrcCol == mSrcCols
322  // For now silently correct limits to avoid crashes
323  // TODO: review
324  if ( *theSrcRow >= mSrcRows ) *theSrcRow = mSrcRows - 1;
325  if ( *theSrcRow < 0 ) *theSrcRow = 0;
326  if ( *theSrcCol >= mSrcCols ) *theSrcCol = mSrcCols - 1;
327  if ( *theSrcCol < 0 ) *theSrcCol = 0;
328 
329  assert( *theSrcRow < mSrcRows );
330  assert( *theSrcCol < mSrcCols );
331 }
332 
333 void QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
334 {
335  int myMatrixRow = matrixRow( theDestRow );
336  int myMatrixCol = matrixCol( theDestCol );
337 
338  if ( myMatrixRow > mHelperTopRow )
339  {
340  // TODO: make it more robust (for random, not sequential reading)
341  nextHelper();
342  }
343 
344  double myDestY = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
345 
346  // See the schema in javax.media.jai.WarpGrid doc (but up side down)
347  // TODO: use some kind of cache of values which can be reused
348  double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
349 
350  destPointOnCPMatrix( myMatrixRow + 1, myMatrixCol, &myDestXMin, &myDestYMin );
351  destPointOnCPMatrix( myMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
352 
353  double yfrac = ( myDestY - myDestYMin ) / ( myDestYMax - myDestYMin );
354 
355  QgsPoint &myTop = pHelperTop[theDestCol];
356  QgsPoint &myBot = pHelperBottom[theDestCol];
357 
358  // Warning: this is very SLOW compared to the following code!:
359  //double mySrcX = myBot.x() + (myTop.x() - myBot.x()) * yfrac;
360  //double mySrcY = myBot.y() + (myTop.y() - myBot.y()) * yfrac;
361 
362  double tx = myTop.x();
363  double ty = myTop.y();
364  double bx = myBot.x();
365  double by = myBot.y();
366  double mySrcX = bx + ( tx - bx ) * yfrac;
367  double mySrcY = by + ( ty - by ) * yfrac;
368 
369  // TODO: check again cell selection (coor is in the middle)
370 
371  *theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - mySrcY ) / mSrcYRes );
372  *theSrcCol = ( int ) floor(( mySrcX - mSrcExtent.xMinimum() ) / mSrcXRes );
373 
374  // For now silently correct limits to avoid crashes
375  // TODO: review
376  if ( *theSrcRow >= mSrcRows ) *theSrcRow = mSrcRows - 1;
377  if ( *theSrcRow < 0 ) *theSrcRow = 0;
378  if ( *theSrcCol >= mSrcCols ) *theSrcCol = mSrcCols - 1;
379  if ( *theSrcCol < 0 ) *theSrcCol = 0;
380  assert( *theSrcRow < mSrcRows );
381  assert( *theSrcCol < mSrcCols );
382 }
383 
385 {
386  for ( int r = 0; r < mCPRows - 1; r++ )
387  {
388  QList<QgsPoint> myRow;
389  for ( int c = 0; c < mCPCols; c++ )
390  {
391  myRow.append( QgsPoint() );
392  }
393  QgsDebugMsgLevel( QString( "insert new row at %1" ).arg( 1 + r*2 ), 3 );
394  mCPMatrix.insert( 1 + r*2, myRow );
395  }
396  mCPRows += mCPRows - 1;
397  for ( int r = 1; r < mCPRows - 1; r += 2 )
398  {
399  calcRow( r );
400  }
401 }
402 
404 {
405  for ( int r = 0; r < mCPRows; r++ )
406  {
407  QList<QgsPoint> myRow;
408  for ( int c = 0; c < mCPCols - 1; c++ )
409  {
410  mCPMatrix[r].insert( 1 + c*2, QgsPoint() );
411  }
412  }
413  mCPCols += mCPCols - 1;
414  for ( int c = 1; c < mCPCols - 1; c += 2 )
415  {
416  calcCol( c );
417  }
418 
419 }
420 
421 void QgsRasterProjector::calcCP( int theRow, int theCol )
422 {
423  double myDestX, myDestY;
424  destPointOnCPMatrix( theRow, theCol, &myDestX, &myDestY );
425  QgsPoint myDestPoint( myDestX, myDestY );
426 
427  mCPMatrix[theRow][theCol] = mCoordinateTransform.transform( myDestPoint );
428 }
429 
430 bool QgsRasterProjector::calcRow( int theRow )
431 {
432  QgsDebugMsgLevel( QString( "theRow = %1" ).arg( theRow ), 3 );
433  for ( int i = 0; i < mCPCols; i++ )
434  {
435  calcCP( theRow, i );
436  }
437 
438  return true;
439 }
440 
441 bool QgsRasterProjector::calcCol( int theCol )
442 {
443  QgsDebugMsgLevel( QString( "theCol = %1" ).arg( theCol ), 3 );
444  for ( int i = 0; i < mCPRows; i++ )
445  {
446  calcCP( i, theCol );
447  }
448 
449  return true;
450 }
451 
453 {
454  for ( int c = 0; c < mCPCols; c++ )
455  {
456  for ( int r = 1; r < mCPRows - 1; r += 2 )
457  {
458  double myDestX, myDestY;
459  destPointOnCPMatrix( r, c, &myDestX, &myDestY );
460  QgsPoint myDestPoint( myDestX, myDestY );
461 
462  QgsPoint mySrcPoint1 = mCPMatrix[r-1][c];
463  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
464  QgsPoint mySrcPoint3 = mCPMatrix[r+1][c];
465 
466  QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
468  double mySqrDist = myDestApprox.sqrDist( myDestPoint );
469  if ( mySqrDist > mSqrTolerance ) { return false; }
470  }
471  }
472  return true;
473 }
474 
476 {
477  for ( int r = 0; r < mCPRows; r++ )
478  {
479  for ( int c = 1; c < mCPCols - 1; c += 2 )
480  {
481  double myDestX, myDestY;
482  destPointOnCPMatrix( r, c, &myDestX, &myDestY );
483 
484  QgsPoint myDestPoint( myDestX, myDestY );
485  QgsPoint mySrcPoint1 = mCPMatrix[r][c-1];
486  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
487  QgsPoint mySrcPoint3 = mCPMatrix[r][c+1];
488 
489  QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
491  double mySqrDist = myDestApprox.sqrDist( myDestPoint );
492  if ( mySqrDist > mSqrTolerance ) { return false; }
493  }
494  }
495  return true;
496 }