QGIS API Documentation  2.0.1-Dufour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 
18 #include "qgsrasterdataprovider.h"
19 #include "qgslogger.h"
20 #include "qgsrasterprojector.h"
21 #include "qgscoordinatetransform.h"
22 
26  QgsRectangle theDestExtent,
27  int theDestRows, int theDestCols,
28  double theMaxSrcXRes, double theMaxSrcYRes,
29  QgsRectangle theExtent )
30  : QgsRasterInterface( 0 )
31  , mSrcCRS( theSrcCRS )
32  , mDestCRS( theDestCRS )
33  , mCoordinateTransform( theDestCRS, theSrcCRS )
34  , mDestExtent( theDestExtent )
35  , mExtent( theExtent )
36  , mDestRows( theDestRows ), mDestCols( theDestCols )
37  , pHelperTop( 0 ), pHelperBottom( 0 )
38  , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
39 {
40  QgsDebugMsg( "Entered" );
41  QgsDebugMsg( "theDestExtent = " + theDestExtent.toString() );
42 
43  calc();
44 }
45 
49  double theMaxSrcXRes, double theMaxSrcYRes,
50  QgsRectangle theExtent )
51  : QgsRasterInterface( 0 )
52  , mSrcCRS( theSrcCRS )
53  , mDestCRS( theDestCRS )
54  , mCoordinateTransform( theDestCRS, theSrcCRS )
55  , mExtent( theExtent )
56  , pHelperTop( 0 ), pHelperBottom( 0 )
57  , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
58 {
59  QgsDebugMsg( "Entered" );
60 }
61 
63  : QgsRasterInterface( 0 )
64  , pHelperTop( 0 ), pHelperBottom( 0 )
65 {
66  QgsDebugMsg( "Entered" );
67 }
68 
70  : QgsRasterInterface( 0 )
71 {
72  mSrcCRS = projector.mSrcCRS;
73  mDestCRS = projector.mDestCRS;
74  mMaxSrcXRes = projector.mMaxSrcXRes;
75  mMaxSrcYRes = projector.mMaxSrcYRes;
76  mExtent = projector.mExtent;
79 }
80 
82 {
83  if ( &projector != this )
84  {
85  mSrcCRS = projector.mSrcCRS;
86  mDestCRS = projector.mDestCRS;
87  mMaxSrcXRes = projector.mMaxSrcXRes;
88  mMaxSrcYRes = projector.mMaxSrcYRes;
89  mExtent = projector.mExtent;
92  }
93  return *this;
94 }
95 
97 {
98  QgsDebugMsg( "Entered" );
100  return projector;
101 }
102 
104 {
105  delete[] pHelperTop;
106  delete[] pHelperBottom;
107 }
108 
110 {
111  if ( mInput ) return mInput->bandCount();
112 
113  return 0;
114 }
115 
117 {
118  if ( mInput ) return mInput->dataType( bandNo );
119 
120  return QGis::UnknownDataType;
121 }
122 
124 {
125  mSrcCRS = theSrcCRS;
126  mDestCRS = theDestCRS;
127  mCoordinateTransform.setSourceCrs( theDestCRS );
128  mCoordinateTransform.setDestCRS( theSrcCRS );
129 }
130 
132 {
133  QgsDebugMsg( "Entered" );
134  mCPMatrix.clear();
135  mCPLegalMatrix.clear();
136  delete[] pHelperTop;
137  pHelperTop = 0;
138  delete[] pHelperBottom;
139  pHelperBottom = 0;
140 
141  // Get max source resolution and extent if possible
142  mMaxSrcXRes = 0;
143  mMaxSrcYRes = 0;
144  if ( mInput )
145  {
146  QgsRasterDataProvider *provider = dynamic_cast<QgsRasterDataProvider*>( mInput->srcInput() );
147  if ( provider && ( provider->capabilities() & QgsRasterDataProvider::Size ) )
148  {
149  mMaxSrcXRes = provider->extent().width() / provider->xSize();
150  mMaxSrcYRes = provider->extent().height() / provider->ySize();
151  }
152  // Get source extent
153  if ( mExtent.isEmpty() )
154  {
155  mExtent = provider->extent();
156  }
157  }
158 
161 
162  // Calculate tolerance
163  // TODO: Think it over better
164  // Note: we are checking on matrix each even point, that means that the real error
165  // in that moment is approximately half size
166  double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes;
167  mSqrTolerance = myDestRes * myDestRes;
168 
169  // Initialize the matrix by corners and middle points
170  mCPCols = mCPRows = 3;
171  for ( int i = 0; i < mCPRows; i++ )
172  {
173  QList<QgsPoint> myRow;
174  myRow.append( QgsPoint() );
175  myRow.append( QgsPoint() );
176  myRow.append( QgsPoint() );
177  mCPMatrix.insert( i, myRow );
178  // And the legal points
179  QList<bool> myLegalRow;
180  myLegalRow.append( bool( false ) );
181  myLegalRow.append( bool( false ) );
182  myLegalRow.append( bool( false ) );
183  mCPLegalMatrix.insert( i, myLegalRow );
184  }
185  for ( int i = 0; i < mCPRows; i++ )
186  {
187  calcRow( i );
188  }
189 
190  while ( true )
191  {
192  bool myColsOK = checkCols();
193  if ( !myColsOK )
194  {
195  insertRows();
196  }
197  bool myRowsOK = checkRows();
198  if ( !myRowsOK )
199  {
200  insertCols();
201  }
202  if ( myColsOK && myRowsOK )
203  {
204  QgsDebugMsg( "CP matrix within tolerance" );
205  mApproximate = true;
206  break;
207  }
208  // What is the maximum reasonable size of transformatio matrix?
209  // TODO: consider better when to break - ratio
210  if ( mCPRows * mCPCols > 0.25 * mDestRows * mDestCols )
211  {
212  QgsDebugMsg( "Too large CP matrix" );
213  mApproximate = false;
214  break;
215  }
216  }
217  QgsDebugMsg( QString( "CPMatrix size: mCPRows = %1 mCPCols = %2" ).arg( mCPRows ).arg( mCPCols ) );
218  mDestRowsPerMatrixRow = ( float )mDestRows / ( mCPRows - 1 );
219  mDestColsPerMatrixCol = ( float )mDestCols / ( mCPCols - 1 );
220 
221  QgsDebugMsgLevel( "CPMatrix:", 5 );
223 
224  // Calculate source dimensions
225  calcSrcExtent();
226  calcSrcRowsCols();
229 
230  // init helper points
233  calcHelper( 0, pHelperTop );
235  mHelperTopRow = 0;
236 }
237 
239 {
240  /* Run around the mCPMatrix and find source extent */
241  // Attention, source limits are not necessarily on destination edges, e.g.
242  // for destination EPSG:32661 Polar Stereographic and source EPSG:4326,
243  // the maximum y may be in the middle of destination extent
244  // TODO: How to find extent exactly and quickly?
245  // For now, we runt through all matrix
246  QgsPoint myPoint = mCPMatrix[0][0];
247  mSrcExtent = QgsRectangle( myPoint.x(), myPoint.y(), myPoint.x(), myPoint.y() );
248  for ( int i = 0; i < mCPRows; i++ )
249  {
250  for ( int j = 0; j < mCPCols ; j++ )
251  {
252  myPoint = mCPMatrix[i][j];
253  if ( mCPLegalMatrix[i][j] )
254  {
255  mSrcExtent.combineExtentWith( myPoint.x(), myPoint.y() );
256  }
257  }
258  }
259  // Expand a bit to avoid possible approx coords falling out because of representation error?
260 
261  // If mMaxSrcXRes, mMaxSrcYRes are defined (fixed src resolution)
262  // align extent to src resolution to avoid jumping of reprojected pixels
263  // when shifting resampled grid.
264  // Important especially if we are over mMaxSrcXRes, mMaxSrcYRes limits
265  // Note however, that preceding filters (like resampler) may read data
266  // on different resolution.
267 
268  QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
269  QgsDebugMsg( "mExtent = " + mExtent.toString() );
270  if ( !mExtent.isEmpty() )
271  {
272  if ( mMaxSrcXRes > 0 )
273  {
274  // with floor/ceil it should work correctly also for mSrcExtent.xMinimum() < mExtent.xMinimum()
275  double col = floor(( mSrcExtent.xMinimum() - mExtent.xMinimum() ) / mMaxSrcXRes );
276  double x = mExtent.xMinimum() + col * mMaxSrcXRes;
277  mSrcExtent.setXMinimum( x );
278 
279  col = ceil(( mSrcExtent.xMaximum() - mExtent.xMinimum() ) / mMaxSrcXRes );
280  x = mExtent.xMinimum() + col * mMaxSrcXRes;
281  mSrcExtent.setXMaximum( x );
282  }
283  if ( mMaxSrcYRes > 0 )
284  {
285  double row = floor(( mExtent.yMaximum() - mSrcExtent.yMaximum() ) / mMaxSrcYRes );
286  double y = mExtent.yMaximum() - row * mMaxSrcYRes;
287  mSrcExtent.setYMaximum( y );
288 
289  row = ceil(( mExtent.yMaximum() - mSrcExtent.yMinimum() ) / mMaxSrcYRes );
290  y = mExtent.yMaximum() - row * mMaxSrcYRes;
291  mSrcExtent.setYMinimum( y );
292  }
293  }
294  QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
295 }
296 
298 {
299  QString myString;
300  for ( int i = 0; i < mCPRows; i++ )
301  {
302  if ( i > 0 )
303  myString += "\n";
304  for ( int j = 0; j < mCPCols; j++ )
305  {
306  if ( j > 0 )
307  myString += " ";
308  QgsPoint myPoint = mCPMatrix[i][j];
309  if ( mCPLegalMatrix[i][j] )
310  {
311  myString += myPoint.toString();
312  }
313  else
314  {
315  myString += "(-,-)";
316  }
317  }
318  }
319  return myString;
320 }
321 
323 {
324  // Wee need to calculate minimum cell size in the source
325  // TODO: Think it over better, what is the right source resolution?
326  // Taking distances between cell centers projected to source along source
327  // axis would result in very high resolution
328  // TODO: different resolution for rows and cols ?
329 
330  // For now, we take cell sizes projected to source but not to source axes
331  double myDestColsPerMatrixCell = mDestCols / mCPCols;
332  double myDestRowsPerMatrixCell = mDestRows / mCPRows;
333  QgsDebugMsg( QString( "myDestColsPerMatrixCell = %1 myDestRowsPerMatrixCell = %2" ).arg( myDestColsPerMatrixCell ).arg( myDestRowsPerMatrixCell ) );
334 
335  double myMinSize = DBL_MAX;
336 
337  for ( int i = 0; i < mCPRows - 1; i++ )
338  {
339  for ( int j = 0; j < mCPCols - 1; j++ )
340  {
341  QgsPoint myPointA = mCPMatrix[i][j];
342  QgsPoint myPointB = mCPMatrix[i][j+1];
343  QgsPoint myPointC = mCPMatrix[i+1][j];
344  if ( mCPLegalMatrix[i][j] && mCPLegalMatrix[i][j+1] && mCPLegalMatrix[i+1][j] )
345  {
346  double mySize = sqrt( myPointA.sqrDist( myPointB ) ) / myDestColsPerMatrixCell;
347  if ( mySize < myMinSize )
348  myMinSize = mySize;
349 
350  mySize = sqrt( myPointA.sqrDist( myPointC ) ) / myDestRowsPerMatrixCell;
351  if ( mySize < myMinSize )
352  myMinSize = mySize;
353  }
354  }
355  }
356 
357  // Make it a bit higher resolution
358  // TODO: find the best coefficient, attention, increasing resolution for WMS
359  // is changing WMS content
360  myMinSize *= 0.75;
361 
362  QgsDebugMsg( QString( "mMaxSrcXRes = %1 mMaxSrcYRes = %2" ).arg( mMaxSrcXRes ).arg( mMaxSrcYRes ) );
363  // mMaxSrcXRes, mMaxSrcYRes may be 0 - no limit (WMS)
364  double myMinXSize = mMaxSrcXRes > myMinSize ? mMaxSrcXRes : myMinSize;
365  double myMinYSize = mMaxSrcYRes > myMinSize ? mMaxSrcYRes : myMinSize;
366  QgsDebugMsg( QString( "myMinXSize = %1 myMinYSize = %2" ).arg( myMinXSize ).arg( myMinYSize ) );
367  QgsDebugMsg( QString( "mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ) );
368 
369  // we have to round to keep alignment set in calcSrcExtent
370  mSrcRows = ( int ) qRound( mSrcExtent.height() / myMinYSize );
371  mSrcCols = ( int ) qRound( mSrcExtent.width() / myMinXSize );
372 
373  QgsDebugMsg( QString( "mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ) );
374 }
375 
376 
377 inline void QgsRasterProjector::destPointOnCPMatrix( int theRow, int theCol, double *theX, double *theY )
378 {
379  *theX = mDestExtent.xMinimum() + theCol * mDestExtent.width() / ( mCPCols - 1 );
380  *theY = mDestExtent.yMaximum() - theRow * mDestExtent.height() / ( mCPRows - 1 );
381 }
382 
383 inline int QgsRasterProjector::matrixRow( int theDestRow )
384 {
385  return ( int )( floor(( theDestRow + 0.5 ) / mDestRowsPerMatrixRow ) );
386 }
387 inline int QgsRasterProjector::matrixCol( int theDestCol )
388 {
389  return ( int )( floor(( theDestCol + 0.5 ) / mDestColsPerMatrixCol ) );
390 }
391 
392 QgsPoint QgsRasterProjector::srcPoint( int theDestRow, int theCol )
393 {
394  Q_UNUSED( theDestRow );
395  Q_UNUSED( theCol );
396  return QgsPoint();
397 }
398 
399 void QgsRasterProjector::calcHelper( int theMatrixRow, QgsPoint *thePoints )
400 {
401  // TODO?: should we also precalc dest cell center coordinates for x and y?
402  for ( int myDestCol = 0; myDestCol < mDestCols; myDestCol++ )
403  {
404  double myDestX = mDestExtent.xMinimum() + ( myDestCol + 0.5 ) * mDestXRes;
405 
406  int myMatrixCol = matrixCol( myDestCol );
407 
408  double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
409 
410  destPointOnCPMatrix( theMatrixRow, myMatrixCol, &myDestXMin, &myDestYMin );
411  destPointOnCPMatrix( theMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
412 
413  double xfrac = ( myDestX - myDestXMin ) / ( myDestXMax - myDestXMin );
414 
415  QgsPoint &mySrcPoint0 = mCPMatrix[theMatrixRow][myMatrixCol];
416  QgsPoint &mySrcPoint1 = mCPMatrix[theMatrixRow][myMatrixCol+1];
417  double s = mySrcPoint0.x() + ( mySrcPoint1.x() - mySrcPoint0.x() ) * xfrac;
418  double t = mySrcPoint0.y() + ( mySrcPoint1.y() - mySrcPoint0.y() ) * xfrac;
419 
420  thePoints[myDestCol].setX( s );
421  thePoints[myDestCol].setY( t );
422  }
423 }
425 {
426  // We just switch pHelperTop and pHelperBottom, memory is not lost
427  QgsPoint *tmp;
428  tmp = pHelperTop;
430  pHelperBottom = tmp;
432  mHelperTopRow++;
433 }
434 
435 void QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
436 {
437  if ( mApproximate )
438  {
439  approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
440  }
441  else
442  {
443  preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
444  }
445 }
446 
447 void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
448 {
449 #ifdef QGISDEBUG
450  QgsDebugMsgLevel( QString( "theDestRow = %1" ).arg( theDestRow ), 5 );
451  QgsDebugMsgLevel( QString( "theDestRow = %1 mDestExtent.yMaximum() = %2 mDestYRes = %3" ).arg( theDestRow ).arg( mDestExtent.yMaximum() ).arg( mDestYRes ), 5 );
452 #endif
453 
454  // Get coordinate of center of destination cell
455  double x = mDestExtent.xMinimum() + ( theDestCol + 0.5 ) * mDestXRes;
456  double y = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
457  double z = 0;
458 
459 #ifdef QGISDEBUG
460  QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
461 #endif
462 
464 
465 #ifdef QGISDEBUG
466  QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
467 #endif
468 
469  // Get source row col
470  *theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - y ) / mSrcYRes );
471  *theSrcCol = ( int ) floor(( x - mSrcExtent.xMinimum() ) / mSrcXRes );
472 #ifdef QGISDEBUG
473  QgsDebugMsgLevel( QString( "mSrcExtent.yMaximum() = %1 mSrcYRes = %2" ).arg( mSrcExtent.yMaximum() ).arg( mSrcYRes ), 5 );
474  QgsDebugMsgLevel( QString( "theSrcRow = %1 theSrcCol = %2" ).arg( *theSrcRow ).arg( *theSrcCol ), 5 );
475 #endif
476 
477  // With epsg 32661 (Polar Stereographic) it was happening that *theSrcCol == mSrcCols
478  // For now silently correct limits to avoid crashes
479  // TODO: review
480  if ( *theSrcRow >= mSrcRows )
481  *theSrcRow = mSrcRows - 1;
482  if ( *theSrcRow < 0 )
483  *theSrcRow = 0;
484  if ( *theSrcCol >= mSrcCols )
485  *theSrcCol = mSrcCols - 1;
486  if ( *theSrcCol < 0 )
487  *theSrcCol = 0;
488 
489  Q_ASSERT( *theSrcRow < mSrcRows );
490  Q_ASSERT( *theSrcCol < mSrcCols );
491 }
492 
493 void QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
494 {
495  int myMatrixRow = matrixRow( theDestRow );
496  int myMatrixCol = matrixCol( theDestCol );
497 
498  if ( myMatrixRow > mHelperTopRow )
499  {
500  // TODO: make it more robust (for random, not sequential reading)
501  nextHelper();
502  }
503 
504  double myDestY = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
505 
506  // See the schema in javax.media.jai.WarpGrid doc (but up side down)
507  // TODO: use some kind of cache of values which can be reused
508  double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
509 
510  destPointOnCPMatrix( myMatrixRow + 1, myMatrixCol, &myDestXMin, &myDestYMin );
511  destPointOnCPMatrix( myMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
512 
513  double yfrac = ( myDestY - myDestYMin ) / ( myDestYMax - myDestYMin );
514 
515  QgsPoint &myTop = pHelperTop[theDestCol];
516  QgsPoint &myBot = pHelperBottom[theDestCol];
517 
518  // Warning: this is very SLOW compared to the following code!:
519  //double mySrcX = myBot.x() + (myTop.x() - myBot.x()) * yfrac;
520  //double mySrcY = myBot.y() + (myTop.y() - myBot.y()) * yfrac;
521 
522  double tx = myTop.x();
523  double ty = myTop.y();
524  double bx = myBot.x();
525  double by = myBot.y();
526  double mySrcX = bx + ( tx - bx ) * yfrac;
527  double mySrcY = by + ( ty - by ) * yfrac;
528 
529  // TODO: check again cell selection (coor is in the middle)
530 
531  *theSrcRow = ( int ) floor(( mSrcExtent.yMaximum() - mySrcY ) / mSrcYRes );
532  *theSrcCol = ( int ) floor(( mySrcX - mSrcExtent.xMinimum() ) / mSrcXRes );
533 
534  // For now silently correct limits to avoid crashes
535  // TODO: review
536  if ( *theSrcRow >= mSrcRows )
537  *theSrcRow = mSrcRows - 1;
538  if ( *theSrcRow < 0 )
539  *theSrcRow = 0;
540  if ( *theSrcCol >= mSrcCols )
541  *theSrcCol = mSrcCols - 1;
542  if ( *theSrcCol < 0 )
543  *theSrcCol = 0;
544  Q_ASSERT( *theSrcRow < mSrcRows );
545  Q_ASSERT( *theSrcCol < mSrcCols );
546 }
547 
549 {
550  for ( int r = 0; r < mCPRows - 1; r++ )
551  {
552  QList<QgsPoint> myRow;
553  QList<bool> myLegalRow;
554  for ( int c = 0; c < mCPCols; c++ )
555  {
556  myRow.append( QgsPoint() );
557  myLegalRow.append( false );
558  }
559  QgsDebugMsgLevel( QString( "insert new row at %1" ).arg( 1 + r*2 ), 3 );
560  mCPMatrix.insert( 1 + r*2, myRow );
561  mCPLegalMatrix.insert( 1 + r*2, myLegalRow );
562  }
563  mCPRows += mCPRows - 1;
564  for ( int r = 1; r < mCPRows - 1; r += 2 )
565  {
566  calcRow( r );
567  }
568 }
569 
571 {
572  for ( int r = 0; r < mCPRows; r++ )
573  {
574  QList<QgsPoint> myRow;
575  QList<bool> myLegalRow;
576  for ( int c = 0; c < mCPCols - 1; c++ )
577  {
578  mCPMatrix[r].insert( 1 + c*2, QgsPoint() );
579  mCPLegalMatrix[r].insert( 1 + c*2, false );
580  }
581  }
582  mCPCols += mCPCols - 1;
583  for ( int c = 1; c < mCPCols - 1; c += 2 )
584  {
585  calcCol( c );
586  }
587 
588 }
589 
590 void QgsRasterProjector::calcCP( int theRow, int theCol )
591 {
592  double myDestX, myDestY;
593  destPointOnCPMatrix( theRow, theCol, &myDestX, &myDestY );
594  QgsPoint myDestPoint( myDestX, myDestY );
595  try
596  {
597  mCPMatrix[theRow][theCol] = mCoordinateTransform.transform( myDestPoint );
598  mCPLegalMatrix[theRow][theCol] = true;
599  }
600  catch ( QgsCsException &e )
601  {
602  Q_UNUSED( e );
603  // Caught an error in transform
604  mCPLegalMatrix[theRow][theCol] = false;
605  }
606 }
607 
608 bool QgsRasterProjector::calcRow( int theRow )
609 {
610  QgsDebugMsgLevel( QString( "theRow = %1" ).arg( theRow ), 3 );
611  for ( int i = 0; i < mCPCols; i++ )
612  {
613  calcCP( theRow, i );
614  }
615 
616  return true;
617 }
618 
619 bool QgsRasterProjector::calcCol( int theCol )
620 {
621  QgsDebugMsgLevel( QString( "theCol = %1" ).arg( theCol ), 3 );
622  for ( int i = 0; i < mCPRows; i++ )
623  {
624  calcCP( i, theCol );
625  }
626 
627  return true;
628 }
629 
631 {
632  for ( int c = 0; c < mCPCols; c++ )
633  {
634  for ( int r = 1; r < mCPRows - 1; r += 2 )
635  {
636  double myDestX, myDestY;
637  destPointOnCPMatrix( r, c, &myDestX, &myDestY );
638  QgsPoint myDestPoint( myDestX, myDestY );
639 
640  QgsPoint mySrcPoint1 = mCPMatrix[r-1][c];
641  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
642  QgsPoint mySrcPoint3 = mCPMatrix[r+1][c];
643 
644  QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
645  if ( !mCPLegalMatrix[r-1][c] || !mCPLegalMatrix[r][c] || !mCPLegalMatrix[r+1][c] )
646  {
647  // There was an error earlier in transform, just abort
648  return false;
649  }
650  try
651  {
653  double mySqrDist = myDestApprox.sqrDist( myDestPoint );
654  if ( mySqrDist > mSqrTolerance )
655  {
656  return false;
657  }
658  }
659  catch ( QgsCsException &e )
660  {
661  Q_UNUSED( e );
662  // Caught an error in transform
663  return false;
664  }
665  }
666  }
667  return true;
668 }
669 
671 {
672  for ( int r = 0; r < mCPRows; r++ )
673  {
674  for ( int c = 1; c < mCPCols - 1; c += 2 )
675  {
676  double myDestX, myDestY;
677  destPointOnCPMatrix( r, c, &myDestX, &myDestY );
678 
679  QgsPoint myDestPoint( myDestX, myDestY );
680  QgsPoint mySrcPoint1 = mCPMatrix[r][c-1];
681  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
682  QgsPoint mySrcPoint3 = mCPMatrix[r][c+1];
683 
684  QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
685  if ( !mCPLegalMatrix[r][c-1] || !mCPLegalMatrix[r][c] || !mCPLegalMatrix[r][c+1] )
686  {
687  // There was an error earlier in transform, just abort
688  return false;
689  }
690  try
691  {
693  double mySqrDist = myDestApprox.sqrDist( myDestPoint );
694  if ( mySqrDist > mSqrTolerance )
695  {
696  return false;
697  }
698  }
699  catch ( QgsCsException &e )
700  {
701  Q_UNUSED( e );
702  // Caught an error in transform
703  return false;
704  }
705  }
706  }
707  return true;
708 }
709 
710 QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & extent, int width, int height )
711 {
712  QgsDebugMsg( QString( "extent:\n%1" ).arg( extent.toString() ) );
713  QgsDebugMsg( QString( "width = %1 height = %2" ).arg( width ).arg( height ) );
714  if ( !mInput )
715  {
716  QgsDebugMsg( "Input not set" );
717  return new QgsRasterBlock();
718  }
719 
720  if ( ! mSrcCRS.isValid() || ! mDestCRS.isValid() || mSrcCRS == mDestCRS )
721  {
722  QgsDebugMsg( "No projection necessary" );
723  return mInput->block( bandNo, extent, width, height );
724  }
725 
727  mDestRows = height;
728  mDestCols = width;
729  calc();
730 
731  QgsDebugMsg( QString( "srcExtent:\n%1" ).arg( srcExtent().toString() ) );
732  QgsDebugMsg( QString( "srcCols = %1 srcRows = %2" ).arg( srcCols() ).arg( srcRows() ) );
733 
734  // If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers
735  if ( srcRows() <= 0 || srcCols() <= 0 )
736  {
737  QgsDebugMsg( "Zero srcRows or srcCols" );
738  return new QgsRasterBlock();
739  }
740 
741  QgsRasterBlock *inputBlock = mInput->block( bandNo, srcExtent(), srcCols(), srcRows() );
742  if ( !inputBlock || inputBlock->isEmpty() )
743  {
744  QgsDebugMsg( "No raster data!" );
745  delete inputBlock;
746  return new QgsRasterBlock();
747  }
748 
749  size_t pixelSize = QgsRasterBlock::typeSize( mInput->dataType( bandNo ) );
750 
751  QgsRasterBlock *outputBlock;
752  if ( inputBlock->hasNoDataValue() )
753  {
754  outputBlock = new QgsRasterBlock( inputBlock->dataType(), width, height, inputBlock->noDataValue() );
755  }
756  else
757  {
758  outputBlock = new QgsRasterBlock( inputBlock->dataType(), width, height );
759  }
760  if ( !outputBlock->isValid() )
761  {
762  QgsDebugMsg( "Cannot create block" );
763  delete inputBlock;
764  return outputBlock;
765  }
766 
767  // No data:
768  // 1) no data value exists (numerical) -> memcpy, not necessary isNoData()/setIsNoData()
769  // 2) no data value does not exist but it may contain no data (numerical no data bitmap)
770  // -> must use isNoData()/setIsNoData()
771  // 3) no data are not used (no no data value, no no data bitmap) -> simple memcpy
772  // 4) image - simple memcpy
773 
774  // To copy no data values stored in bitmaps we have to use isNoData()/setIsNoData(),
775  // we cannot fill output block with no data because we use memcpy for data, not setValue().
776  bool doNoData = inputBlock->hasNoData() && !inputBlock->hasNoDataValue();
777 
778  int srcRow, srcCol;
779  for ( int i = 0; i < height; ++i )
780  {
781  for ( int j = 0; j < width; ++j )
782  {
783  srcRowCol( i, j, &srcRow, &srcCol );
784  size_t srcIndex = ( size_t )srcRow * mSrcCols + srcCol;
785  QgsDebugMsgLevel( QString( "row = %1 col = %2 srcRow = %3 srcCol = %4" ).arg( i ).arg( j ).arg( srcRow ).arg( srcCol ), 5 );
786 
787  // isNoData() may be slow so we check doNoData first
788  if ( doNoData && inputBlock->isNoData( srcRow, srcCol ) )
789  {
790  outputBlock->setIsNoData( srcRow, srcCol );
791  continue ;
792  }
793 
794  size_t destIndex = ( size_t )i * width + j;
795  char *srcBits = inputBlock->bits( srcIndex );
796  char *destBits = outputBlock->bits( destIndex );
797  if ( !srcBits )
798  {
799  QgsDebugMsg( QString( "Cannot get input block data: row = %1 col = %2" ).arg( i ).arg( j ) );
800  continue;
801  }
802  if ( !destBits )
803  {
804  QgsDebugMsg( QString( "Cannot set output block data: srcRow = %1 srcCol = %2" ).arg( srcRow ).arg( srcCol ) );
805  continue;
806  }
807  memcpy( destBits, srcBits, pixelSize );
808  }
809  }
810 
811  delete inputBlock;
812 
813  return outputBlock;
814 }