void object::test<6>()
    {
        // Index of test file being tested
        const std::size_t fileIdx = 1;

        std::string file(data_ + SEP);
        file += grids_.at(fileIdx).file_;

        GDALDatasetH ds = GDALOpen(file.c_str(), GA_ReadOnly);
        ensure("Can't open dataset: " + file, nullptr != ds);

        GDALRasterBandH band = GDALGetRasterBand(ds, grids_.at(fileIdx).band_);
        ensure("Can't get raster band", nullptr != band);

        const double noData = GDALGetRasterNoDataValue(band, nullptr);
        ensure_equals("Grid NODATA value wrong or missing", noData, -99999);

        ensure_equals("Data type is not GDT_Float32", GDALGetRasterDataType(band), GDT_Float32);

        GDALClose(ds);
    }

void generic_read(struct FI *f, char *filename)
{
	if (filename_corresponds_to_tiffo(filename)) {
#ifdef FANCY_TIFF
		f->tiffo = true;
		tiff_octaves_init0(f->t, filename, f->megabytes,f->max_octaves);
		if (f->option_write) f->t->option_write = true;
		f->w = f->t->i->w;
		f->h = f->t->i->h;
		f->pd = f->t->i->spp;
		f->no = f->t->noctaves;
#else
		assert(false);
#endif
	} else if (!f->option_write && FORCE_GDAL()) {
#ifdef FANCY_GDAL
		f->gdal = true;
		GDALAllRegister();
		char buf[2*FILENAME_MAX];
		snprintf(buf, 2*FILENAME_MAX, has_prefix(filename, "http://") ||
				has_prefix(filename, "https://") ?
				"/vsicurl/%s" : "%s", filename);
		f->gdal_img = GDALOpen(buf, GA_ReadOnly);
		fprintf(stderr, "gdal_dataset = %p\n", f->gdal_img);
		f->pd = GDALGetRasterCount(f->gdal_img);
		f->w = GDALGetRasterXSize(f->gdal_img);
		f->h = GDALGetRasterYSize(f->gdal_img);
		f->no = 1;
		for (int i = 0; i < f->pd; i++)
			f->gdal_band[i] = GDALGetRasterBand(f->gdal_img, i+1);
#else
		assert(false);
#endif
	} else {
		f->x = iio_read_image_float_vec(filename, &f->w, &f->h, &f->pd);
		f->no = build_pyramid(f, f->max_octaves);
		snprintf(f->x_filename, FILENAME_MAX, "%s", filename);
		f->x_changed = false;
	}
}

HRESULT CImage::SetPixel(int nRows, int nCols, BYTE* pPixel)
{
	if(m_pGdalImage==NULL)
	{
		return S_FALSE;
	}

	if(m_bFlip)
	{
		nRows=m_nRows-nRows-1;
	}

	if(nRows<0||nRows>=m_nRows||nCols<0||nCols>=m_nCols)
	{
		return S_OK;
	}

	for(int i=0;i<m_nBandNum;i++)
	{
		int band=i+1;

		GDALDataType datatype=m_datatype;
		int nBPB=m_nBPB;
		
		BYTE* pGray=new BYTE[nBPB];
		memcpy(pGray,pPixel+i*nBPB,nBPB);
		((GDALRasterBand*)GDALGetRasterBand(m_pGdalImage,band))->
		RasterIO(GF_Write,
				 nCols,nRows,
				 1,1,
				 pGray,
				 1,1,
				 datatype,
				 0,0);

		delete [] pGray; pGray=NULL;
	}

	return S_OK;
}

void writeGeoTiffF(char * fileName, float * result, int nRow, int nCol, double xMin, double yMax, double cellSize)
{
	
	GDALAllRegister();
	OGRRegisterAll();

	GDALDatasetH hDstDS;
	GDALDriverH hDriver;
	GDALRasterBandH hBand;
	double adfGeoTransform[6];

	char *papszOptions[] = {"COMPRESS=LZW",NULL};
	const char *pszFormat="GTiff";

	if(NULL == (hDriver = GDALGetDriverByName(pszFormat)))
	{
		printf("ERROR: hDriver is null cannot output using GDAL\n");
		exit(1);
	}
	
	hDstDS = GDALCreate(hDriver, fileName, nCol, nRow, 1, GDT_Float32, papszOptions);

	adfGeoTransform[0] = xMin;
	adfGeoTransform[1] = cellSize;
	adfGeoTransform[2] = 0;
	adfGeoTransform[3] = yMax;
	adfGeoTransform[4] = 0;
	adfGeoTransform[5] = -cellSize;

	GDALSetGeoTransform(hDstDS,adfGeoTransform);

	hBand=GDALGetRasterBand(hDstDS,1);
	GDALSetRasterNoDataValue(hBand,-1);
	GDALRasterIO(hBand, GF_Write, 0, 0, nCol, nRow, result, nCol, nRow, GDT_Float32, 0, 0 );
	
	GDALClose(hDstDS);

	return;
}

GDALDatasetH CPL_STDCALL
GDALCreateWarpedVRT( GDALDatasetH hSrcDS, 
                     int nPixels, int nLines, double *padfGeoTransform,
                     GDALWarpOptions *psOptions )
    
{
    VALIDATE_POINTER1( hSrcDS, "GDALCreateWarpedVRT", NULL );

/* -------------------------------------------------------------------- */
/*      Create the VRTDataset and populate it with bands.               */
/* -------------------------------------------------------------------- */
    VRTWarpedDataset *poDS = new VRTWarpedDataset( nPixels, nLines );
    int i;

    psOptions->hDstDS = (GDALDatasetH) poDS;

    poDS->SetGeoTransform( padfGeoTransform );

    for( i = 0; i < psOptions->nBandCount; i++ )
    {
        VRTWarpedRasterBand *poBand;
        GDALRasterBand *poSrcBand = (GDALRasterBand *) 
            GDALGetRasterBand( hSrcDS, i+1 );

        poDS->AddBand( poSrcBand->GetRasterDataType(), NULL );

        poBand = (VRTWarpedRasterBand *) poDS->GetRasterBand( i+1 );
        poBand->CopyCommonInfoFrom( poSrcBand );
    }

/* -------------------------------------------------------------------- */
/*      Initialize the warp on the VRTWarpedDataset.                    */
/* -------------------------------------------------------------------- */
    poDS->Initialize( psOptions );
    
    return (GDALDatasetH) poDS;
}

void CDialog3D::OnBnClickedBtnOpen()
{
	// TODO: 在此添加控件通知处理程序代码
	CFileDialog fileDlg(TRUE);
	if(fileDlg.DoModal()!=IDOK)
		return;
	CString strExt	   =fileDlg.GetFileExt();
	CString strPathName=fileDlg.GetPathName();
	if(strExt=="txt")
	{
		ifstream ifs(strPathName,ios_base::in);
		char tmpchr[2048];
		ifs.getline(tmpchr,2048);
		do 
		{
			PNT3D tmpPnt;
			ifs.getline(tmpchr,2048);
			sscanf(tmpchr,"%f,%f,%f",&tmpPnt.DX,&tmpPnt.DY,&tmpPnt.DZ);
			m_vec_PNT3Ds.push_back(tmpPnt);
		} while (!ifs.eof());
		ifs.close();
	}
	if(strExt=="BMP"||strExt=="bmp"||strExt=="JPG"||strExt=="jpg"||strExt=="TIF"||strExt=="tif")
	{
		GDALAllRegister();
		GDALDatasetH hSrcDS=GDALOpen(strPathName,GA_ReadOnly);
		double adfGeoTrans[6];
		double scalex=1,scaley=1;
		GDALGetGeoTransform(hSrcDS,adfGeoTrans);
		int xsize=GDALGetRasterXSize(hSrcDS);
		int ysize=GDALGetRasterYSize(hSrcDS);
		int tmpxsize=xsize,tmpysize=ysize;
		if(xsize>800)
		{
			tmpxsize=800;
			scalex=xsize/tmpxsize;
		}
		if(ysize>600)
		{
			tmpysize=600;
			scaley=ysize/tmpysize;
		}
		float *dataIn=new float[tmpxsize*tmpysize];
		GDALRasterIO(GDALGetRasterBand(hSrcDS,1),GF_Read,0,0,xsize,ysize,dataIn,tmpxsize,tmpysize,GDT_Float32,0,0);
		for(int i=0;i<tmpxsize;i++)
		{
			for (int j=0;j<tmpysize;j++)
			{
				PNT3D tmpPnt;
				tmpPnt.DX=adfGeoTrans[0]+adfGeoTrans[1]*i*scalex+adfGeoTrans[2]*j*scaley;
				tmpPnt.DY=adfGeoTrans[3]+adfGeoTrans[4]*i*scalex+adfGeoTrans[5]*j*scaley;
				//tmpPnt.DX=i;
				//tmpPnt.DY=j;
				tmpPnt.DZ=dataIn[j*tmpxsize+i];
				m_vec_PNT3Ds.push_back(tmpPnt);
			}
		}
		delete[]dataIn;
		GDALClose(hSrcDS);
	}
	if(!m_vec_PNT3Ds.empty())
	{
		m_min_DX=m_max_DX=m_vec_PNT3Ds[0].DX;
		m_min_DY=m_max_DY=m_vec_PNT3Ds[0].DY;
		m_min_DZ=m_max_DZ=m_vec_PNT3Ds[0].DZ;
		for (int i=0;i<m_vec_PNT3Ds.size();i++)
		{
			m_max_DX=max(m_vec_PNT3Ds[i].DX,m_max_DX);
			m_min_DX=min(m_vec_PNT3Ds[i].DX,m_min_DX);
			m_max_DY=max(m_vec_PNT3Ds[i].DY,m_max_DY);
			m_min_DY=min(m_vec_PNT3Ds[i].DY,m_min_DY);
			m_max_DZ=max(m_vec_PNT3Ds[i].DZ,m_max_DZ);
			m_min_DZ=min(m_vec_PNT3Ds[i].DZ,m_min_DZ);
		}
		AfxMessageBox("数据读取成功!\n");
		InvalidateRect(NULL,FALSE);
	}
}

int main( int argc, char ** argv )

{
    GDALDatasetH hSrcDS;
    int         iY, iX, nOutLevel=0, nXSize, nYSize, iArg, nFillDist=0;
    void        *pStream;
    GInt16      *panData;
    const char  *pszFilename = NULL;
    GDALRasterBandH hSrcBand;
    double       adfGeoTransform[6];
    int          bEnableTrim = FALSE;
    GInt16       noDataValue = 0;
    int          bHasNoData;

/* -------------------------------------------------------------------- */
/*      Identify arguments.                                             */
/* -------------------------------------------------------------------- */

    for( iArg = 1; iArg < argc; iArg++ )
    {
        if( EQUAL(argv[iArg],"-trim") )
            bEnableTrim = TRUE;

        else if( EQUAL(argv[iArg],"-fill") )
            nFillDist = atoi(argv[++iArg]);

        else if( EQUAL(argv[iArg],"-level") )
            nOutLevel = atoi(argv[++iArg]);
        else
        {
            if( pszFilename != NULL )
                Usage();
            pszFilename = argv[iArg];
        }
    }

    if( pszFilename == NULL )
        Usage();

/* -------------------------------------------------------------------- */
/*      Open input file.                                                */
/* -------------------------------------------------------------------- */
    GDALAllRegister();
    hSrcDS = GDALOpen( pszFilename, GA_ReadOnly );
    if( hSrcDS == NULL )
        exit(1);

    hSrcBand = GDALGetRasterBand( hSrcDS, 1 );

    noDataValue = (GInt16)GDALGetRasterNoDataValue(hSrcBand, &bHasNoData);

    nXSize = GDALGetRasterXSize( hSrcDS );
    nYSize = GDALGetRasterYSize( hSrcDS );

    GDALGetGeoTransform( hSrcDS, adfGeoTransform );

/* -------------------------------------------------------------------- */
/*      Create output stream.                                           */
/* -------------------------------------------------------------------- */
    pStream = DTEDCreatePtStream( ".", nOutLevel );

    if( pStream == NULL )
        exit( 1 );

/* -------------------------------------------------------------------- */
/*      Process all the profiles.                                       */
/* -------------------------------------------------------------------- */
    panData = (GInt16 *) malloc(sizeof(GInt16) * nXSize);

    for( iY = 0; iY < nYSize; iY++ )
    {
        GDALRasterIO( hSrcBand, GF_Read, 0, iY, nXSize, 1,
                      panData, nXSize, 1, GDT_Int16, 0, 0 );

        if (bHasNoData)
        {
            for( iX = 0; iX < nXSize; iX++ )
            {
                if (panData[iX] == noDataValue)
                    panData[iX] = DTED_NODATA_VALUE;
            }
        }

        for( iX = 0; iX < nXSize; iX++ )
        {
            DTEDWritePt( pStream,
                         adfGeoTransform[0]
                         + adfGeoTransform[1] * (iX + 0.5)
                         + adfGeoTransform[2] * (iY + 0.5),
                         adfGeoTransform[3]
                         + adfGeoTransform[4] * (iX + 0.5)
                         + adfGeoTransform[5] * (iY + 0.5),
                         panData[iX] );
        }
    }

    free( panData );

/* -------------------------------------------------------------------- */
/*      Cleanup.                                                        */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        printf( "Subdatasets:\n" );
        for( i = 0; papszMetadata[i] != NULL; i++ )
        {
            printf( "  %s\n", papszMetadata[i] );
        }
    }

/* -------------------------------------------------------------------- */
/*      Report corners.                                                 */
/* -------------------------------------------------------------------- */
    printf( "Corner Coordinates:\n" );
    GDALInfoReportCorner( hDataset, "Upper Left", 
                          0.0, 0.0 );
    GDALInfoReportCorner( hDataset, "Lower Left", 
                          0.0, GDALGetRasterYSize(hDataset));
    GDALInfoReportCorner( hDataset, "Upper Right", 
                          GDALGetRasterXSize(hDataset), 0.0 );
    GDALInfoReportCorner( hDataset, "Lower Right", 
                          GDALGetRasterXSize(hDataset), 
                          GDALGetRasterYSize(hDataset) );
    GDALInfoReportCorner( hDataset, "Center", 
                          GDALGetRasterXSize(hDataset)/2.0, 
                          GDALGetRasterYSize(hDataset)/2.0 );

/* ==================================================================== */
/*      Loop over bands.                                                */
/* ==================================================================== */
    for( iBand = 0; iBand < GDALGetRasterCount( hDataset ); iBand++ )
    {
        double      dfMin, dfMax, adfCMinMax[2], dfNoData;
        int         bGotMin, bGotMax, bGotNodata;
        int         nBlockXSize, nBlockYSize;

        hBand = GDALGetRasterBand( hDataset, iBand+1 );
        GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
        printf( "Band %d Block=%dx%d Type=%d, ColorInterp=%d\n", iBand+1,
                nBlockXSize, nBlockYSize,
                GDALGetRasterDataType(hBand),
                GDALGetRasterColorInterpretation(hBand) );

        dfMin = GDALGetRasterMinimum( hBand, &bGotMin );
        dfMax = GDALGetRasterMaximum( hBand, &bGotMax );
        printf( "  Min=%.3f/%d, Max=%.3f/%d",  dfMin, bGotMin, dfMax, bGotMax);
        
        if( bComputeMinMax )
        {
            GDALComputeRasterMinMax( hBand, TRUE, adfCMinMax );
            printf( ", Computed Min/Max=%.3f,%.3f", 
                    adfCMinMax[0], adfCMinMax[1] );
        }
        printf( "\n" );

        dfNoData = GDALGetRasterNoDataValue( hBand, &bGotNodata );
        if( bGotNodata )
        {
            printf( "  NoData Value=%g\n", dfNoData );
        }

        if( GDALGetOverviewCount(hBand) > 0 )
        {
            int		iOverview;

            printf( "  Overviews: " );
            for( iOverview = 0; 
                 iOverview < GDALGetOverviewCount(hBand);
                 iOverview++ )

int main( int argc, char ** argv )

{
    const char *pszSrcFilename = NULL;
    int anReqOverviews[1000];
    int nReqOverviewCount = 0;
    bool bMasks = false;

    GDALAllRegister();

    argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
    if( argc < 1 )
        exit( -argc );

/* -------------------------------------------------------------------- */
/*      Process arguments.                                              */
/* -------------------------------------------------------------------- */
    for( int iArg = 1; iArg < argc; iArg++ )
    {
        if( EQUAL(argv[iArg],"-masks") )
        {
            bMasks = true;
        }
        else if( pszSrcFilename == NULL )
        {
            pszSrcFilename = argv[iArg];
        }
        else if( atoi(argv[iArg]) > 0 || EQUAL(argv[iArg],"0") )
        {
            anReqOverviews[nReqOverviewCount++] = atoi(argv[iArg]);
        }
        else
        {
            Usage();
        }
    }

    if( pszSrcFilename == NULL )
        Usage();

/* -------------------------------------------------------------------- */
/*      Open the input file.                                            */
/* -------------------------------------------------------------------- */
    GDALDatasetH hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );

    if( hSrcDS == NULL )
        exit( 1 );

/* ==================================================================== */
/*      Process all bands.                                              */
/* ==================================================================== */
    int iBand;
    int nBandCount = GDALGetRasterCount( hSrcDS );

    for( iBand = 0; iBand < nBandCount; iBand++ )
    {
        GDALRasterBandH hBaseBand = GDALGetRasterBand( hSrcDS, iBand+1 );

/* -------------------------------------------------------------------- */
/*      Process all overviews.                                          */
/* -------------------------------------------------------------------- */
        int iOverview;
        int nOverviewCount = GDALGetOverviewCount( hBaseBand );

        for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
        {
            GDALRasterBandH hSrcOver = GDALGetOverview( hBaseBand, iOverview );

            if (hSrcOver == NULL)
            {
                fprintf(stderr, "skipping overview %d as being null\n", iOverview);
                continue;
            }

/* -------------------------------------------------------------------- */
/*      Is this a requested overview?                                   */
/* -------------------------------------------------------------------- */
            if( nReqOverviewCount > 0 )
            {
                int i;

                for( i = 0; i < nReqOverviewCount; i++ )
                {
                    if( anReqOverviews[i] == iOverview )
                        break;
                }

                if( i == nReqOverviewCount )
                    continue;
            }

/* -------------------------------------------------------------------- */
/*      Create matching output file.                                    */
/* -------------------------------------------------------------------- */
            CPLString osFilename;
            osFilename.Printf( "%s_%d_%d.tif",
                               CPLGetBasename(pszSrcFilename),
                               iBand+1, iOverview );
            DumpBand( hSrcDS, hSrcOver, osFilename );

//.........这里部分代码省略.........

static gint print_handler( void * cb_data, void * scanline_in )

{
    unsigned char *scanline = (unsigned char *) scanline_in;
    static GDALDatasetH  working_ds = NULL;
    static int next_scanline = 0;
    static GDALRasterBandH red_band, green_band, blue_band;
    gint   width;

    if( cb_data == NULL && scanline_in == NULL )
    {
        next_scanline = 0;
        working_ds = NULL;
        return -1;
    }

    if( working_ds != cb_data )
    {
        working_ds = (GDALDatasetH) cb_data;
        next_scanline = 0;

        red_band = GDALGetRasterBand( working_ds, 1 );
        if( GDALGetRasterCount( working_ds ) >= 3 )
        {
            green_band = GDALGetRasterBand( working_ds, 2 );
            blue_band = GDALGetRasterBand( working_ds, 3 );
        }
        else
        {
            green_band = blue_band = NULL;
        }

        if( red_band == NULL )
            return -1;
    }

    width = GDALGetRasterXSize( working_ds );

    if( green_band == NULL )
    {
        GByte	*grey;
        int     i, value;

        grey = g_new( GByte, width );
        
        for( i = 0; i < width; i++ )
        {
            value = *(scanline++);
            value += *(scanline++);
            value += *(scanline++);
            value = (value+1) / 3;
            grey[i] = value;
        }
        
        GDALRasterIO( red_band, GF_Write, 0, next_scanline, width, 1, 
                      grey, width, 1, GDT_Byte, 1, 0 );
        g_free( grey );
    }
    else
    {
        GDALRasterIO( red_band, GF_Write, 0, next_scanline, width, 1, 
                      scanline+0, width, 1, GDT_Byte, 3, 0 );
        GDALRasterIO( green_band, GF_Write, 0, next_scanline, width, 1, 
                      scanline+1, width, 1, GDT_Byte, 3, 0 );
        GDALRasterIO( blue_band, GF_Write, 0, next_scanline, width, 1, 
                      scanline+2, width, 1, GDT_Byte, 3, 0 );
    }
    
    next_scanline++;

    return 0;
}

//.........这里部分代码省略.........
            dfXMax = sEnvelope.MaxX;
            bIsXExtentSet = TRUE;
        }

        if ( !bIsYExtentSet )
        {
            dfYMin = sEnvelope.MinY;
            dfYMax = sEnvelope.MaxY;
            bIsYExtentSet = TRUE;
        }
    }

/* -------------------------------------------------------------------- */
/*      Perform gridding.                                               */
/* -------------------------------------------------------------------- */

    const double    dfDeltaX = ( dfXMax - dfXMin ) / nXSize;
    const double    dfDeltaY = ( dfYMax - dfYMin ) / nYSize;

    if ( !bQuiet )
    {
        printf( "Grid data type is \"%s\"\n", GDALGetDataTypeName(eType) );
        printf( "Grid size = (%lu %lu).\n",
                (unsigned long)nXSize, (unsigned long)nYSize );
        printf( "Corner coordinates = (%f %f)-(%f %f).\n",
                dfXMin - dfDeltaX / 2, dfYMax + dfDeltaY / 2,
                dfXMax + dfDeltaX / 2, dfYMin - dfDeltaY / 2 );
        printf( "Grid cell size = (%f %f).\n", dfDeltaX, dfDeltaY );
        printf( "Source point count = %lu.\n", (unsigned long)adfX.size() );
        PrintAlgorithmAndOptions( eAlgorithm, pOptions );
        printf("\n");
    }

    GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, nBand );

    if (adfX.size() == 0)
    {
        // FIXME: Shoulda' set to nodata value instead
        GDALFillRaster( hBand, 0.0 , 0.0 );
        return CE_None;
    }

    GUInt32 nXOffset, nYOffset;
    int     nBlockXSize, nBlockYSize;
    int     nDataTypeSize = GDALGetDataTypeSize(eType) / 8;

    // Try to grow the work buffer up to 16 MB if it is smaller
    GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
    const GUInt32 nDesiredBufferSize = 16*1024*1024;
    if( (GUInt32)nBlockXSize < nXSize && (GUInt32)nBlockYSize < nYSize &&
        (GUInt32)nBlockXSize < nDesiredBufferSize / (nBlockYSize * nDataTypeSize) )
    {
        int nNewBlockXSize  = nDesiredBufferSize / (nBlockYSize * nDataTypeSize);
        nBlockXSize = (nNewBlockXSize / nBlockXSize) * nBlockXSize;
        if( (GUInt32)nBlockXSize > nXSize )
            nBlockXSize = nXSize;
    }
    else if( (GUInt32)nBlockXSize == nXSize && (GUInt32)nBlockYSize < nYSize &&
             (GUInt32)nBlockYSize < nDesiredBufferSize / (nXSize * nDataTypeSize) )
    {
        int nNewBlockYSize = nDesiredBufferSize / (nXSize * nDataTypeSize);
        nBlockYSize = (nNewBlockYSize / nBlockYSize) * nBlockYSize;
        if( (GUInt32)nBlockYSize > nYSize )
            nBlockYSize = nYSize;
    }
    CPLDebug("GDAL_GRID", "Work buffer: %d * %d", nBlockXSize, nBlockYSize);

bool CUtils::isFloat32DataType(GDALDatasetH hDataset)
{
	return ( GDALGetRasterDataType( GDALGetRasterBand(hDataset, 1) ) == GDT_Float32 );
}

void CUtils::calculateByteGeoTIFFStatistics(GDALDatasetH hDataset, int userBandNumber, byte flNoDataValueAsBackground, byte NoDataValue)
{
	fputs("\nCalculate statistics...\n", stderr);

	GDALRasterBandH hBand =  GDALGetRasterBand(hDataset, 1);
	int cols  = GDALGetRasterBandXSize(hBand);
	int rows  = GDALGetRasterBandYSize(hBand);
	int bands = GDALGetRasterCount(hDataset);
	
	byte * pbuf = NULL;
	pbuf = (byte *)CPLMalloc(sizeof(byte)*cols);
	
	byte min = 0, max = 0, mean = 0;
	double stddev = 0;
    double summ = 0;
    int count = 0;
	
	for(int band=1; band<=bands; band++)
	{
		if(userBandNumber != -1) fprintf(stderr, "Band %d...\n", userBandNumber);
		else fprintf(stderr, "Band %d...\n", band);
		
		hBand =  GDALGetRasterBand(hDataset, band);		
		if(flNoDataValueAsBackground) NoDataValue = getFloatNoDataValueAsBackground(hBand);
		
		min = max = mean = stddev = summ = 0;
		count = 0;
		bool flFirst = true;
		
		int pr = CUtils::progress_ln_ex(stderr, 0, 0, START_PROGRESS);
		for(int i=0; i<rows; i++)
		{
			GDALRasterIO(hBand, GF_Read, 0, i, cols, 1, pbuf, cols, 1, GDT_Byte, 0, 0 );
			for(int j=0; j<cols; j++) if(pbuf[j]!=NoDataValue) 
			{
				if(flFirst)
				{
					mean = pbuf[j];
					min = max = mean;
					flFirst = false;
				}
				else
				{
					mean += pbuf[j];
					if( min > pbuf[j] ) min = pbuf[j];
					if( max < pbuf[j] ) max = pbuf[j];
				}
				count++;
			}
			pr = CUtils::progress_ln_ex(stderr, i, rows, pr);
		}
		CUtils::progress_ln_ex(stderr, 0, 0, END_PROGRESS);
		
		double dmean = 0;
		if(count > 0) dmean = mean / (double)count;
		
		pr = CUtils::progress_ln_ex(stderr, 0, 0, START_PROGRESS);
		for(int i=0; i<rows; i++)
		{
			GDALRasterIO(hBand, GF_Read, 0, i, cols, 1, pbuf, cols, 1, GDT_Byte, 0, 0 );
			for(int j=0; j<cols; j++) if(pbuf[j]!=NoDataValue) summ += ((double)pbuf[j]-dmean)*((double)pbuf[j]-dmean);
						
			pr = CUtils::progress_ln_ex(stderr, i, rows, pr);
		}
		CUtils::progress_ln_ex(stderr, 0, 0, END_PROGRESS);
		
		summ = 0; stddev = 0;
		if((count-1)>0)
		{
			summ /= (double)(count-1);
			if(summ!=0) stddev = sqrt(summ);
		}
			
		GDALSetRasterStatistics(hBand, min, max, mean, stddev);
		GDALSetRasterNoDataValue(hBand, NoDataValue);
	}
	
	CPLFree(pbuf);
}

//.........这里部分代码省略.........
                if (eErr != CE_None)
                {
                    break;
                }
            }

            GDALDatasetH hMemDS = GDALCreate(hMemDriver, "MEM:::",
                                              nReqXSize, nReqYSize, 0, 
                                              eDataType, NULL);
            if (hMemDS == NULL)
            {
                eErr = CE_Failure;
                break;
            }
            
            int iBand;
            for(iBand = 0; iBand < nBands; iBand ++)
            {
                char** papszOptions = NULL;
                char szTmp[64];
                memset(szTmp, 0, sizeof(szTmp));
                CPLPrintPointer(szTmp,
                                pabyMEMDSBuffer + iBand * nDataTypeSize *
                                nReqXSize * nReqYSize, sizeof(szTmp));
                papszOptions = CSLSetNameValue(papszOptions, "DATAPOINTER", szTmp);
                GDALAddBand(hMemDS, eDataType, papszOptions);
                CSLDestroy(papszOptions);
            }

            if( hPrevOvrMemDS != NULL )
            {
                for(iBand = 0; iBand < nBands; iBand ++)
                {
                    GDALRasterBandH hDstOvrBand = GDALGetRasterBand(hMemDS, iBand+1);

                    eErr = GDALRegenerateOverviews( GDALGetRasterBand(hPrevOvrMemDS, iBand+1),
                                                    1, &hDstOvrBand,
                                                    pszResampling,
                                                    NULL, NULL );
                    if( eErr != CE_None )
                        break;
                }

                GDALClose(hPrevOvrMemDS);
            }

            GDALDatasetH hOutDS = GDALCreateCopy(hTileDriver,
                                        osTempFileName.c_str(), hMemDS, FALSE,
                                        papszTileDriverOptions, NULL, NULL);

            GDALClose(hMemDS);
            if (hOutDS)
                GDALClose(hOutDS);
            else
            {
                eErr = CE_Failure;
                break;
            }

/* -------------------------------------------------------------------- */
/*      Insert new entry into raster table                              */
/* -------------------------------------------------------------------- */

            vsi_l_offset nDataLength;
            GByte *pabyData = VSIGetMemFileBuffer( osTempFileName.c_str(),
                                                   &nDataLength, FALSE);

void convert(const input_arguments& args)
{
// determine the subdataset pattern name using the first available file
  if(args.verbose)
    std::cout << "\tlooking for the names of subdatasets (or variables)... " << std::flush;
  
  std::vector<std::string> band_names = modis2scidb::extract_subdatasets_pattern_names(args.source_file_name);

  if(args.verbose)
    std::cout << "OK!" << std::endl;
  
// check if band numbers are in the valid range
  if(args.verbose)
    std::cout << "\tchecking the range for choosed band numbers... " << std::flush;
  
  std::size_t num_bands = args.bands.size();
  
  for(std::size_t i = 0; i != num_bands; ++i)
    if(args.bands[i] >= band_names.size())
      throw modis2scidb::invalid_arg_value() << modis2scidb::error_description("band number is invalid!");
  
  if(args.verbose)
    std::cout << "OK!" << std::endl;

// let's buffer each subdataset
  if(args.verbose)
    std::cout << "\tbuffering data... " << std::flush;

  std::vector<boost::shared_array<unsigned char> > data_buffers;
  std::vector<unsigned char*> aux_data_buffers;
  std::vector<std::size_t> band_datatype_size;
  
  int64_t ncols = 0;
  int64_t nrows = 0;

  for(std::size_t i = 0; i != num_bands; ++i)
  {
    if(args.verbose)
      std::cout << "\n\t\tband #" << args.bands[i] << "... " << std::flush;

    boost::format subdataset(band_names[args.bands[i]]);

    subdataset.bind_arg(1, args.source_file_name);
    
    GDALDatasetH dataset = GDALOpen(subdataset.str().c_str(), GA_ReadOnly);
    
    if(dataset == 0)
    {
      boost::format err_msg("could not open subdataset: '%1%', for input hdf file: '%2%'!");
      throw modis2scidb::gdal_error() << modis2scidb::error_description((err_msg % subdataset.str() % args.source_file_name).str());
    }
    
    GDALRasterBandH band = GDALGetRasterBand(dataset, 1);
    
    if(band == 0)
    {
      GDALClose(dataset);
      boost::format err_msg("could not access band: %1%!");
      throw modis2scidb::gdal_error() << modis2scidb::error_description((err_msg % args.bands[i]).str());
    }
    
    if(i == 0)
    {
      ncols = GDALGetRasterBandXSize(band);
      nrows = GDALGetRasterBandYSize(band);
    }
    else
    {
      if((GDALGetRasterBandXSize(band) != ncols) ||
         (GDALGetRasterBandYSize(band) != nrows))
      {
        GDALClose(dataset);
        throw modis2scidb::gdal_error() << modis2scidb::error_description("selected bands must have the same dimension (rows and cols)!");
      }
    }
    
    GDALDataType pixel_type = GDALGetRasterDataType(band);
    
    std::size_t pixel_size = modis2scidb::num_bytes(pixel_type);
    
    band_datatype_size.push_back(pixel_size);
    
    boost::shared_array<unsigned char> buffer(new unsigned char[ncols * nrows * pixel_size]);
    
    data_buffers.push_back(buffer);
    
    aux_data_buffers.push_back(buffer.get());
    
    CPLErr result = GDALRasterIO(band, GF_Read, 0, 0, static_cast<int>(ncols), static_cast<int>(nrows), buffer.get(), static_cast<int>(ncols), static_cast<int>(nrows), pixel_type, 0, 0);

    if(result == CE_Failure)
    {
      GDALClose(dataset);
      boost::format err_msg("could not read subdataset: '%1%', for input hdf file: '%2%'!");
      throw modis2scidb::gdal_error() << modis2scidb::error_description((err_msg % subdataset.str() % args.source_file_name).str());
    }

    GDALClose(dataset);

    if(args.verbose)
//.........这里部分代码省略.........

/**
 * @param theBandNumber the number of the band for which you want a color table
 * @param theList a pointer the object that will hold the color table
 * @return true of a color table was able to be read, false otherwise
 */
QList<QgsColorRampShader::ColorRampItem> QgsGdalProviderBase::colorTable( GDALDatasetH gdalDataset, int theBandNumber )const
{
  QgsDebugMsg( "entered." );
  QList<QgsColorRampShader::ColorRampItem> ct;

  //Invalid band number, segfault prevention
  if ( 0 >= theBandNumber )
  {
    QgsDebugMsg( "Invalid parameter" );
    return ct;
  }

  GDALRasterBandH myGdalBand = GDALGetRasterBand( gdalDataset, theBandNumber );
  GDALColorTableH myGdalColorTable = GDALGetRasterColorTable( myGdalBand );

  if ( myGdalColorTable )
  {
    QgsDebugMsg( "Color table found" );

    // load category labels
    char ** categoryNames = GDALGetRasterCategoryNames( myGdalBand );
    QVector<QString> labels;
    if ( categoryNames )
    {
      int i = 0;
      while ( categoryNames[i] )
      {
        labels.append( QString( categoryNames[i] ) );
        i++;
      }
    }

    int myEntryCount = GDALGetColorEntryCount( myGdalColorTable );
    GDALColorInterp myColorInterpretation =  GDALGetRasterColorInterpretation( myGdalBand );
    QgsDebugMsg( "Color Interpretation: " + QString::number(( int )myColorInterpretation ) );
    GDALPaletteInterp myPaletteInterpretation  = GDALGetPaletteInterpretation( myGdalColorTable );
    QgsDebugMsg( "Palette Interpretation: " + QString::number(( int )myPaletteInterpretation ) );

    const GDALColorEntry* myColorEntry = 0;
    for ( int myIterator = 0; myIterator < myEntryCount; myIterator++ )
    {
      myColorEntry = GDALGetColorEntry( myGdalColorTable, myIterator );

      if ( !myColorEntry )
      {
        continue;
      }
      else
      {
        QString label = labels.value( myIterator );
        if ( label.isEmpty() )
        {
          label = QString::number( myIterator );
        }
        //Branch on the color interpretation type
        if ( myColorInterpretation == GCI_GrayIndex )
        {
          QgsColorRampShader::ColorRampItem myColorRampItem;
          myColorRampItem.value = ( double )myIterator;
          myColorRampItem.label = label;
          myColorRampItem.color = QColor::fromRgb( myColorEntry->c1, myColorEntry->c1, myColorEntry->c1, myColorEntry->c4 );
          ct.append( myColorRampItem );
        }
        else if ( myColorInterpretation == GCI_PaletteIndex )
        {
          QgsColorRampShader::ColorRampItem myColorRampItem;
          myColorRampItem.value = ( double )myIterator;
          myColorRampItem.label = label;
          //Branch on palette interpretation
          if ( myPaletteInterpretation  == GPI_RGB )
          {
            myColorRampItem.color = QColor::fromRgb( myColorEntry->c1, myColorEntry->c2, myColorEntry->c3, myColorEntry->c4 );
          }
          else if ( myPaletteInterpretation  == GPI_CMYK )
          {
            myColorRampItem.color = QColor::fromCmyk( myColorEntry->c1, myColorEntry->c2, myColorEntry->c3, myColorEntry->c4 );
          }
          else if ( myPaletteInterpretation  == GPI_HLS )
          {
            myColorRampItem.color = QColor::fromHsv( myColorEntry->c1, myColorEntry->c3, myColorEntry->c2, myColorEntry->c4 );
          }
          else
          {
            myColorRampItem.color = QColor::fromRgb( myColorEntry->c1, myColorEntry->c1, myColorEntry->c1, myColorEntry->c4 );
          }
          ct.append( myColorRampItem );
        }
        else
        {
          QgsDebugMsg( "Color interpretation type not supported yet" );
          return ct;
        }
      }
    }
  }
//.........这里部分代码省略.........

CPLErr CPL_STDCALL
GDALFillNodata( GDALRasterBandH hTargetBand,
                GDALRasterBandH hMaskBand,
                double dfMaxSearchDist,
                CPL_UNUSED int bDeprecatedOption,
                int nSmoothingIterations,
                char **papszOptions,
                GDALProgressFunc pfnProgress,
                void * pProgressArg )

{
    VALIDATE_POINTER1( hTargetBand, "GDALFillNodata", CE_Failure );

    const int nXSize = GDALGetRasterBandXSize(hTargetBand);
    const int nYSize = GDALGetRasterBandYSize(hTargetBand);

    if( dfMaxSearchDist == 0.0 )
        dfMaxSearchDist = std::max(nXSize, nYSize) + 1;

    const int nMaxSearchDist = static_cast<int>(floor(dfMaxSearchDist));

    // Special "x" pixel values identifying pixels as special.
    GDALDataType eType = GDT_UInt16;
    GUInt32 nNoDataVal = 65535;

    if( nXSize > 65533 || nYSize > 65533 )
    {
        eType = GDT_UInt32;
        nNoDataVal = 4000002;
    }

    if( hMaskBand == nullptr )
        hMaskBand = GDALGetMaskBand( hTargetBand );

    // If there are smoothing iterations, reserve 10% of the progress for them.
    const double dfProgressRatio = nSmoothingIterations > 0 ? 0.9 : 1.0;

    const char* pszNoData = CSLFetchNameValue(papszOptions, "NODATA");
    bool bHasNoData = false;
    float fNoData = 0.0f;
    if( pszNoData )
    {
        bHasNoData = true;
        fNoData = static_cast<float>(CPLAtof(pszNoData));
    }

/* -------------------------------------------------------------------- */
/*      Initialize progress counter.                                    */
/* -------------------------------------------------------------------- */
    if( pfnProgress == nullptr )
        pfnProgress = GDALDummyProgress;

    if( !pfnProgress( 0.0, "Filling...", pProgressArg ) )
    {
        CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
        return CE_Failure;
    }

/* -------------------------------------------------------------------- */
/*      Determine format driver for temp work files.                    */
/* -------------------------------------------------------------------- */
    CPLString osTmpFileDriver = CSLFetchNameValueDef(
            papszOptions, "TEMP_FILE_DRIVER", "GTiff");
    GDALDriverH hDriver = GDALGetDriverByName(osTmpFileDriver.c_str());

    if( hDriver == nullptr )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Given driver is not registered");
        return CE_Failure;
    }

    if( GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) == nullptr )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Given driver is incapable of creating temp work files");
        return CE_Failure;
    }

    char **papszWorkFileOptions = nullptr;
    if( osTmpFileDriver == "GTiff" )
    {
        papszWorkFileOptions = CSLSetNameValue(
                papszWorkFileOptions, "COMPRESS", "LZW");
        papszWorkFileOptions = CSLSetNameValue(
                papszWorkFileOptions, "BIGTIFF", "IF_SAFER");
    }

/* -------------------------------------------------------------------- */
/*      Create a work file to hold the Y "last value" indices.          */
/* -------------------------------------------------------------------- */
    const CPLString osTmpFile = CPLGenerateTempFilename("");
    const CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";

    GDALDatasetH hYDS =
        GDALCreate( hDriver, osYTmpFile, nXSize, nYSize, 1,
                    eType, papszWorkFileOptions );

    if( hYDS == nullptr )
    {
//.........这里部分代码省略.........

static void DumpBand( GDALDatasetH hBaseDS, GDALRasterBandH hSrcOver,
                      const char *pszName )

{
/* -------------------------------------------------------------------- */
/*      Get base ds info.                                               */
/* -------------------------------------------------------------------- */
    double adfGeoTransform[6];
    bool bHaveGT = GDALGetGeoTransform( hBaseDS, adfGeoTransform ) == CE_None;

    int nOrigXSize = GDALGetRasterXSize( hBaseDS );
    int nOrigYSize = GDALGetRasterYSize( hBaseDS );

/* -------------------------------------------------------------------- */
/*      Create matching output file.                                    */
/* ---------------------------------