QGIS API Documentation  2.0.1-Dufour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgsrelief.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrelief.cpp - description
3  ---------------------------
4  begin : November 2011
5  copyright : (C) 2011 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
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 "qgsrelief.h"
19 #include "qgsaspectfilter.h"
20 #include "qgshillshadefilter.h"
21 #include "qgsslopefilter.h"
22 #include "qgis.h"
23 #include "cpl_string.h"
24 #include <QProgressDialog>
25 #include <cfloat>
26 
27 #include <QFile>
28 #include <QTextStream>
29 
30 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
31 #define TO8F(x) (x).toUtf8().constData()
32 #else
33 #define TO8F(x) QFile::encodeName( x ).constData()
34 #endif
35 
36 QgsRelief::QgsRelief( const QString& inputFile, const QString& outputFile, const QString& outputFormat ): \
37  mInputFile( inputFile ), mOutputFile( outputFile ), mOutputFormat( outputFormat ), mZFactor( 1.0 )
38 {
39  mSlopeFilter = new QgsSlopeFilter( inputFile, outputFile, outputFormat );
40  mAspectFilter = new QgsAspectFilter( inputFile, outputFile, outputFormat );
41  mHillshadeFilter285 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 285, 30 );
42  mHillshadeFilter300 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 300, 30 );
43  mHillshadeFilter315 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 315, 30 );
44 
45  /*mReliefColors = calculateOptimizedReliefClasses();
46  setDefaultReliefColors();*/
47 }
48 
50 {
51  delete mSlopeFilter;
52  delete mAspectFilter;
53  delete mHillshadeFilter285;
54  delete mHillshadeFilter300;
55  delete mHillshadeFilter315;
56 }
57 
59 {
60  mReliefColors.clear();
61 }
62 
64 {
65  mReliefColors.push_back( color );
66 }
67 
69 {
71  addReliefColorClass( ReliefColor( QColor( 9, 176, 76 ), 0, 200 ) );
72  addReliefColorClass( ReliefColor( QColor( 20, 228, 128 ), 200, 500 ) );
73  addReliefColorClass( ReliefColor( QColor( 167, 239, 153 ), 500, 1000 ) );
74  addReliefColorClass( ReliefColor( QColor( 218, 188, 143 ), 1000, 2000 ) );
75  addReliefColorClass( ReliefColor( QColor( 233, 158, 91 ), 2000, 4000 ) );
76  addReliefColorClass( ReliefColor( QColor( 255, 255, 255 ), 4000, 9000 ) );
77 }
78 
79 int QgsRelief::processRaster( QProgressDialog* p )
80 {
81  //open input file
82  int xSize, ySize;
83  GDALDatasetH inputDataset = openInputFile( xSize, ySize );
84  if ( inputDataset == NULL )
85  {
86  return 1; //opening of input file failed
87  }
88 
89  //output driver
90  GDALDriverH outputDriver = openOutputDriver();
91  if ( outputDriver == 0 )
92  {
93  return 2;
94  }
95 
96  GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
97  if ( outputDataset == NULL )
98  {
99  return 3; //create operation on output file failed
100  }
101 
102  //initialize dependency filters with cell sizes
118 
119  //open first raster band for reading (operation is only for single band raster)
120  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
121  if ( rasterBand == NULL )
122  {
123  GDALClose( inputDataset );
124  GDALClose( outputDataset );
125  return 4;
126  }
127  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
133 
134  GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset, 1 );
135  GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset, 2 );
136  GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset, 3 );
137 
138  if ( outputRedBand == NULL || outputGreenBand == NULL || outputBlueBand == NULL )
139  {
140  GDALClose( inputDataset );
141  GDALClose( outputDataset );
142  return 5;
143  }
144  //try to set -9999 as nodata value
145  GDALSetRasterNoDataValue( outputRedBand, -9999 );
146  GDALSetRasterNoDataValue( outputGreenBand, -9999 );
147  GDALSetRasterNoDataValue( outputBlueBand, -9999 );
148  mOutputNodataValue = GDALGetRasterNoDataValue( outputRedBand, NULL );
154 
155  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
156  {
157  GDALClose( inputDataset );
158  GDALClose( outputDataset );
159  return 6;
160  }
161 
162  //keep only three scanlines in memory at a time
163  float* scanLine1 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
164  float* scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
165  float* scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
166 
167  unsigned char* resultRedLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize );
168  unsigned char* resultGreenLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize );
169  unsigned char* resultBlueLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize );
170 
171  if ( p )
172  {
173  p->setMaximum( ySize );
174  }
175 
176  bool resultOk;
177 
178  //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
179  for ( int i = 0; i < ySize; ++i )
180  {
181  if ( p )
182  {
183  p->setValue( i );
184  }
185 
186  if ( p && p->wasCanceled() )
187  {
188  break;
189  }
190 
191  if ( i == 0 )
192  {
193  //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
194  for ( int a = 0; a < xSize; ++a )
195  {
196  scanLine1[a] = mInputNodataValue;
197  }
198  GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 );
199  }
200  else
201  {
202  //normally fetch only scanLine3 and release scanline 1 if we move forward one row
203  CPLFree( scanLine1 );
204  scanLine1 = scanLine2;
205  scanLine2 = scanLine3;
206  scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
207  }
208 
209  if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
210  {
211  for ( int a = 0; a < xSize; ++a )
212  {
213  scanLine3[a] = mInputNodataValue;
214  }
215  }
216  else
217  {
218  GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 );
219  }
220 
221  for ( int j = 0; j < xSize; ++j )
222  {
223  if ( j == 0 )
224  {
225  resultOk = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j+1], &mInputNodataValue, &scanLine2[j], \
226  &scanLine2[j+1], &mInputNodataValue, &scanLine3[j], &scanLine3[j+1], \
227  &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
228  }
229  else if ( j == xSize - 1 )
230  {
231  resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &mInputNodataValue, &scanLine2[j-1], &scanLine2[j], \
232  &mInputNodataValue, &scanLine3[j-1], &scanLine3[j], &mInputNodataValue, \
233  &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
234  }
235  else
236  {
237  resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &scanLine1[j+1], &scanLine2[j-1], &scanLine2[j], \
238  &scanLine2[j+1], &scanLine3[j-1], &scanLine3[j], &scanLine3[j+1], \
239  &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
240  }
241 
242  if ( !resultOk )
243  {
244  resultRedLine[j] = mOutputNodataValue;
245  resultGreenLine[j] = mOutputNodataValue;
246  resultBlueLine[j] = mOutputNodataValue;
247  }
248  }
249 
250  GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 );
251  GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 );
252  GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 );
253  }
254 
255  if ( p )
256  {
257  p->setValue( ySize );
258  }
259 
260  CPLFree( resultRedLine );
261  CPLFree( resultBlueLine );
262  CPLFree( resultGreenLine );
263  CPLFree( scanLine1 );
264  CPLFree( scanLine2 );
265  CPLFree( scanLine3 );
266 
267  GDALClose( inputDataset );
268 
269  if ( p && p->wasCanceled() )
270  {
271  //delete the dataset without closing (because it is faster)
272  GDALDeleteDataset( outputDriver, TO8F( mOutputFile ) );
273  return 7;
274  }
275  GDALClose( outputDataset );
276 
277  return 0;
278 }
279 
280 bool QgsRelief::processNineCellWindow( float* x1, float* x2, float* x3, float* x4, float* x5, float* x6, float* x7, float* x8, float* x9,
281  unsigned char* red, unsigned char* green, unsigned char* blue )
282 {
283  //1. component: color and hillshade from 300 degrees
284  int r = 0;
285  int g = 0;
286  int b = 0;
287 
288  float hillShadeValue300 = mHillshadeFilter300->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
289  if ( hillShadeValue300 != mOutputNodataValue )
290  {
291  if ( !setElevationColor( *x5, &r, &g, &b ) )
292  {
293  r = hillShadeValue300;
294  g = hillShadeValue300;
295  b = hillShadeValue300;
296  }
297  else
298  {
299  r = r / 2.0 + hillShadeValue300 / 2.0;
300  g = g / 2.0 + hillShadeValue300 / 2.0;
301  b = b / 2.0 + hillShadeValue300 / 2.0;
302  }
303  }
304 
305  //2. component: hillshade and slope
306  float hillShadeValue315 = mHillshadeFilter315->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
307  float slope = mSlopeFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
308  if ( hillShadeValue315 != mOutputNodataValue && slope != mOutputNodataValue )
309  {
310  int r2, g2, b2;
311  if ( slope > 15 )
312  {
313  r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
314  g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
315  b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
316  }
317  else if ( slope >= 1 )
318  {
319  int slopeValue = 255 - ( slope / 15.0 * 255.0 );
320  r2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
321  g2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
322  b2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
323  }
324  else
325  {
326  r2 = hillShadeValue315; g2 = hillShadeValue315; b2 = hillShadeValue315;
327  }
328 
329  //combine with r,g,b with 70 percentage coverage
330  r = r * 0.7 + r2 * 0.3;
331  g = g * 0.7 + g2 * 0.3;
332  b = b * 0.7 + b2 * 0.3;
333  }
334 
335  //3. combine yellow aspect with 10% transparency, illumination from 285 degrees
336  float hillShadeValue285 = mHillshadeFilter285->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
337  float aspect = mAspectFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
338  if ( hillShadeValue285 != mOutputNodataValue && aspect != mOutputNodataValue )
339  {
340  double angle_diff = qAbs( 285 - aspect );
341  if ( angle_diff > 180 )
342  {
343  angle_diff -= 180;
344  }
345 
346  int r3, g3, b3;
347  if ( angle_diff < 90 )
348  {
349  int aspectVal = ( 1 - cos( angle_diff * M_PI / 180 ) ) * 255;
350  r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
351  g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
352  b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5;
353  }
354  else //white
355  {
356  r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
357  g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
358  b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
359  }
360 
361  r = r3 * 0.1 + r * 0.9;
362  g = g3 * 0.1 + g * 0.9;
363  b = b3 * 0.1 + b * 0.9;
364  }
365 
366  *red = ( unsigned char )r;
367  *green = ( unsigned char )g;
368  *blue = ( unsigned char )b;
369  return true;
370 }
371 
372 bool QgsRelief::setElevationColor( double elevation, int* red, int* green, int* blue )
373 {
374  QList< ReliefColor >::const_iterator reliefColorIt = mReliefColors.constBegin();
375  for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
376  {
377  if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
378  {
379  const QColor& c = reliefColorIt->color;
380  *red = c.red();
381  *green = c.green();
382  *blue = c.blue();
383 
384  return true;
385  }
386  }
387  return false;
388 }
389 
390 //duplicated from QgsNineCellFilter. Todo: make common base class
391 GDALDatasetH QgsRelief::openInputFile( int& nCellsX, int& nCellsY )
392 {
393  GDALDatasetH inputDataset = GDALOpen( TO8F( mInputFile ), GA_ReadOnly );
394  if ( inputDataset != NULL )
395  {
396  nCellsX = GDALGetRasterXSize( inputDataset );
397  nCellsY = GDALGetRasterYSize( inputDataset );
398 
399  //we need at least one band
400  if ( GDALGetRasterCount( inputDataset ) < 1 )
401  {
402  GDALClose( inputDataset );
403  return NULL;
404  }
405  }
406  return inputDataset;
407 }
408 
410 {
411  char **driverMetadata;
412 
413  //open driver
414  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
415 
416  if ( outputDriver == NULL )
417  {
418  return outputDriver; //return NULL, driver does not exist
419  }
420 
421  driverMetadata = GDALGetMetadata( outputDriver, NULL );
422  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
423  {
424  return NULL; //driver exist, but it does not support the create operation
425  }
426 
427  return outputDriver;
428 }
429 
430 GDALDatasetH QgsRelief::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
431 {
432  if ( inputDataset == NULL )
433  {
434  return NULL;
435  }
436 
437  int xSize = GDALGetRasterXSize( inputDataset );
438  int ySize = GDALGetRasterYSize( inputDataset );;
439 
440  //open output file
441  char **papszOptions = NULL;
442 
443  //use PACKBITS compression for tiffs by default
444  papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
445 
446  //create three band raster (reg, green, blue)
447  GDALDatasetH outputDataset = GDALCreate( outputDriver, TO8F( mOutputFile ), xSize, ySize, 3, GDT_Byte, papszOptions );
448  if ( outputDataset == NULL )
449  {
450  return outputDataset;
451  }
452 
453  //get geotransform from inputDataset
454  double geotransform[6];
455  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
456  {
457  GDALClose( outputDataset );
458  return NULL;
459  }
460  GDALSetGeoTransform( outputDataset, geotransform );
461 
462  //make sure mCellSizeX and mCellSizeY are always > 0
463  mCellSizeX = geotransform[1];
464  if ( mCellSizeX < 0 )
465  {
467  }
468  mCellSizeY = geotransform[5];
469  if ( mCellSizeY < 0 )
470  {
472  }
473 
474  const char* projection = GDALGetProjectionRef( inputDataset );
475  GDALSetProjection( outputDataset, projection );
476 
477  return outputDataset;
478 }
479 
480 //this function is mainly there for debugging
482 {
483  int nCellsX, nCellsY;
484  GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
485  if ( inputDataset == NULL )
486  {
487  return false;
488  }
489 
490  //open first raster band for reading (elevation raster is always single band)
491  GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
492  if ( elevationBand == NULL )
493  {
494  GDALClose( inputDataset );
495  return false;
496  }
497 
498  //1. get minimum and maximum of elevation raster -> 252 elevation classes
499  int minOk, maxOk;
500  double minMax[2];
501  minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
502  minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
503 
504  if ( !minOk || !maxOk )
505  {
506  GDALComputeRasterMinMax( elevationBand, TRUE, minMax );
507  }
508 
509  //2. go through raster cells and get frequency of classes
510 
511  //store elevation frequency in 256 elevation classes
512  double frequency[252];
513  double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
514  //initialize to zero
515  for ( int i = 0; i < 252; ++i )
516  {
517  frequency[i] = 0;
518  }
519 
520  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
521  int elevationClass = -1;
522 
523  for ( int i = 0; i < nCellsY; ++i )
524  {
525  GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
526  scanLine, nCellsX, 1, GDT_Float32,
527  0, 0 );
528  for ( int j = 0; j < nCellsX; ++j )
529  {
530  elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
531  if ( elevationClass >= 0 )
532  {
533  frequency[elevationClass] += 1.0;
534  }
535  }
536  }
537 
538  CPLFree( scanLine );
539 
540  //log10 transformation for all frequency values
541  for ( int i = 0; i < 252; ++i )
542  {
543  frequency[i] = log10( frequency[i] );
544  }
545 
546  //write out frequency values to csv file for debugging
547  QFile outFile( file );
548  if ( !outFile.open( QIODevice::WriteOnly ) )
549  {
550  return false;
551  }
552 
553  QTextStream outstream( &outFile );
554  for ( int i = 0; i < 252; ++i )
555  {
556  outstream << QString::number( i ) + "," + QString::number( frequency[i] ) << endl;
557  }
558  outFile.close();
559  return true;
560 }
561 
562 QList< QgsRelief::ReliefColor > QgsRelief::calculateOptimizedReliefClasses()
563 {
564  QList< QgsRelief::ReliefColor > resultList;
565 
566  int nCellsX, nCellsY;
567  GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
568  if ( inputDataset == NULL )
569  {
570  return resultList;
571  }
572 
573  //open first raster band for reading (elevation raster is always single band)
574  GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
575  if ( elevationBand == NULL )
576  {
577  GDALClose( inputDataset );
578  return resultList;
579  }
580 
581  //1. get minimum and maximum of elevation raster -> 252 elevation classes
582  int minOk, maxOk;
583  double minMax[2];
584  minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
585  minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
586 
587  if ( !minOk || !maxOk )
588  {
589  GDALComputeRasterMinMax( elevationBand, TRUE, minMax );
590  }
591 
592  //2. go through raster cells and get frequency of classes
593 
594  //store elevation frequency in 256 elevation classes
595  double frequency[252];
596  double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
597  //initialize to zero
598  for ( int i = 0; i < 252; ++i )
599  {
600  frequency[i] = 0;
601  }
602 
603  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
604  int elevationClass = -1;
605 
606  for ( int i = 0; i < nCellsY; ++i )
607  {
608  GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
609  scanLine, nCellsX, 1, GDT_Float32,
610  0, 0 );
611  for ( int j = 0; j < nCellsX; ++j )
612  {
613  elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
614  if ( elevationClass < 0 )
615  {
616  elevationClass = 0;
617  }
618  else if ( elevationClass >= 252 )
619  {
620  elevationClass = 251;
621  }
622  frequency[elevationClass] += 1.0;
623  }
624  }
625 
626  CPLFree( scanLine );
627 
628  //log10 transformation for all frequency values
629  for ( int i = 0; i < 252; ++i )
630  {
631  frequency[i] = log10( frequency[i] );
632  }
633 
634  //start with 9 uniformly distributed classes
635  QList<int> classBreaks;
636  classBreaks.append( 0 );
637  classBreaks.append( 28 );
638  classBreaks.append( 56 );
639  classBreaks.append( 84 );
640  classBreaks.append( 112 );
641  classBreaks.append( 140 );
642  classBreaks.append( 168 );
643  classBreaks.append( 196 );
644  classBreaks.append( 224 );
645  classBreaks.append( 252 );
646 
647  for ( int i = 0; i < 10; ++i )
648  {
649  optimiseClassBreaks( classBreaks, frequency );
650  }
651 
652  //debug, print out all the classbreaks
653  for ( int i = 0; i < classBreaks.size(); ++i )
654  {
655  qWarning( "%d", classBreaks[i] );
656  }
657 
658  //set colors according to optimised class breaks
659  QList<QColor> colorList;
660  colorList.push_back( QColor( 7, 165, 144 ) );
661  colorList.push_back( QColor( 12, 221, 162 ) );
662  colorList.push_back( QColor( 33, 252, 183 ) );
663  colorList.push_back( QColor( 247, 252, 152 ) );
664  colorList.push_back( QColor( 252, 196, 8 ) );
665  colorList.push_back( QColor( 252, 166, 15 ) );
666  colorList.push_back( QColor( 175, 101, 15 ) );
667  colorList.push_back( QColor( 255, 133, 92 ) );
668  colorList.push_back( QColor( 204, 204, 204 ) );
669 
670  for ( int i = 1; i < classBreaks.size(); ++i )
671  {
672  double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
673  double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
674  resultList.push_back( QgsRelief::ReliefColor( colorList.at( i - 1 ), minElevation, maxElevation ) );
675  }
676 
677  return resultList;
678 }
679 
680 void QgsRelief::optimiseClassBreaks( QList<int>& breaks, double* frequencies )
681 {
682  int nClasses = breaks.size() - 1;
683  double* a = new double[nClasses]; //slopes
684  double* b = new double[nClasses]; //y-offsets
685 
686  for ( int i = 0; i < nClasses; ++i )
687  {
688  //get all the values between the class breaks into input
689  QList< QPair < int, double > > regressionInput;
690  for ( int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j )
691  {
692  regressionInput.push_back( qMakePair( j, frequencies[j] ) );
693  }
694 
695  double aParam, bParam;
696  if ( regressionInput.size() > 0 && calculateRegression( regressionInput, aParam, bParam ) )
697  {
698  a[i] = aParam;
699  b[i] = bParam;
700  }
701  else
702  {
703  a[i] = 0;
704  b[i] = 0; //better default value
705  }
706  }
707 
708  QList<int> classesToRemove;
709 
710  //shift class boundaries or eliminate classes which fall together
711  for ( int i = 1; i < nClasses ; ++i )
712  {
713  if ( breaks[i] == breaks[ i - 1 ] )
714  {
715  continue;
716  }
717 
718  if ( qgsDoubleNear( a[i - 1 ], a[i] ) )
719  {
720  continue;
721  }
722  else
723  {
724  int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
725 
726  if ( newX <= breaks[i - 1] )
727  {
728  newX = breaks[i - 1];
729  // classesToRemove.push_back( i );//remove this class later as it falls together with the preceding one
730  }
731  else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
732  {
733  newX = breaks[i + 1];
734  // classesToRemove.push_back( i );//remove this class later as it falls together with the next one
735  }
736 
737  breaks[i] = newX;
738  }
739  }
740 
741  for ( int i = classesToRemove.size() - 1; i >= 0; --i )
742  {
743  breaks.removeAt( classesToRemove.at( i ) );
744  }
745 
746  delete[] a;
747  delete[] b;
748 }
749 
750 int QgsRelief::frequencyClassForElevation( double elevation, double minElevation, double elevationClassRange )
751 {
752  return ( elevation - minElevation ) / elevationClassRange;
753 }
754 
755 bool QgsRelief::calculateRegression( const QList< QPair < int, double > >& input, double& a, double& b )
756 {
757  double xMean, yMean;
758  double xSum = 0;
759  double ySum = 0;
760  QList< QPair < int, double > >::const_iterator inputIt = input.constBegin();
761  for ( ; inputIt != input.constEnd(); ++inputIt )
762  {
763  xSum += inputIt->first;
764  ySum += inputIt->second;
765  }
766  xMean = xSum / input.size();
767  yMean = ySum / input.size();
768 
769  double sumCounter = 0;
770  double sumDenominator = 0;
771  inputIt = input.constBegin();
772  for ( ; inputIt != input.constEnd(); ++inputIt )
773  {
774  sumCounter += (( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
775  sumDenominator += (( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
776  }
777 
778  a = sumCounter / sumDenominator;
779  b = yMean - a * xMean;
780 
781  return true;
782 }
783