C++ GDALRasterBand::GetMinimum方法代码示例

bool getRawValuesFromFile(string fname,vector<vector<float>>& vecs)
{
   
   //vector<float> temp = vector<float>()
   GDALDataset *poDataset;
   GDALAllRegister();
   poDataset= (GDALDataset*) GDALOpen(fname.c_str(),GA_ReadOnly);
   if(poDataset == NULL)
   {
      cout << "OUCH!" << endl;
      return false;
   }
   cout << "Data size: " << GDALGetRasterXSize(poDataset) << " " << GDALGetRasterYSize(poDataset) << endl;

   GDALRasterBand  *poBand;
   int             nBlockXSize, nBlockYSize;
   int             bGotMin, bGotMax;
   double          adfMinMax[2];
        
   poBand = poDataset->GetRasterBand( 1 );
   poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
   printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
            nBlockXSize, nBlockYSize,
            GDALGetDataTypeName(poBand->GetRasterDataType()),
            GDALGetColorInterpretationName(
            poBand->GetColorInterpretation()) );

   adfMinMax[0] = poBand->GetMinimum( &bGotMin );
   adfMinMax[1] = poBand->GetMaximum( &bGotMax );
   if( ! (bGotMin && bGotMax) )
      GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
   int width = poBand->GetXSize(); 
   int height = poBand->GetYSize();
   int bands = poDataset->GetRasterCount();
   float *pafScanline;
   std::cout << "Before allocation" << adfMinMax[0] << " " << adfMinMax[1] << endl;
   int dsize = 256;
   pafScanline = (float *) CPLMalloc(sizeof(float)*width*height);
   vector<vector<float>> out = vector<vector<float>>(height,vector<float> (width,0));
   poBand->RasterIO(GF_Read,0,0,width,height,pafScanline,width,height,GDT_Float32,0,0);
   cout << "After allocation" << endl;
   for(int i = 0; i < height; i++)
   {
    for(int j = 0; j < width; j++)
    {
      //cout << i << j << endl << pafS;
      out[i][j] = pafScanline[i*width+j];
    }
   }
   CPLFree(pafScanline);
   //for(auto i : out)
   //for(auto j : i)
   //		cout << j << endl;
  cout << "After allocation" << endl;
  vecs = out;
   return true;
}

int main()
{
    GDALDataset  *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset *) GDALOpen( "GE01.tif", GA_ReadOnly );
    printf("Working! \n");

    if( poDataset != NULL ){
    	//Get Dataset Information
        double adfGeoTransform[6];

		printf( "Driver: %s/%s\n", poDataset->GetDriver()->GetDescription(), poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );

		printf( "Size is %dx%dx%d\n", poDataset->GetRasterXSize(), poDataset->GetRasterYSize(), poDataset->GetRasterCount() );

		if( poDataset->GetProjectionRef()  != NULL )
    		printf( "Projection is `%s'\n", poDataset->GetProjectionRef() );

		if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None ){
    		printf( "Origin = (%.6f,%.6f)\n",
            adfGeoTransform[0], adfGeoTransform[3] );
    		printf( "Pixel Size = (%.6f,%.6f)\n",
            adfGeoTransform[1], adfGeoTransform[5] );
		}

		//Fetch Raster Band
		GDALRasterBand  *poBand;
		int             nBlockXSize, nBlockYSize;
		int             bGotMin, bGotMax;
		double          adfMinMax[2];
		poBand = poDataset->GetRasterBand( 1 );
		poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
		printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
		        nBlockXSize, nBlockYSize,
		        GDALGetDataTypeName(poBand->GetRasterDataType()),
		        GDALGetColorInterpretationName(
		            poBand->GetColorInterpretation()) );
		adfMinMax[0] = poBand->GetMinimum( &bGotMin );
		adfMinMax[1] = poBand->GetMaximum( &bGotMax );
		if( ! (bGotMin && bGotMax) )
		    GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
		printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
		if( poBand->GetOverviewCount() > 0 )
		    printf( "Band has %d overviews.\n", poBand->GetOverviewCount() );
		if( poBand->GetColorTable() != NULL )
		    printf( "Band has a color table with %d entries.\n", 
		             poBand->GetColorTable()->GetColorEntryCount() );

		//Close Dataset
		GDALClose(poDataset);
		//Exit      
		return 0;
	}


}

void HypSpecImage::bounds(int band, float *min, float *max)
{
    GDALRasterBand *raster = m_data->GetRasterBand(band);
    int bGotMin, bGotMax;
    *min = raster->GetMinimum(&bGotMin);
    *max = raster->GetMaximum(&bGotMax);
    if(!(bGotMin && bGotMax)) {
        double adfMinMax[2];
        GDALComputeRasterMinMax((GDALRasterBandH) raster, TRUE, adfMinMax);
        *min = adfMinMax[0];
        *max = adfMinMax[1];
    }
}

bool Slic::_InitData()
{
    // Check the params
    if (_poSrcDS == NULL|| _poDstDS == NULL)
    {
        std::cerr<<"Input image or output image invalid !"<<std::endl;
        return false;
    }

    if (_regionSize<0 || _regularizer<0)
    {
        std::cerr<<"Parameter regionSize and regularizer must bigger than 0!"<<std::endl;
        return false;
    }

    // Init some vars
    _M = static_cast<int> (static_cast<double>(_width)  / _regionSize + 0.5);
    _N = static_cast<int> (static_cast<double>(_height) / _regionSize + 0.5);

    // Init normalization params
    for (int k=0;k<_bandCount;++k)
    {
        GDALRasterBand* poBand =  _poSrcDS->GetRasterBand(k+1);
        int     bGotMin, bGotMax;
        double  adfMinMax[2];
        adfMinMax[0] = poBand->GetMinimum( &bGotMin );
        adfMinMax[1] = poBand->GetMaximum( &bGotMax );
        if( ! (bGotMin && bGotMax) )
            GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
        std::vector<double>  adfNormalizationParam(2);
        adfNormalizationParam[0] = 1./(adfMinMax[1]-adfMinMax[0]);
        adfNormalizationParam[1] = adfMinMax[0]/(adfMinMax[1]-adfMinMax[0]);
        _normalizationParam.push_back(adfNormalizationParam);
    }

    // Init centers
    for (int i=_regionSize/2;i<_height ;i+=_regionSize)
    {
        for (int j=_regionSize/2;j<_width ;j += _regionSize)
        {
            FeatureVector featureVec;
            _InitCenterFeature(i,j,featureVec);
            _centerVector.push_back(featureVec);
        }
    }

    return true;
}

void Create8BitImage(const char *srcfile,const char *dstfile)
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");
	GDALDataset *pDataset=(GDALDataset *) GDALOpen( srcfile, GA_ReadOnly );
	int bandNum=pDataset->GetRasterCount();

	GDALDriver *pDriver=GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset *dstDataset=pDriver->Create(dstfile,800,800,3,GDT_Byte,NULL);
	GDALRasterBand *pBand;
	GDALRasterBand *dstBand;
	//写入光栅数据
	ushort *sbuf= new ushort[800*800];
	uchar *cbuf=new uchar[800*800];
	for (int i=bandNum,j=1;i>=2;i--,j++)
	{
		pBand=pDataset->GetRasterBand(i);
		pBand->RasterIO(GF_Read,0,0,800,800,sbuf,800,800,GDT_UInt16,0,0);
		int bGotMin, bGotMax;
		double adfMinMax[2];
		adfMinMax[0] = pBand->GetMinimum( &bGotMin );
		adfMinMax[1] = pBand->GetMaximum( &bGotMax );
		if( ! (bGotMin && bGotMax) )
			GDALComputeRasterMinMax((GDALRasterBandH)pBand, TRUE, adfMinMax);
		MinMaxStretch(sbuf,cbuf,800,800,adfMinMax[0],adfMinMax[1]);
		/*double min,max;
		HistogramAccumlateMinMax16S(sbuf,800,800,&min,&max);
		MinMaxStretchNew(sbuf,cbuf,800,800,min,max);*/
		dstBand=dstDataset->GetRasterBand(j);
		dstBand->RasterIO(GF_Write,0,0,800,800,cbuf,800,800,GDT_Byte,0,0);
	}
	delete []cbuf;
	delete []sbuf;
	GDALClose(pDataset);
	GDALClose(dstDataset);
}

//.........这里部分代码省略.........
            GDALClose((GDALDatasetH)poDataset);
            return false;
        }
    }

    qDebug( "Driver: %s/%s\n",
            poDataset->GetDriver()->GetDescription(),
            poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );

    qDebug( "Size is %dx%dx%d\n",
            poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
            poDataset->GetRasterCount() );

    GdalAdapter::ImgType theType = GdalAdapter::Unknown;
    int bandCount = poDataset->GetRasterCount();
    int ixA = -1;
    int ixR, ixG, ixB;
    int ixH, ixS, ixL;
    int ixC, ixM, ixY, ixK;
    int ixYuvY, ixYuvU, ixYuvV;
    double adfMinMax[2];
    double UnknownUnit;
    GDALColorTable* colTable = NULL;
    for (int i=0; i<bandCount; ++i) {
        GDALRasterBand  *poBand = poDataset->GetRasterBand( i+1 );
        GDALColorInterp bandtype = poBand->GetColorInterpretation();
        qDebug() << "Band " << i+1 << " Color: " <<  GDALGetColorInterpretationName(poBand->GetColorInterpretation());

        switch (bandtype)
        {
        case GCI_Undefined:
            theType = GdalAdapter::Unknown;
            int             bGotMin, bGotMax;
            adfMinMax[0] = poBand->GetMinimum( &bGotMin );
            adfMinMax[1] = poBand->GetMaximum( &bGotMax );
            if( ! (bGotMin && bGotMax) )
                GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
            UnknownUnit = (adfMinMax[1] - adfMinMax[0]) / 256;
            break;
        case GCI_GrayIndex:
            theType = GdalAdapter::GrayScale;
            break;
        case GCI_RedBand:
            theType = GdalAdapter::Rgb;
            ixR = i;
            break;
        case GCI_GreenBand:
            theType = GdalAdapter::Rgb;
            ixG = i;
            break;
        case GCI_BlueBand :
            theType = GdalAdapter::Rgb;
            ixB = i;
            break;
        case GCI_HueBand:
            theType = GdalAdapter::Hsl;
            ixH = i;
            break;
        case GCI_SaturationBand:
            theType = GdalAdapter::Hsl;
            ixS = i;
            break;
        case GCI_LightnessBand:
            theType = GdalAdapter::Hsl;
            ixL = i;
            break;

std::tuple<boost::shared_ptr<Map_Matrix<DataFormat> >, std::string, GeoTransform>  read_in_map(fs::path file_path, GDALDataType data_type, const bool doCategorise) throw(std::runtime_error)
{

    std::string projection;
    GeoTransform transformation;
    GDALDriver driver;
    
	//Check that the file name is valid
	if (!(fs::is_regular_file(file_path)))
	{
		throw std::runtime_error("Input file is not a regular file");
	}	

	// Get GDAL to open the file - code is based on the tutorial at http://www.gdal.org/gdal_tutorial.html
	GDALDataset *poDataset;
	GDALAllRegister(); //This registers all availble raster file formats for use with this utility. How neat is that. We can input any GDAL supported rater file format.

	//Open the Raster by calling GDALOpen. http://www.gdal.org/gdal_8h.html#a6836f0f810396c5e45622c8ef94624d4
	//char pszfilename[] = file_path.c_str(); //Set this to the file name, as GDALOpen requires the standard C char pointer as function parameter.
	poDataset = (GDALDataset *) GDALOpen (file_path.string().c_str(), GA_ReadOnly);
	if (poDataset == NULL)
	{
		throw std::runtime_error("Unable to open file");
	}
	
//	Print some general information about the raster
	double        adfGeoTransform[6]; //An array of doubles that will be used to save information about the raster - where the origin is, what the raster pizel size is.
    
    
    printf( "Driver: %s/%s\n",
            poDataset->GetDriver()->GetDescription(),
            poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );

    printf( "Size is %dx%dx%d\n", 
            poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
            poDataset->GetRasterCount() );

    if( poDataset->GetProjectionRef()  != NULL )
    {
        printf( "Projection is `%s'\n", poDataset->GetProjectionRef() );
        projection = poDataset->GetProjectionRef();
    }

    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
        printf( "Origin = (%.6f,%.6f)\n",
                adfGeoTransform[0], adfGeoTransform[3] );

        printf( "Pixel Size = (%.6f,%.6f)\n",
                adfGeoTransform[1], adfGeoTransform[5] );
        
        transformation.x_origin = adfGeoTransform[0];
        transformation.pixel_width = adfGeoTransform[1];
        transformation.x_line_space = adfGeoTransform[2];
        transformation.y_origin = adfGeoTransform[3];
        transformation.pixel_height = adfGeoTransform[4];
        transformation.y_line_space = adfGeoTransform[5];
        
    }


	/// Some raster file formats allow many layers of data (called a 'band', with each having the same pixel size and origin location and spatial extent). We will get the data for the first layer into a Boost Array.
	//Get the data from the first band, 
	// TODO implement method with input to specify what band.
	    GDALRasterBand  *poBand;
        int             nBlockXSize, nBlockYSize;
        int             bGotMin, bGotMax;
        double          adfMinMax[2];
        
        poBand = poDataset->GetRasterBand( 1 );
        poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
        printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
                nBlockXSize, nBlockYSize,
                GDALGetDataTypeName(poBand->GetRasterDataType()),
                GDALGetColorInterpretationName(
                    poBand->GetColorInterpretation()) );

        adfMinMax[0] = poBand->GetMinimum( &bGotMin );
        adfMinMax[1] = poBand->GetMaximum( &bGotMax );
        if( ! (bGotMin && bGotMax) )
            GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);

        printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
        
        if( poBand->GetOverviewCount() > 0 )
            printf( "Band has %d overviews.\n", poBand->GetOverviewCount() );

        if( poBand->GetColorTable() != NULL )
            printf( "Band has a color table with %d entries.\n", 
                     poBand->GetColorTable()->GetColorEntryCount() );

		DataFormat * pafScanline;
        int   nXSize = poBand->GetXSize();
		int   nYSize = poBand->GetYSize();
		boost::shared_ptr<Map_Matrix<DataFormat> > in_map(new Map_Matrix<DataFormat>(nYSize, nXSize));
		//get a c array of this size and read into this.
		
		
		//pafScanline = new DataFormat[nXSize];
		//for (int i = 0; i < nYSize; i++) //rows
//.........这里部分代码省略.........

bool GDALImageFileType::read(      Image       *OSG_GDAL_ARG(pImage), 
                             const Char8       *OSG_GDAL_ARG(fileName)) 
{
#ifdef OSG_WITH_GDAL
    bool returnValue = false;

    GDALDataset *pDataset;

    pDataset = static_cast<GDALDataset *>(GDALOpen(fileName, GA_ReadOnly));

    if(pDataset != NULL)
    {
        GeoReferenceAttachmentUnrecPtr pGeoRef = 
            GeoReferenceAttachment::create();

        pImage->addAttachment(pGeoRef);

        double        adfGeoTransform[6];
        
        if(pDataset->GetGeoTransform(adfGeoTransform) == CE_None)
        {
            pGeoRef->editOrigin().setValues(adfGeoTransform[0], 
                                            adfGeoTransform[3]);

            pGeoRef->editPixelSize().setValues(adfGeoTransform[1], 
                                               adfGeoTransform[5]);

            if(GDALGetProjectionRef(pDataset) != NULL)
            {
                OGRSpatialReferenceH  hSRS;

                Char8 *szProjection = 
                    const_cast<char *>(GDALGetProjectionRef(pDataset));
        
                hSRS = OSRNewSpatialReference(NULL);

                if(OSRImportFromWkt(hSRS, &szProjection) == CE_None)
                {
                    pGeoRef->editEllipsoidAxis().setValues(
                        OSRGetSemiMajor(hSRS, NULL),
                        OSRGetSemiMinor(hSRS, NULL));

                    const Char8 *szDatum = OSRGetAttrValue(hSRS, "DATUM", 0);

                    if(szDatum != NULL && 0 == strcmp(szDatum, "WGS_1984"))
                    {
                        pGeoRef->editDatum() = 
                            GeoReferenceAttachment::WGS84;
                    }
                    else
                    {
                        fprintf(stderr, "Unknow datum %s\n",
                                szDatum);

                        pGeoRef->editDatum() = 
                            GeoReferenceAttachment::UnknownDatum;
                    }
                }            

                OSRDestroySpatialReference(hSRS);
            }
        }

        GDALRasterBand *pBand;
        int             nBlockXSize, nBlockYSize;
        int             bGotMin, bGotMax;
        double          adfMinMax[2];
        
        pBand = pDataset->GetRasterBand( 1 );
        pBand->GetBlockSize( &nBlockXSize, &nBlockYSize );

        adfMinMax[0] = pBand->GetMinimum( &bGotMin );
        adfMinMax[1] = pBand->GetMaximum( &bGotMax );

        if(!(bGotMin && bGotMax))
            GDALComputeRasterMinMax(GDALRasterBandH(pBand), TRUE, adfMinMax);

        pBand = pDataset->GetRasterBand(1);

        if(pBand != NULL)
        {
            Image::PixelFormat ePF = Image::OSG_INVALID_PF;

            switch(pDataset->GetRasterCount())
            {
                case 1:
                    ePF = Image::OSG_L_PF;
                    break;
                case 2:
                    ePF = Image::OSG_LA_PF;
                    break;
                case 3:
                    ePF = Image::OSG_RGB_PF;
                    break;
                case 4:
                    ePF = Image::OSG_RGBA_PF;
                    break;
            }

            Image::Type eDT = Image::OSG_INVALID_IMAGEDATATYPE;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  size_t nBlockXSize, nBlockYSize, nBands;
  int bGotMin, bGotMax;
  double adfMinMax[2];
  float *pafScanline;

  while( true ) {
    pfs::FrameFile ff = it.getNextFrameFile();
    if( ff.fh == NULL ) break; // No more frames
    it.closeFrameFile( ff );

    VERBOSE_STR << "reading file '" << ff.fileName << "'" << std::endl;
    if( !( poDataset = (GDALDataset *) GDALOpen( ff.fileName, GA_ReadOnly ) ) ) {
      std::cerr << "input does not seem to be in a format supported by GDAL" << std::endl;
      throw QuietException();
    }
    VERBOSE_STR << "GDAL driver: " << poDataset->GetDriver()->GetDescription()
	<< " / " << poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) << std::endl;
    nBlockXSize = poDataset->GetRasterXSize();
    nBlockYSize = poDataset->GetRasterYSize();
    nBands = poDataset->GetRasterCount();
    VERBOSE_STR << "Data size " << nBlockXSize << "x" << nBlockYSize << "x" << nBands << std::endl;
    if( poDataset->GetProjectionRef() ) {
      VERBOSE_STR << "Projection " << poDataset->GetProjectionRef() << std::endl;
    }
    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None ) {
      VERBOSE_STR << "Origin = (" << adfGeoTransform[0] << ", " << adfGeoTransform[3] << ")" << std::endl;
      VERBOSE_STR << "Pixel Size = (" << adfGeoTransform[1] << ", " << adfGeoTransform[5] << ")" << std::endl;
    }

    if( nBlockXSize==0  || nBlockYSize==0 
	    || ( SIZE_MAX / nBlockYSize < nBlockXSize ) 
	    || ( SIZE_MAX / (nBlockXSize * nBlockYSize ) < 4 ) ) {
      std::cerr << "input data has invalid size" << std::endl;
      throw QuietException();
    }
    if( !(pafScanline = (float *) CPLMalloc( sizeof(float) * nBlockXSize ) ) ) {
      std::cerr << "not enough memory" << std::endl;
      throw QuietException();
    }

    pfs::Frame *frame = pfsio.createFrame( nBlockXSize, nBlockYSize );
    pfs::Channel *C[nBands];
    char channel_name[32];

    frame->getTags()->setString( "X-GDAL_DRIVER_SHORTNAME", poDataset->GetDriver()->GetDescription() );
    frame->getTags()->setString( "X-GDAL_DRIVER_LONGNAME", poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );
    frame->getTags()->setString( "X-PROJECTION", poDataset->GetProjectionRef() );
    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None ) {
      frame->getTags()->setString( "X-ORIGIN_X", stringify(adfGeoTransform[0]).c_str() );
      frame->getTags()->setString( "X-ORIGIN_Y", stringify(adfGeoTransform[3]).c_str() );
      frame->getTags()->setString( "X-PIXEL_WIDTH", stringify(adfGeoTransform[1]).c_str() );
      frame->getTags()->setString( "X-PIXEL_HEIGHT", stringify(adfGeoTransform[5]).c_str() );
    }

    for ( size_t band = 1; band <= nBands; band++) {
      size_t nBandXSize, nBandYSize;
      VERBOSE_STR << "Band " << band << ": " << std::endl;
      snprintf( channel_name, 32, "X-GDAL%zu", band );
      C[band - 1] = frame->createChannel( channel_name );
      poBand = poDataset->GetRasterBand( band );
      nBandXSize = poBand->GetXSize();
      nBandYSize = poBand->GetYSize();
      VERBOSE_STR << "    " << nBandXSize << "x" << nBandYSize << std::endl;
      if( nBandXSize != (int)nBlockXSize || nBandYSize != (int)nBlockYSize ) {
        std::cerr << "data in band " << band << " has different size" << std::endl;
        throw QuietException();
      }
      VERBOSE_STR << "    Type " << GDALGetDataTypeName( poBand->GetRasterDataType() ) << ", "
	  << "Color Interpretation " << GDALGetColorInterpretationName( poBand->GetColorInterpretation() ) << std::endl;
      adfMinMax[0] = poBand->GetMinimum( &bGotMin );
      adfMinMax[1] = poBand->GetMaximum( &bGotMax );
      if( ! (bGotMin && bGotMax) ) {
      	  GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
      }
      VERBOSE_STR << "    Min " << adfMinMax[0] << ", Max " << adfMinMax[1] << std::endl;
      C[band - 1]->getTags()->setString( "X-TYPE", 
	      GDALGetDataTypeName( poBand->GetRasterDataType() ) );
      C[band - 1]->getTags()->setString( "X-COLOR_INTERPRETATION", 
	      GDALGetColorInterpretationName( poBand->GetColorInterpretation() ) );
      C[band - 1]->getTags()->setString( "X-MIN", stringify(adfMinMax[0]).c_str() );
      C[band - 1]->getTags()->setString( "X-MAX", stringify(adfMinMax[1]).c_str() );
      for( size_t y = 0; y < nBlockYSize; y++ ) {
        if( poBand->RasterIO( GF_Read, 0, y, nBlockXSize, 1, pafScanline, 
		    nBlockXSize, 1, GDT_Float32, 0, 0) != CE_None ) {
          std::cerr << "input error" << std::endl;
          throw QuietException();
	}
	memcpy( C[band - 1]->getRawData() + y * nBlockXSize, pafScanline, nBlockXSize * sizeof(float) );
      }
    }
    CPLFree( pafScanline );
    GDALClose( poDataset );

    const char *fileNameTag = strcmp( "-", ff.fileName )==0 ? "stdin" : ff.fileName;
    frame->getTags()->setString( "FILE_NAME", fileNameTag );
        
    pfsio.writeFrame( frame, stdout );
    pfsio.freeFrame( frame );
  }
}

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

			QVariant xVar = QVariant::fromValue<int>(rasterX);
			QVariant yVar = QVariant::fromValue<int>(rasterY);
			pc->setMetaData("raster_width",xVar);
			pc->setMetaData("raster_height",yVar);
		}

		//fetch raster bands
		bool zRasterProcessed = false;
		unsigned zInvalid = 0;
		double zMinMax[2] = {0, 0};

		for (int i=1; i<=rasterCount; ++i)
		{
			ccLog::Print( "Reading band #%i", i);
			GDALRasterBand* poBand = poDataset->GetRasterBand(i);

			GDALColorInterp colorInterp = poBand->GetColorInterpretation();
			GDALDataType bandType = poBand->GetRasterDataType();

			int nBlockXSize, nBlockYSize;
			poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
			ccLog::Print( "Block=%dx%d Type=%s, ColorInterp=%s", nBlockXSize, nBlockYSize, GDALGetDataTypeName(poBand->GetRasterDataType()), GDALGetColorInterpretationName(colorInterp) );

			//fetching raster scan-line
			int nXSize = poBand->GetXSize();
			int nYSize = poBand->GetYSize();
			assert(nXSize == rasterX);
			assert(nYSize == rasterY);
			
			int bGotMin, bGotMax;
			double adfMinMax[2] = {0, 0};
			adfMinMax[0] = poBand->GetMinimum( &bGotMin );
			adfMinMax[1] = poBand->GetMaximum( &bGotMax );
			if (!bGotMin || !bGotMax )
				//DGM FIXME: if the file is corrupted (e.g. ASCII ArcGrid with missing rows) this method will enter in a infinite loop!
				GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
			ccLog::Print( "Min=%.3fd, Max=%.3f", adfMinMax[0], adfMinMax[1] );

			GDALColorTable* colTable = poBand->GetColorTable();
			if( colTable != NULL )
				printf( "Band has a color table with %d entries", colTable->GetColorEntryCount() );

			if( poBand->GetOverviewCount() > 0 )
				printf( "Band has %d overviews", poBand->GetOverviewCount() );

			if (colorInterp == GCI_Undefined && !zRasterProcessed/*&& !colTable*/) //probably heights?
			{
				zRasterProcessed = true;
				zMinMax[0] = adfMinMax[0];
				zMinMax[1] = adfMinMax[1];

				double* scanline = (double*) CPLMalloc(sizeof(double)*nXSize);
				//double* scanline = new double[nXSize];
				memset(scanline,0,sizeof(double)*nXSize);

				for (int j=0; j<nYSize; ++j)
				{
					if (poBand->RasterIO( GF_Read, /*xOffset=*/0, /*yOffset=*/j, /*xSize=*/nXSize, /*ySize=*/1, /*buffer=*/scanline, /*bufferSizeX=*/nXSize, /*bufferSizeY=*/1, /*bufferType=*/GDT_Float64, /*x_offset=*/0, /*y_offset=*/0 ) != CE_None)
					{
						delete pc;
						CPLFree(scanline);
						GDALClose(poDataset);
						return CC_FERR_READING;
					}

void generateTexture(string fname, GLuint& tex, int bandnum)
{
   if(bandnum <= 0 )
   {
     bandnum = 1;
   }
   GDALDataset *poDataset;
   GDALAllRegister();
   poDataset= (GDALDataset*) GDALOpen(fname.c_str(),GA_ReadOnly);
   if(poDataset == NULL)
   {
      cout << "OUCH!" << endl;
      //exit(0);
      return;
   }
   cout << "Data size: " << GDALGetRasterXSize(poDataset) << " " << GDALGetRasterYSize(poDataset) << endl;

   GDALRasterBand  *poBand;
   int             nBlockXSize, nBlockYSize;
   int             bGotMin, bGotMax;
   double          adfMinMax[2];
  int bands = poDataset->GetRasterCount();
   bandnum = bandnum % bands + 1;
   if(bandnum > bands)
   {
      bandnum = 1;
   }
   poBand = poDataset->GetRasterBand( bandnum );
   poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
   printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
            nBlockXSize, nBlockYSize,
            GDALGetDataTypeName(poBand->GetRasterDataType()),
            GDALGetColorInterpretationName(
            poBand->GetColorInterpretation()) );

   float max = adfMinMax[0] = poBand->GetMinimum( &bGotMin );
   float min = adfMinMax[1] = poBand->GetMaximum( &bGotMax );
   if( ! (bGotMin && bGotMax) )
      GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
   int width = poBand->GetXSize(); 
   int height = poBand->GetYSize();

   float *pafScanline;
   std::cout << "Before allocation" << adfMinMax[0] << " " << adfMinMax[1] << endl;
   min = adfMinMax[0];
   max = adfMinMax[1];
   int dsize = 256;
   pafScanline = (float *) CPLMalloc(sizeof(float)*512*512);
   vector<vector<float>> out = vector<vector<float>>(height,vector<float> (width,0));
   //vector<vector<unsigned char>> texs = vector<vector<unsigned char>>(height,vector<unsigned char> (width,0));
   unsigned char texs[512*512];
   poBand->RasterIO(GF_Read,0,0,width,height,pafScanline,512,512,GDT_Float32,0,0);
   float no = poBand->GetNoDataValue();
   cout << "After allocation" << endl;
   for(int i = 0; i < 512; i++)
   {
    for(int j = 0; j < 512; j++)
    {
      //cout << i << j << endl << pafS;
      if(pafScanline[i*width+j] != no)
      {
        // set tex val
        texs[i*512+j] = (unsigned char)(255*((pafScanline[i*512+j] - min)/(max-min)));
        //if((int)texs[i*width] < 0)
        //cout << (int)texs[i*512 +j] << " " << pafScanline[i*512+j] << " " << no << " " << fname << " " << min << " " << max << endl;
      }
      else
      {
        // Set zero val
        texs[i*512+j] = 0;
        //cout << (int)texs[i*512 +j] << fname << endl;
      }
      //texs[i*512+j] = 255;

      //ut[i][j] = pafScanline[i*width+j];
    }
   }
   CPLFree(pafScanline);
   //exit(0);
   // Create a texture
   glGenTextures(1,&tex);
   glBindTexture(GL_TEXTURE_2D,tex);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
   glTexImage2D(GL_TEXTURE_2D, 0, GL_RED, 512,512, 0, GL_RED, GL_UNSIGNED_BYTE,texs);
   GDALClose( (GDALDatasetH) poDataset);
   return;
} 

//.........这里部分代码省略.........
        const char *str = CSLFetchNameValue(papszOptions,"PAGESIZE");
        int newBlockX, newBlockY;
        sscanf(str, "%i,%i", &newBlockX, &newBlockY);
        if (newBlockX && newBlockY) {
            blockX = newBlockX;
            blockY = newBlockY;
        }
        else {
            CPLError(CE_Failure, CPLE_OpenFailed,
                     "FIT - Unable to parse option PAGESIZE values [%s]", str);
        }
    }

    // XXX - need to do lots of checking of block size
    // * provide ability to override block size with options
    // * handle non-square block size (like scanline)
    //   - probably default from non-tiled image - have default block size
    // * handle block size bigger than image size
    // * undesirable block size (non power of 2, others?)
    // * mismatched block sizes for different bands
    // * image that isn't even pages (i.e. partially empty pages at edge)
    CPLDebug("FIT write", "using block size %ix%i", blockX, blockY);

    head->xPageSize = blockX;
    gst_swapb(head->xPageSize);
    head->yPageSize = blockY;
    gst_swapb(head->yPageSize);
    head->zPageSize = 1;
    gst_swapb(head->zPageSize);
    head->cPageSize = nBands;
    gst_swapb(head->cPageSize);

    // XXX - need to check all bands
    head->minValue = firstBand->GetMinimum();
    gst_swapb(head->minValue);
    // XXX - need to check all bands
    head->maxValue = firstBand->GetMaximum();
    gst_swapb(head->maxValue);
    head->dataOffset = static_cast<unsigned int>(size);
    gst_swapb(head->dataOffset);

    CPL_IGNORE_RET_VAL(VSIFWriteL(head, size, 1, fpImage));

/* -------------------------------------------------------------------- */
/*      Loop over image, copying image data.                            */
/* -------------------------------------------------------------------- */
    unsigned long bytesPerPixel = nBands * nDTSize;

    size_t pageBytes = blockX * blockY * bytesPerPixel;
    char *output = (char *) malloc(pageBytes);
    if (! output)
    {
        CPLError(CE_Failure, CPLE_OutOfMemory,
                 "FITRasterBand couldn't allocate %lu bytes",
                 static_cast<unsigned long>(pageBytes));
        CPL_IGNORE_RET_VAL(VSIFCloseL(fpImage));
        return nullptr;
    }
    FreeGuard<char> guardOutput( output );

    long maxx = (long) ceil(poSrcDS->GetRasterXSize() / (double) blockX);
    long maxy = (long) ceil(poSrcDS->GetRasterYSize() / (double) blockY);
    long maxx_full = (long) floor(poSrcDS->GetRasterXSize() / (double) blockX);
    long maxy_full = (long) floor(poSrcDS->GetRasterYSize() / (double) blockY);

    CPLDebug("FIT", "about to write %ld x %ld blocks", maxx, maxy);

T* GISToFloatArray(char* fname, int interpWidth, int interpHeight)
{
  // Important note ------ Gdal considers images to be north up
  // the origin of datasets takes place in the upper-left or North-West corner.
  // Now to create a GDAL dataset
  // auto ds = ((GDALDataset*) GDALOpen(fname,GA_ReadOnly));
  GDALDataset* ds = ((GDALDataset*) GDALOpen(fname,GA_ReadOnly));
  if(ds == NULL)
  {
    return NULL;
  }
  
  // Creating a Raster band variable
  // A band represents one whole dataset within a dataset
  // in your case your files have one band.
  GDALRasterBand  *poBand;
  int             nBlockXSize, nBlockYSize;
  int             bGotMin, bGotMax;
  double          adfMinMax[2];
  
  // Assign the band      
  poBand = ds->GetRasterBand( 1 );
  poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );

  // find the min and max
  adfMinMax[0] = poBand->GetMinimum( &bGotMin );
  adfMinMax[1] = poBand->GetMaximum( &bGotMax );
  if( ! (bGotMin && bGotMax) )
    GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
  int min = adfMinMax[0];
  int max = adfMinMax[1];

  // get the width and height of the band or dataset
  int width = poBand->GetXSize();
  int height = poBand->GetYSize();

  // GDAL can handle files that have multiple datasets jammed witin it
  int bands = ds->GetRasterCount();

  // the float variable to hold the DEM!
  T *pafScanline;
  // std::cout << "Min: " << adfMinMax[0] << " Max: " << adfMinMax[1] << endl;
  int dsize = 256;
  // pafScanline = (T *) CPLMalloc(sizeof(T)*width*height);
  pafScanline = (T *) CPLMalloc(sizeof(T)*interpWidth*interpHeight);

  // Lets acquire the data.  ..... this funciton will interpolate for you
  // poBand->RasterIO(GF_Read,0,0,width,height,pafScanline,width,height,GDT_Float32,0,0);
  poBand->RasterIO(GF_Read,0,0,width,height,pafScanline,interpWidth,interpHeight,GDT_Float32,0,0);
  //        chage these two to interpolate automatically ^      ^

  // The Geotransform gives information on where a dataset is located in the world
  // and the resolution.
  // for more information look at http://www.gdal.org/gdal_datamodel.html
  double geot[6];
  ds->GetGeoTransform(geot);

  // Get the x resolution per pixel(south and west) and y resolution per pixel (north and south)
  // float xres = geot[1];
  // float yres = geot[5];
  // string proj;
  // proj = string(ds->GetProjectionRef());

  // You can get the projection string
  // The projection gives information about the coordinate system a dataset is in
  // This is important to allow other GIS softwares to place datasets into the same
  // coordinate space 
  // char* test = &proj[0];

  // The origin of the dataset 
  // float startx = geot[0]; // east - west coord.
  // float starty = geot[3]; // north - south coord.

  
  // here is some code that I used to push that 1D array into a 2D array
  // I believe this puts everything in the correct order....
  /*for(int i = 0; i < hsize; i++)
  {
    for(int j = 0; j < wsize; j++)
    {
      //cout << i << j << endl << pafS;
      vecs[i][j] = pafScanline[((int)i)*(int)wsize+((int)j)];
      if(vecs[i][j]>0 && vecs[i][j] > max)
      {
          max = vecs[i][j];
      }
      if(vecs[i][j]>0 && vecs[i][j] < min)
      {
          min = vecs[i][j];
      }
    }
   }*/
   //CPLFree(pafScanline);
   return pafScanline;

}

	bool DEM::load(const std::string& filename) {
		GDALAllRegister();

		GDALDataset* poDS;
		poDS = (GDALDataset*)GDALOpenEx(filename.c_str(), GDAL_OF_RASTER, NULL, NULL, NULL);
		if (poDS == NULL) return false;

		double adfGeoTransform[6];
		if (poDS->GetGeoTransform(adfGeoTransform) == CE_None) {
			origin.x = adfGeoTransform[0];
			origin.y = adfGeoTransform[3];

			pixelSize.x = adfGeoTransform[1];
			pixelSize.y = abs(adfGeoTransform[5]);
		}

		width = poDS->GetRasterXSize() * pixelSize.x;
		height = poDS->GetRasterYSize() * pixelSize.y;
		data.resize(width * height);
		min_val = std::numeric_limits<float>::max();
		max_val = -std::numeric_limits<float>::max();

		// bandが存在しない場合は、エラー
		if (poDS->GetRasterCount() == 0) return false;

		// 最初のbandのみを読み込む。複数bandは未対応
		GDALRasterBand* poBand = poDS->GetRasterBand(1);

		int nBlockXSize, nBlockYSize;
		poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
		//printf("Block=%dx%d Type=%s, ColorInterp=%s\n", nBlockXSize, nBlockYSize, GDALGetDataTypeName(poBand->GetRasterDataType()), GDALGetColorInterpretationName(poBand->GetColorInterpretation()));

		// 最低、最高の値を取得
		int bGotMin, bGotMax;
		double adfMinMax[2];
		adfMinMax[0] = poBand->GetMinimum(&bGotMin);
		adfMinMax[1] = poBand->GetMaximum(&bGotMax);
		if (!(bGotMin && bGotMax)) GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
		//printf("Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1]);
		min_val = adfMinMax[0];
		max_val = adfMinMax[1];
			
		//int nXSize = poBand->GetXSize();
		//int nYSize = poBand->GetYSize();

		int nXBlocks = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
		int nYBlocks = (poBand->GetYSize() + nBlockYSize - 1) / nBlockYSize;
		float *pData = (float *)CPLMalloc(sizeof(float*) * nBlockXSize * nBlockYSize);
		for (int iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
			for (int iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
				int nXValid, nYValid;

				poBand->ReadBlock(iXBlock, iYBlock, pData);

				// Compute the portion of the block that is valid
				// for partial edge blocks.
				if ((iXBlock + 1) * nBlockXSize > poBand->GetXSize())
					nXValid = poBand->GetXSize() - iXBlock * nBlockXSize;
				else
					nXValid = nBlockXSize;
				if ((iYBlock + 1) * nBlockYSize > poBand->GetYSize())
					nYValid = poBand->GetYSize() - iYBlock * nBlockYSize;
				else
					nYValid = nBlockYSize;

				for (int iY = 0; iY < nYValid; iY++) {
					for (int iX = 0; iX < nXValid; iX++) {
						float val;
						if (pData[iY * nBlockXSize + iX] > max_val) {
							val = max_val;
						}
						else if (pData[iY * nBlockXSize + iX] < min_val) {
							val = min_val;
						}
						else {
							val = pData[iY * nBlockXSize + iX];
						}

						for (int y = 0; y < pixelSize.y; ++y) {
							for (int x = 0; x < pixelSize.x; ++x) {
								data[((iYBlock * nBlockYSize + iY) * pixelSize.y + y) * width + (iXBlock * nBlockXSize + iX) * pixelSize.x + x] = val;
							}
						}
						
					}
				}
			}
		}

		GDALClose(poDS);

		return true;
	}

int readRaster(const char* pszFilename)
{
    // Set config options for GDAL (needs >= 2.0). Setting GTIFF_VIRTUAL_MEM_IO to "YES" can cause things bleed to
    // swap if enough RAM is not available. Use "IF_ENOUGH_RAM" for safer performance if unsure. NOTE that faster
    // mem io only works with *uncompressed* GeoTIFFs.
    // New in GDAL 2.0, from https://2015.foss4g-na.org/sites/default/files/slides/GDAL%202.0%20overview.pdf
    // GeoTIFF driver (with i7-4700 HQ (8 vCPUs)):with i7-4700 HQ (8 vCPUs)
    //   - Default: time ./testblockcache -ondisk: 7.5s
    //   - GTIFF_DIRECT_IO=YES: short circuit the block cache&libtiff for most RasterIO() operations (restricted to
    //     uncompressed stripped GeoTIFF). time ./testblockcache -ondisk: 2s
    //   - GTIFF_VIRTUAL_MEM_IO=YES: same as above, with tiled GeoTIFF as well. Uses memory-mapped file access. Linux
    //     only for now, 64bit recommended (could be extended to other platforms possible).
    //     time ./testblockcache -ondisk: 0.3s
    //
    CPLSetConfigOption("GTIFF_VIRTUAL_MEM_IO", "YES" );
    // Open the dataset
    GDALDataset *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );

    if( poDataset == NULL )
    {
        std::cout << "Dataset " << pszFilename << " could not be opened." << std::endl;
        return -1;
    }
    else
    {
        GDALRasterBand  *poBand;
        int             nBlockXSize, nBlockYSize;
        int             bGotMin, bGotMax;
        double          adfMinMax[2];

        //  Get raster band and its size
        poBand = poDataset->GetRasterBand(1);
        poBand->GetBlockSize( &nBlockXSize, &nBlockYSize);

        std::cout << "Dataset: " << pszFilename << std::endl;

        std::cout << "Block=" << nBlockXSize << "x" << nBlockYSize << " Type=" <<
                GDALGetDataTypeName(poBand->GetRasterDataType()) << " ColorInterp=" <<
                GDALGetColorInterpretationName(poBand->GetColorInterpretation()) << std::endl;

        // Calculate some stats
        adfMinMax[0] = poBand->GetMinimum(&bGotMin);
        adfMinMax[1] = poBand->GetMaximum(&bGotMax);

        if(!(bGotMin && bGotMax)) {
            GDALComputeRasterMinMax((GDALRasterBandH) poBand, TRUE, adfMinMax);
        }
        std::cout << "Min=" << adfMinMax[0] << " Max=" << adfMinMax[1] << std::endl;

        if(poBand->GetOverviewCount() > 0) {
            std::cout << "Band has " << poBand->GetOverviewCount() << " overviews." << std::endl;
        }
        if( poBand->GetColorTable() != NULL ) {
            std::cout << "Band has a color table with " << poBand->GetColorTable()->GetColorEntryCount() <<
                    " entries." << std::endl;
        }

        // Get the actual data
        float *pafScanline;
        int   nXSize = poBand->GetXSize();
        pafScanline = (float *) CPLMalloc(sizeof(float) * nXSize);

        // RasterIO has a new argument psExtraArg in GDAL > 2.0. NOTE: that  GDALRasterBand::ReadBlock() probably has
        // better performance for reading the whole data at one go.
        #ifdef USE_GDAL_2
            GDALRasterIOExtraArg* arg = NULL;
            poBand->RasterIO(GF_Read, 0, 0, nXSize, 1, pafScanline, nXSize, 1, GDT_Float3GDALRasterBand::ReadBlock 2,
               0, 0, arg);
        #else
            poBand->RasterIO(GF_Read, 0, 0, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
        #endif

        // ... do something with the data ...

        // Free resources
        CPLFree(pafScanline);
        GDALClose((GDALDatasetH) poDataset);
        std::cout << std::endl;
    }
    return 0;
}