QGIS中的金字塔机制优化

QGIS中可以通过图层的属性窗口来给栅格数据创建金字塔，如下图：

右侧的Pyramidresolutions表示可以建立的金字塔的级别，数字前面的符号有个小叉就表示该级别的金字塔还没有被建立。可以通过点击右下方的Build pyramid按钮来为影像建立金字塔。

QGIS内部也是调用了GDAL的GDALBuildOverviews方法来为影像建立金字塔。但在使用Build pyramid建立金字塔的过程中发现速度非常慢，与使用gdaladdo.exe命令来建立金字塔相比，要耗费它几倍的时间，对此感到非常疑惑，于是看了一下内部的实现机制。

建立金字塔的实现函数位于QgsRasterLayer中，函数名为buildPyramids，定义如下：

/** /brief Create GDAL pyramid overviews */

    QStringbuildPyramids( constRasterPyramidList &,

                           const QString&  theResamplingMethod= "NEAREST",

                           bool theTryInternalFlag= false );

RasterPyramidList参数表示要建立金字塔的级别。

theResamplingMethod参数表示重采样方法。

theTryInternalFlag参数表示是否将金字塔嵌入影像内部，该操作好像只支持TIFF格式。

找到了函数，查看其内部实现，很快就找到了原因所在。

以下代码是QGIS建立金字塔的核心代码

RasterPyramidList::const_iteratormyRasterPyramidIterator;

  for( myRasterPyramidIterator = theRasterPyramidList.begin();

        myRasterPyramidIterator!= theRasterPyramidList.end();

        ++myRasterPyramidIterator)

  {

    if(( *myRasterPyramidIterator ).build )

    {

      QgsDebugMsg("Building....." );

      emitdrawingProgress( myCount,myTotal );

      intmyOverviewLevelsArray[1] = {( *myRasterPyramidIterator ).level};

      /*From : http://remotesensing.org/gdal/classGDALDataset.html#a23

       * pszResampling: one of "NEAREST", "AVERAGE" or "MODE"controlling the downsampling method applied.

       * nOverviews :number of overviews to build.

       *panOverviewList : the list of overview decimation factors to build.

       * nBand : numberof bands to build overviews for in panBandList. Build for all bands if this is0.

       * panBandList :list of band numbers.

       * pfnProgress :a function to call to report progress, or NULL.

       * pProgressData: application data to pass to the progress function.

       */

      //buildthe pyramid and show progress to console

      try

      {

 

        //buildthe pyramid and show progress to console

        //NOTEthis (magphase) is disabled in teh gui since it tends

        //tocreate corrupted images. The images can be repaired

        //byrunning one of the other resampling strategies below.

        //seeticket #284

        if( theResamplingMethod == tr( "AverageMagphase" ) )

        {

          myError= GDALBuildOverviews( mGdalBaseDataset, "MODE",1, myOverviewLevelsArray, 0, NULL,

                                       progressCallback, this ); //this is the argfor the gdal progress callback

        }

        elseif ( theResamplingMethod== tr( "Average") )

 

        {

          myError= GDALBuildOverviews( mGdalBaseDataset, "AVERAGE",1, myOverviewLevelsArray, 0, NULL,

                                       progressCallback, this ); //this is the argfor the gdal progress callback

        }

        else// fall back to nearest neighbor

        {

          myError= GDALBuildOverviews( mGdalBaseDataset, "NEAREST",1, myOverviewLevelsArray, 0, NULL,

                                       progressCallback, this ); //this is the argfor the gdal progress callback

        }

 

        if( myError == CE_Failure|| CPLGetLastErrorNo() == CPLE_NotSupported )

        {

          //somethingbad happenend

          //QStringmyString = QString (CPLGetLastError());

          GDALClose(mGdalBaseDataset );

          mGdalBaseDataset= GDALOpen( QFile::encodeName( mDataSource).constData(), GA_ReadOnly);

          //Sincewe are not a virtual warped dataset, mGdalDataSet and mGdalBaseDataset aresupposed to be the same

          mGdalDataset= mGdalBaseDataset;

 

          emitdrawingProgress( 0, 0 );

          return"FAILED_NOT_SUPPORTED";

        }

        else

        {

          //makesure the raster knows it has pyramids

          mHasPyramids= true;

        }

        myCount++;

 

      }

      catch( CPLErr )

      {

        QgsLogger::warning( "Pyramidoverview building failed!" );

      }

    }

可以看到，QGIS将建立金字塔的关键函数GDALBuildOverviews放到了一个for循环中。也就是说每建立一级金字塔，都会调用一次GDALBuildOverviews函数。

对此我使用gdaladdo命令做了一个验证，使用GDALBuildOverviews函数一次建立4级金字塔和分四次建立一个金字塔耗费的时间差不多是1:4。这就是速度慢的问题所在，不清楚QGIS的开发团队为什么要用这种方法来实现。

我对代码进行了一些修改，将GDALBuildOverviews函数移到了for循环之外。然后重新在属性面板中建立金字塔，速度提高很多。

修改的时侯要注意的核心代码如下：

     原代码：

int myOverviewLevelsArray[1]= {( *myRasterPyramidIterator ).level };

myError = GDALBuildOverviews(mGdalBaseDataset, "NEAREST",1, myOverviewLevelsArray, 0, NULL, progressCallback,this );

修改后：假如5级金字塔

Int iLevels = 5;

int myOverviewLevelsArray[5];

在循环中对myOverviewLevelsArray分别赋值：( *myRasterPyramidIterator ).level

myError = GDALBuildOverviews( mGdalBaseDataset, "NEAREST",iLevels, myOverviewLevelsArray,0, NULL，progressCallback, this );  

重点就是红色部分的两个变量修改，希望能对大家有所帮助。