C++ GDALColorTable类代码示例

void ReadLUT(GDALRasterBand* rasterBand, TeRasterParams& params)
{
	GDALColorTable* colorTable = rasterBand->GetColorTable();
	if(colorTable == 0)
		return;

	int ncolors = colorTable->GetColorEntryCount();
	for(int i = 0; i < ncolors; ++i)
	{
		const GDALColorEntry* ce = colorTable->GetColorEntry(i);
		if(ce == 0)
			continue;

		params.lutr_.push_back(ce->c1); 
		params.lutg_.push_back(ce->c2);
		params.lutb_.push_back(ce->c3);
	}
}

GDALColorTable *GDALProxyPoolRasterBand::GetColorTable()
{
    GDALRasterBand* poUnderlyingRasterBand = RefUnderlyingRasterBand();
    if (poUnderlyingRasterBand == NULL)
        return NULL;

    if (poColorTable)
        delete poColorTable;
    poColorTable = NULL;

    GDALColorTable* poUnderlyingColorTable = poUnderlyingRasterBand->GetColorTable();
    if (poUnderlyingColorTable)
        poColorTable = poUnderlyingColorTable->Clone();

    UnrefUnderlyingRasterBand(poUnderlyingRasterBand);

    return poColorTable;
}

 CALSWrapperSrcBand(GDALDataset* poSrcDSIn)
 {
     this->poSrcDS = poSrcDSIn;
     SetMetadataItem("NBITS", "1", "IMAGE_STRUCTURE");
     poSrcDS->GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
     eDataType = GDT_Byte;
     bInvertValues = TRUE;
     GDALColorTable* poCT = poSrcDS->GetRasterBand(1)->GetColorTable();
     if( poCT != NULL && poCT->GetColorEntryCount() >= 2 )
     {
         const GDALColorEntry* psEntry1 = poCT->GetColorEntry(0);
         const GDALColorEntry* psEntry2 = poCT->GetColorEntry(1);
         if( psEntry1->c1 == 255 && psEntry1->c2 == 255 && psEntry1->c3 == 255 &&
             psEntry2->c1 == 0 && psEntry2->c2 == 0 && psEntry2->c3 == 0 )
         {
             bInvertValues = FALSE;
         }
     }
 }

GDALDataset *
GIFDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
               int bStrict, char ** papszOptions, 
               GDALProgressFunc pfnProgress, void * pProgressData )

{
    int  nBands = poSrcDS->GetRasterCount();
    int  nXSize = poSrcDS->GetRasterXSize();
    int  nYSize = poSrcDS->GetRasterYSize();
    int	 bInterlace = FALSE;

/* -------------------------------------------------------------------- */
/*      Check for interlaced option.                                    */
/* -------------------------------------------------------------------- */
    bInterlace = CSLFetchBoolean(papszOptions, "INTERLACING", FALSE);

/* -------------------------------------------------------------------- */
/*      Some some rudimentary checks                                    */
/* -------------------------------------------------------------------- */
    if( nBands != 1 )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "GIF driver only supports one band images.\n" );

        return NULL;
    }

    if (nXSize > 65535 || nYSize > 65535)
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "GIF driver only supports datasets up to 65535x65535 size.\n" );

        return NULL;
    }

    if( poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte 
        && bStrict )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "GIF driver doesn't support data type %s. "
                  "Only eight bit bands supported.\n", 
                  GDALGetDataTypeName( 
                      poSrcDS->GetRasterBand(1)->GetRasterDataType()) );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Open the output file.                                           */
/* -------------------------------------------------------------------- */
    GifFileType *hGifFile;
    VSILFILE *fp;

    fp = VSIFOpenL( pszFilename, "wb" );
    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "Failed to create %s:\n%s", 
                  pszFilename, VSIStrerror( errno ) );
        return NULL;
    }

#if defined(GIFLIB_MAJOR) && GIFLIB_MAJOR >= 5
    int nError;
    hGifFile = EGifOpen( fp, VSIGIFWriteFunc, &nError );
#else
    hGifFile = EGifOpen( fp, VSIGIFWriteFunc );
#endif
    if( hGifFile == NULL )
    {
        VSIFCloseL( fp );
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "EGifOpenFilename(%s) failed.  Does file already exist?",
                  pszFilename );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Prepare colortable.                                             */
/* -------------------------------------------------------------------- */
    GDALRasterBand	*poBand = poSrcDS->GetRasterBand(1);
    ColorMapObject	*psGifCT;
    int			iColor;

    if( poBand->GetColorTable() == NULL )
    {
        psGifCT = GifMakeMapObject( 256, NULL );
        for( iColor = 0; iColor < 256; iColor++ )
        {
            psGifCT->Colors[iColor].Red = (GifByteType) iColor;
            psGifCT->Colors[iColor].Green = (GifByteType) iColor;
            psGifCT->Colors[iColor].Blue = (GifByteType) iColor;
        }
    }
    else
    {
        GDALColorTable	*poCT = poBand->GetColorTable();
        int nFullCount = 1;

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

//.........这里部分代码省略.........
    poDS->adfGeoTransform[0] = fTmp;
    memcpy(&fTmp, poDS->pachHeader + 120, 4);
    CPL_LSBPTR32(&fTmp);
    poDS->adfGeoTransform[1] = fTmp;
    poDS->adfGeoTransform[2] = 0.0;
    memcpy(&fTmp, poDS->pachHeader + 116, 4);
    CPL_LSBPTR32(&fTmp);
    poDS->adfGeoTransform[3] = fTmp;
    poDS->adfGeoTransform[4] = 0.0;
    memcpy(&fTmp, poDS->pachHeader + 124, 4);
    CPL_LSBPTR32(&fTmp);
    poDS->adfGeoTransform[5] = - fTmp;

    // adjust for center of pixel vs. top left corner of pixel.
    poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[1] * 0.5;
    poDS->adfGeoTransform[3] -= poDS->adfGeoTransform[5] * 0.5;

/* -------------------------------------------------------------------- */
/*      If we didn't get any georeferencing, try for a worldfile.       */
/* -------------------------------------------------------------------- */
    if( poDS->adfGeoTransform[1] == 0.0
        || poDS->adfGeoTransform[5] == 0.0 )
    {
        if( !GDALReadWorldFile( poOpenInfo->pszFilename, NULL,
                                poDS->adfGeoTransform ) )
            GDALReadWorldFile( poOpenInfo->pszFilename, ".wld",
                               poDS->adfGeoTransform );
    }

/* -------------------------------------------------------------------- */
/*      Try to come up with something for the coordinate system.        */
/* -------------------------------------------------------------------- */
    memcpy(&nTmp16, poDS->pachHeader + 88, 2);
    CPL_LSBPTR16(&nTmp16);
    int nCoordSys = nTmp16;

    if( nCoordSys == 0 )
    {
        poDS->pszProjection = CPLStrdup(SRS_WKT_WGS84);
    }
    else if( nCoordSys == 1 )
    {
        poDS->pszProjection =
            CPLStrdup( "LOCAL_CS[\"UTM - Zone Unknown\",UNIT[\"Meter\",1]]" );
    }
    else if( nCoordSys == 2 )
    {
        poDS->pszProjection =
            CPLStrdup(
                "LOCAL_CS[\"State Plane - Zone Unknown\","
                "UNIT[\"US survey foot\",0.3048006096012192]]" );
    }
    else
    {
        poDS->pszProjection =
            CPLStrdup( "LOCAL_CS[\"Unknown\",UNIT[\"Meter\",1]]" );
    }

/* -------------------------------------------------------------------- */
/*      Check for a trailer file with a colormap in it.                 */
/* -------------------------------------------------------------------- */
    char *pszPath = CPLStrdup( CPLGetPath(poOpenInfo->pszFilename) );
    char *pszBasename = CPLStrdup( CPLGetBasename(poOpenInfo->pszFilename) );
    const char *pszTRLFilename =
        CPLFormCIFilename( pszPath, pszBasename, "trl" );
    VSILFILE *fpTRL = VSIFOpenL( pszTRLFilename, "rb" );
    if( fpTRL != NULL )
    {
        char szTRLData[896] = { '\0' };

        CPL_IGNORE_RET_VAL(VSIFReadL( szTRLData, 1, 896, fpTRL ));
        CPL_IGNORE_RET_VAL(VSIFCloseL( fpTRL ));

        GDALColorTable *poCT = new GDALColorTable();
        for( int iColor = 0; iColor < 256; iColor++ )
        {
            GDALColorEntry sEntry = { 0, 0, 0, 0};

            sEntry.c2 = reinterpret_cast<GByte *>(szTRLData)[iColor+128];
            sEntry.c1 = reinterpret_cast<GByte *>(szTRLData)[iColor+128+256];
            sEntry.c3 = reinterpret_cast<GByte *>(szTRLData)[iColor+128+512];
            sEntry.c4 = 255;
            poCT->SetColorEntry( iColor, &sEntry );

            // Only 16 colors in 4bit files.
            if( nPixelOffset == -1 && iColor == 15 )
                break;
        }

        poDS->GetRasterBand(1)->SetColorTable( poCT );
        poDS->GetRasterBand(1)->SetColorInterpretation( GCI_PaletteIndex );

        delete poCT;
    }

    CPLFree( pszPath );
    CPLFree( pszBasename );

    return poDS;
}

void output_pcolor_geotiff ( image<short> & img, image<float> & dist,
                             image<unsigned char> & v )
{
    int r, c;
    int val;
    OGRSpatialReference ref;
    GDALDataset  *df;
    char *wkt = NULL;
    GDALRasterBand *bd;
    double        trans[6];
    GDALDriver   *gt;
    char         **options = NULL;
    int          ov[] = { 2, 4, 8, 16, 32 };
    int          nov;
    int          i, n;
    char         file[1000];
    GDALColorTable ct;
    GDALColorEntry color;

    color.c1 = 0;
    color.c2 = 0;
    color.c3 = 0;
    color.c4 = 255;
    ct.SetColorEntry ( 0, &color );
    for ( i = 1; i < 256; i++ ) {
        color.c1 = 30 + rand() % 240 - rand()%30;
        color.c2 = 20 + rand() % 200 - rand()%20;
        color.c3 = 0 + rand() % 160;
        ct.SetColorEntry ( i, &color );
    }

    options = CSLSetNameValue ( options, "TILED", "NO" );
    options = CSLSetNameValue ( options, "COMPRESS", "LZW" );

    gt = GetGDALDriverManager()->GetDriverByName("GTiff");
    if ( !gt ) {
        fprintf(stderr,"Could not get GTiff driver\n");
        exit(1);
    }

    strcpy ( file, p.value("blob_file").c_str() );
    df = gt->Create( file, dist.cols, dist.rows, 1, GDT_Byte, options );
    if( df == NULL ) {
        fprintf(stderr,"Could not create %s\n", file );
        exit(1);
    }

    trans[0] = p.dvalue("easting_left");
    trans[1] = p.dvalue("output_cell_size");
    trans[2] = 0.0;
    trans[3] = p.dvalue("northing_top");
    trans[4] = 0.0;
    trans[5] = -p.dvalue("output_cell_size");
    df->SetGeoTransform ( trans );
    ref.SetUTM ( p.ivalue("utm_zone") );
    ref.SetWellKnownGeogCS ( "NAD27" );
    ref.exportToWkt ( &wkt );
    df->SetProjection(wkt);
    CPLFree ( wkt );

    for ( r=0; r < dist.rows; r++ ) {
        for ( c=0; c < dist.cols; c++ ) {
            if ( dist[r][c] == 0 ) v[r][c] = 0;
            else v[r][c] = img[r][c] % 256;
        }
    }

    bd = df->GetRasterBand(1);
    bd->SetColorTable ( &ct );
    bd->RasterIO( GF_Write, 0, 0, v.cols, v.rows, v.data,
                  v.cols, v.rows, GDT_Byte, 0, 0 );

    delete df;

    df = (GDALDataset *)GDALOpen ( file, GA_Update );
    if( df == NULL ) {
        fprintf(stderr,"Could not open for update %s\n", file );
        exit(1);
    }
    nov = p.ivalue("overviews");
    if ( nov > 5 ) nov = 5;
    if ( nov > 0 ) {
        n = 1;
        df->BuildOverviews("NEAREST", nov, ov, 1, &n, NULL, NULL );
    }

}

static GDALDataset *
BSBCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
               int bStrict, char ** papszOptions, 
               GDALProgressFunc pfnProgress, void * pProgressData )

{
    int  nBands = poSrcDS->GetRasterCount();
    int  nXSize = poSrcDS->GetRasterXSize();
    int  nYSize = poSrcDS->GetRasterYSize();

/* -------------------------------------------------------------------- */
/*      Some some rudimentary checks                                    */
/* -------------------------------------------------------------------- */
    if( nBands != 1 )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "BSB driver only supports one band images.\n" );

        return NULL;
    }

    if( poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte 
        && bStrict )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "BSB driver doesn't support data type %s. "
                  "Only eight bit bands supported.\n", 
                  GDALGetDataTypeName( 
                      poSrcDS->GetRasterBand(1)->GetRasterDataType()) );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Open the output file.                                           */
/* -------------------------------------------------------------------- */
    BSBInfo *psBSB;

    psBSB = BSBCreate( pszFilename, 0, 200, nXSize, nYSize );
    if( psBSB == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Prepare initial color table.colortable.                         */
/* -------------------------------------------------------------------- */
    GDALRasterBand	*poBand = poSrcDS->GetRasterBand(1);
    int			iColor;
    unsigned char       abyPCT[771];
    int                 nPCTSize;
    int                 anRemap[256];

    abyPCT[0] = 0;
    abyPCT[1] = 0;
    abyPCT[2] = 0;

    if( poBand->GetColorTable() == NULL )
    {
        /* map greyscale down to 63 grey levels. */
        for( iColor = 0; iColor < 256; iColor++ )
        {
            int nOutValue = (int) (iColor / 4.1) + 1;

            anRemap[iColor] = nOutValue;
            abyPCT[nOutValue*3 + 0] = (unsigned char) iColor;
            abyPCT[nOutValue*3 + 1] = (unsigned char) iColor;
            abyPCT[nOutValue*3 + 2] = (unsigned char) iColor;
        }
        nPCTSize = 64;
    }
    else
    {
        GDALColorTable	*poCT = poBand->GetColorTable();
        int nColorTableSize = poCT->GetColorEntryCount();
        if (nColorTableSize > 255)
            nColorTableSize = 255;

        for( iColor = 0; iColor < nColorTableSize; iColor++ )
        {
            GDALColorEntry	sEntry;

            poCT->GetColorEntryAsRGB( iColor, &sEntry );

            anRemap[iColor] = iColor + 1;
            abyPCT[(iColor+1)*3 + 0] = (unsigned char) sEntry.c1;
            abyPCT[(iColor+1)*3 + 1] = (unsigned char) sEntry.c2;
            abyPCT[(iColor+1)*3 + 2] = (unsigned char) sEntry.c3;
        }

        nPCTSize = nColorTableSize + 1;

        // Add entries for pixel values which apparently will not occur.
        for( iColor = nPCTSize; iColor < 256; iColor++ )
            anRemap[iColor] = 1;
    }

/* -------------------------------------------------------------------- */
/*      Boil out all duplicate entries.                                 */
/* -------------------------------------------------------------------- */
    int  i;

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

//.........这里部分代码省略.........
            QString srsProj = QString(theProj4);
            if (!srsProj.isEmpty() && theProjection != srsProj) {
                cleanup();
                theProjection = srsProj;
            }
            isLatLon = (theSrs->IsGeographic() == TRUE);
        }
    }
    if (theProjection.isEmpty()) {
        theProjection = ProjectionChooser::getProjection(QCoreApplication::translate("ImportExportGdal", "Unable to set projection; please specify one"));
        if (theProjection.isEmpty()) {
            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:

//.........这里部分代码省略.........
    AttachMetadata((GDALDatasetH) poVDS, papszMetadataOptions);

    const char *pszInterleave = GDALGetMetadataItem(hDataset, "INTERLEAVE", "IMAGE_STRUCTURE");
    if (pszInterleave)
        poVDS->SetMetadataItem("INTERLEAVE", pszInterleave, "IMAGE_STRUCTURE");

/* -------------------------------------------------------------------- */
/*      Transfer metadata that remains valid if the spatial             */
/*      arrangement of the data is unaltered.                           */
/* -------------------------------------------------------------------- */
    if (bSpatialArrangementPreserved)
    {
        char **papszMD;

        papszMD = ((GDALDataset*)hDataset)->GetMetadata("RPC");
        if (papszMD != NULL)
            poVDS->SetMetadata(papszMD, "RPC");

        papszMD = ((GDALDataset*)hDataset)->GetMetadata("GEOLOCATION");
        if (papszMD != NULL)
            poVDS->SetMetadata(papszMD, "GEOLOCATION");
    }

    int nSrcBandCount = nBandCount;

    if (nRGBExpand != 0)
    {
        GDALRasterBand *poSrcBand;
        poSrcBand = ((GDALDataset*)
                     hDataset)->GetRasterBand(ABS(panBandList[0]));
        if (panBandList[0] < 0)
            poSrcBand = poSrcBand->GetMaskBand();

        GDALColorTable *poColorTable = poSrcBand->GetColorTable();
        if (poColorTable == NULL)
        {
            fprintf(stderr, "Error : band %d has no color table\n", ABS(panBandList[0]));
            GDALClose(hDataset);
            CPLFree(panBandList);
            GDALDestroyDriverManager();
            CSLDestroy(argv);
            CSLDestroy(papszCreateOptions);
            exit(1);
        }

        /* Check that the color table only contains gray levels */
        /* when using -expand gray */
        if (nRGBExpand == 1)
        {
            int nColorCount = poColorTable->GetColorEntryCount();
            int nColor;

            for (nColor = 0; nColor < nColorCount; nColor++)
            {
                const GDALColorEntry *poEntry = poColorTable->GetColorEntry(nColor);
                if (poEntry->c1 != poEntry->c2 || poEntry->c1 != poEntry->c2)
                {
                    fprintf(stderr, "Warning : color table contains non gray levels colors\n");
                    break;
                }
            }
        }

        if (nBandCount == 1)
            nBandCount = nRGBExpand;
        else if (nBandCount == 2 && (nRGBExpand == 3 || nRGBExpand == 4))

/**
* read data
*/
bool KGDAL2CV::readData(cv::Mat img){
	// make sure the image is the proper size
	if (img.size().height != m_height){
		return false;
	}
	if (img.size().width != m_width){
		return false;
	}

	// make sure the raster is alive
	if (m_dataset == NULL || m_driver == NULL){
		return false;
	}

	// set the image to zero
	img = 0;

	// iterate over each raster band
	// note that OpenCV does bgr rather than rgb
	int nChannels = m_dataset->GetRasterCount();

	GDALColorTable* gdalColorTable = NULL;
	if (m_dataset->GetRasterBand(1)->GetColorTable() != NULL){
		gdalColorTable = m_dataset->GetRasterBand(1)->GetColorTable();
	}

	const GDALDataType gdalType = m_dataset->GetRasterBand(1)->GetRasterDataType();
	int nRows, nCols;

	//if (nChannels > img.channels()){
	//	nChannels = img.channels();
	//}

	for (int c = 0; c < img.channels(); c++){

		int realBandIndex = c;
		
		// get the GDAL Band
		GDALRasterBand* band = m_dataset->GetRasterBand(c + 1);

		//if (hasColorTable == false){
		if (GCI_RedBand == band->GetColorInterpretation()) realBandIndex = 2;
		if (GCI_GreenBand == band->GetColorInterpretation()) realBandIndex = 1;
		if (GCI_BlueBand == band->GetColorInterpretation()) realBandIndex = 0;
		//}
		if (hasColorTable && gdalColorTable->GetPaletteInterpretation() == GPI_RGB) c = img.channels() - 1;
		// make sure the image band has the same dimensions as the image
		if (band->GetXSize() != m_width || band->GetYSize() != m_height){ return false; }

		// grab the raster size
		nRows = band->GetYSize();
		nCols = band->GetXSize();

		// create a temporary scanline pointer to store data
		double* scanline = new double[nCols];

		// iterate over each row and column
		for (int y = 0; y<nRows; y++){

			// get the entire row
			band->RasterIO(GF_Read, 0, y, nCols, 1, scanline, nCols, 1, GDT_Float64, 0, 0);

			// set inside the image
			for (int x = 0; x<nCols; x++){

				// set depending on image types
				// given boost, I would use enable_if to speed up.  Avoid for now.
				if (hasColorTable == false){
					write_pixel(scanline[x], gdalType, nChannels, img, y, x, realBandIndex);
				}
				else{
					write_ctable_pixel(scanline[x], gdalType, gdalColorTable, img, y, x, c);
				}
			}
		}
		// delete our temp pointer
		delete[] scanline;
	}

	return true;
}

cv::Mat KGDAL2CV::ImgReadByGDAL(cv::String filename, int xStart, int yStart, int xWidth, int yWidth, bool beReadFourth)
{
	m_filename = filename;
	if (!readHeader()) return cv::Mat();
	
	int tempType = m_type;

	if (xStart < 0 || yStart < 0 || xWidth < 1 || yWidth < 1 || xStart > m_width - 1 || yStart > m_height - 1) return cv::Mat();

	if (xStart + xWidth > m_width)
	{
		std::cout << "The specified width is invalid, Automatic optimization is executed!" << std::endl;
		xWidth = m_width - xStart;
	}

	if (yStart + yWidth > m_height)
	{
		std::cout << "The specified height is invalid, Automatic optimization is executed!" << std::endl;
		yWidth = m_height - yStart;
	}

	if (!beReadFourth && 4 == m_nBand)
	{
		for (unsigned int index = CV_8S; index < CV_USRTYPE1; ++index)
		{
			if (CV_MAKETYPE(index, m_nBand) == m_type)
			{
				std::cout << "We won't read the fourth band unless it's datatype is GDT_Byte!" << std::endl;
				//tempType = -1;
				tempType = tempType - ((3 << CV_CN_SHIFT) - (2 << CV_CN_SHIFT));
				break;
				//tempType = CV_MAKETYPE(index, m_nBand);
			}
		}
	}

	cv::Mat img(yWidth, xWidth, tempType, cv::Scalar::all(0.f));
	// iterate over each raster band
	// note that OpenCV does bgr rather than rgb
	int nChannels = m_dataset->GetRasterCount();

	GDALColorTable* gdalColorTable = NULL;
	if (m_dataset->GetRasterBand(1)->GetColorTable() != NULL){
		gdalColorTable = m_dataset->GetRasterBand(1)->GetColorTable();
	}

	const GDALDataType gdalType = m_dataset->GetRasterBand(1)->GetRasterDataType();

	//if (nChannels > img.channels()){
	//	nChannels = img.channels();
	//}
	for (int c = 0; c < img.channels(); c++){

		int realBandIndex = c;

		// get the GDAL Band
		GDALRasterBand* band = m_dataset->GetRasterBand(c + 1);

		//if (hasColorTable == false){
		if (GCI_RedBand == band->GetColorInterpretation()) realBandIndex = 2;
		if (GCI_GreenBand == band->GetColorInterpretation()) realBandIndex = 1;
		if (GCI_BlueBand == band->GetColorInterpretation()) realBandIndex = 0;
		//}
		if (hasColorTable && gdalColorTable->GetPaletteInterpretation() == GPI_RGB) c = img.channels() - 1;
		// make sure the image band has the same dimensions as the image
		if (band->GetXSize() != m_width || band->GetYSize() != m_height){ return cv::Mat(); }

		// create a temporary scanline pointer to store data
		double* scanline = new double[xWidth];

		// iterate over each row and column
		for (int y = 0; y<yWidth; y++){

			// get the entire row
			band->RasterIO(GF_Read, xStart, y + yStart, xWidth, 1, scanline, xWidth, 1, GDT_Float64, 0, 0);
			// set inside the image
			for (int x = 0; x<xWidth; x++){
				// set depending on image types
				// given boost, I would use enable_if to speed up.  Avoid for now.
				if (hasColorTable == false){
					write_pixel(scanline[x], gdalType, nChannels, img, y, x, realBandIndex);
				}
				else{
					write_ctable_pixel(scanline[x], gdalType, gdalColorTable, img, y, x, c);
				}
			}
		}
		// delete our temp pointer
		delete[] scanline;
	}
	return img;
}

static unsigned char *
ParseXPM( const char *pszInput, int *pnXSize, int *pnYSize, 
          GDALColorTable **ppoRetTable )

{
/* ==================================================================== */
/*      Parse input into an array of strings from within the first C    */
/*      initializer (list os comma separated strings in braces).        */
/* ==================================================================== */
    char **papszXPMList = NULL;
    const char *pszNext = pszInput;
    int  i;

    // Skip till after open brace.
    while( *pszNext != '\0' && *pszNext != '{' ) 
        pszNext++;

    if( *pszNext == '\0' )
        return NULL;

    pszNext++;

    // Read lines till close brace.
    
    while( *pszNext != '\0' && *pszNext != '}' )
    {
        // skip whole comment. 
        if( EQUALN(pszNext,"/*",2) )
        {
            pszNext += 2;
            while( *pszNext != '\0' && !EQUALN(pszNext,"*/",2) )
                pszNext++;
        }

        // reading string constants
        else if( *pszNext == '"' )
        {
            char   *pszLine;

            pszNext++;
            i = 0;

            while( pszNext[i] != '\0' && pszNext[i] != '"' )
                i++;

            if( pszNext[i] == '\0' )
            {
                CSLDestroy( papszXPMList );
                return NULL;
            }

            pszLine = (char *) CPLMalloc(i+1);
            strncpy( pszLine, pszNext, i );
            pszLine[i] = '\0';

            papszXPMList = CSLAddString( papszXPMList, pszLine );
            CPLFree( pszLine );
            pszNext = pszNext + i + 1;
        }

        // just ignore everything else (whitespace, commas, newlines, etc).
        else
            pszNext++;
    }

    if( CSLCount(papszXPMList) < 3 || *pszNext != '}' )
    {
        CSLDestroy( papszXPMList );
        return NULL;
    }
    
/* -------------------------------------------------------------------- */
/*      Get the image information.                                      */
/* -------------------------------------------------------------------- */
    int nColorCount, nCharsPerPixel;

    if( sscanf( papszXPMList[0], "%d %d %d %d", 
                pnXSize, pnYSize, &nColorCount, &nCharsPerPixel ) != 4 )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Image definition (%s) not well formed.",
                  papszXPMList[0] );
        CSLDestroy( papszXPMList );
        return NULL;
    }

    if( nCharsPerPixel != 1 )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Only one character per pixel XPM images supported by GDAL at this time." );
        CSLDestroy( papszXPMList );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Parse out colors.                                               */
/* -------------------------------------------------------------------- */
    int iColor;
    int anCharLookup[256];
    GDALColorTable oCTable;
//.........这里部分代码省略.........

static GDALDataset *
XPMCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
               int bStrict, char ** papszOptions, 
               GDALProgressFunc pfnProgress, void * pProgressData )

{
    int  nBands = poSrcDS->GetRasterCount();
    int  nXSize = poSrcDS->GetRasterXSize();
    int  nYSize = poSrcDS->GetRasterYSize();
    GDALColorTable *poCT;

/* -------------------------------------------------------------------- */
/*      Some some rudimentary checks                                    */
/* -------------------------------------------------------------------- */
    if( nBands != 1 )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "XPM driver only supports one band images.\n" );

        return NULL;
    }

    if( poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte 
        && bStrict )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "XPM driver doesn't support data type %s. "
                  "Only eight bit bands supported.\n", 
                  GDALGetDataTypeName( 
                      poSrcDS->GetRasterBand(1)->GetRasterDataType()) );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      If there is no colortable on the source image, create a         */
/*      greyscale one with 64 levels of grey.                           */
/* -------------------------------------------------------------------- */
    GDALRasterBand	*poBand = poSrcDS->GetRasterBand(1);
    int                 i;
    GDALColorTable      oGreyTable;

    poCT = poBand->GetColorTable();
    if( poCT == NULL )
    {
        poCT = &oGreyTable;

        for( i = 0; i < 256; i++ )
        {
            GDALColorEntry sColor;

            sColor.c1 = (short) i;
            sColor.c2 = (short) i;
            sColor.c3 = (short) i;
            sColor.c4 = 255;

            poCT->SetColorEntry( i, &sColor );
        }
    }

/* -------------------------------------------------------------------- */
/*      Build list of active colors, and the mapping from pixels to     */
/*      our active colormap.                                            */
/* -------------------------------------------------------------------- */
    const char *pszColorCodes = " [email protected]#$%^&*()-+=[]|:;,.<>?/";

    int  anPixelMapping[256];
    GDALColorEntry asPixelColor[256];
    int  nActiveColors = MIN(poCT->GetColorEntryCount(),256);

    // Setup initial colortable and pixel value mapping. 
    memset( anPixelMapping+0, 0, sizeof(int) * 256 );
    for( i = 0; i < nActiveColors; i++ )
    {
        poCT->GetColorEntryAsRGB( i, asPixelColor + i );
        anPixelMapping[i] = i;
    }

/* ==================================================================== */
/*      Iterate merging colors until we are under our limit (about 85). */
/* ==================================================================== */
    while( nActiveColors > (int) strlen(pszColorCodes) )
    {
        int nClosestDistance = 768;
        int iClose1 = -1, iClose2 = -1;
        int iColor1, iColor2;

        // Find the closest pair of colors. 
        for( iColor1 = 0; iColor1 < nActiveColors; iColor1++ )
        {
            for( iColor2 = iColor1+1; iColor2 < nActiveColors; iColor2++ )
            {
                int	nDistance;

                if( asPixelColor[iColor1].c4 < 128 
                    && asPixelColor[iColor2].c4 < 128 )
                    nDistance = 0;
                else
                    nDistance = 
                        ABS(asPixelColor[iColor1].c1-asPixelColor[iColor2].c1)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    }
/* -------------------------------------------------------------------- */
/*      Otherwise just open it for update access.                       */
/* -------------------------------------------------------------------- */
    else
    {
        fpL = VSIFOpenL( pszFilename, "r+" );
        if( fpL == nullptr )
            hOTIFF = nullptr;
        else
            hOTIFF = VSI_TIFFOpen( pszFilename, "r+", fpL );
        if( hOTIFF == nullptr )
        {
            if( CPLGetLastErrorNo() == 0 )
                CPLError( CE_Failure, CPLE_OpenFailed,
                          "Attempt to create new tiff file `%s' "
                          "failed in VSI_TIFFOpen().",
                          pszFilename );
            if( fpL != nullptr )
                CPL_IGNORE_RET_VAL(VSIFCloseL(fpL));
            return CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      Do we have a palette?  If so, create a TIFF compatible version. */
/* -------------------------------------------------------------------- */
    unsigned short *panRed = nullptr;
    unsigned short *panGreen = nullptr;
    unsigned short *panBlue = nullptr;

    if( nPhotometric == PHOTOMETRIC_PALETTE )
    {
        GDALColorTable *poCT = papoBandList[0]->GetColorTable();
        int nColorCount = 65536;

        if( nBitsPerPixel <= 8 )
            nColorCount = 256;

        panRed = static_cast<unsigned short *>(
            CPLCalloc(nColorCount, sizeof(unsigned short)) );
        panGreen = static_cast<unsigned short *>(
            CPLCalloc(nColorCount, sizeof(unsigned short)) );
        panBlue = static_cast<unsigned short *>(
            CPLCalloc(nColorCount, sizeof(unsigned short)) );

        for( int iColor = 0; iColor < nColorCount; iColor++ )
        {
          GDALColorEntry sRGB = { 0, 0, 0, 0 };

            if( poCT->GetColorEntryAsRGB( iColor, &sRGB ) )
            {
                // TODO(schwehr): Check for underflow.
                // Going from signed short to unsigned short.
                panRed[iColor] = static_cast<unsigned short>(257 * sRGB.c1);
                panGreen[iColor] = static_cast<unsigned short>(257 * sRGB.c2);
                panBlue[iColor] = static_cast<unsigned short>(257 * sRGB.c3);
            }
        }
    }

/* -------------------------------------------------------------------- */
/*      Do we need some metadata for the overviews?                     */
/* -------------------------------------------------------------------- */
    CPLString osMetadata;
    GDALDataset *poBaseDS = papoBandList[0]->GetDataset();

//.........这里部分代码省略.........
    // RGB or RGBA are same
    // 

    if ((image.pagesize.c==1)&&(node=CPLGetXMLNode(defimage,"Palette"))) {
         int entries=static_cast<int>(getXMLNum(node,"Size",255));
         GDALPaletteInterp eInterp=GPI_RGB;
         // A flag to convert from HLS to HSV
         bool is_hsv=false;
         CPLString pModel=CPLGetXMLValue(node,"Model","RGB");
         if (!pModel.empty()) {
             if (pModel.find("HSV")!=string::npos) {
                 eInterp=GPI_HLS;
                 is_hsv=true;
             } else if (pModel.find("HLS")!=string::npos)
                 eInterp=GPI_HLS;
             else if (pModel.find("CMYK")!=string::npos) eInterp=GPI_CMYK;
             // Can it do LuminanceAlpha?
             else if (pModel.find("L")!=string::npos) eInterp=GPI_Gray;
             // RGBA and RGB are the same
             else if (pModel.find("RGB")!=string::npos) eInterp=GPI_RGB;
             else {
                 CPLError(CE_Failure, CPLE_AppDefined,
                     "GDAL MRF: Palette Model %s is unknown, use RGB,RGBA,HSV,HLS,CMYK or L",
                     pModel.c_str());
                 return CE_Failure;
             }
         }

         if ((entries>0)&&(entries<257)) {
             int start_idx, end_idx;
             GDALColorEntry ce_start={0,0,0,255},ce_end={0,0,0,255};

             // Create it and initialize it to nothing
             GDALColorTable *poColorTable = new GDALColorTable(eInterp);
             poColorTable->CreateColorRamp(0,&ce_start,entries-1,&ce_end);
             // Read the values
             CPLXMLNode *p=CPLGetXMLNode(node,"Entry");
             if (p) {
                 // Initialize the first entry, just in case
                 ce_start=GetXMLColorEntry(p);
                 if (is_hsv) ce_start=HSVSwap(ce_start);
                 start_idx=static_cast<int>(getXMLNum(p,"idx",0));
                 if (start_idx<0) {
                    CPLError(CE_Failure, CPLE_AppDefined,
                        "GDAL MRF: Palette index %d not allowed",start_idx);
                    delete poColorTable;
                    return CE_Failure;
                 }
                 poColorTable->SetColorEntry(start_idx,&ce_start);
                 while (NULL!=(p=SearchXMLSiblings(p,"Entry"))) {
                     // For every entry, create a ramp
                     ce_end=GetXMLColorEntry(p);
                     if (is_hsv) ce_end=HSVSwap(ce_end);
                     end_idx=static_cast<int>(getXMLNum(p,"idx",start_idx+1));
                     if ((end_idx<=start_idx)||(start_idx>=entries)) {
                         CPLError(CE_Failure, CPLE_AppDefined,
                             "GDAL MRF: Index Error at index %d",end_idx);
                         delete poColorTable;
                         return CE_Failure;
                     }
                     poColorTable->CreateColorRamp(start_idx,&ce_start,
                                                   end_idx,&ce_end);
                     ce_start=ce_end;
                     start_idx=end_idx;
                 }
             }