GDAL映像統合
15599 ワード
// merge.cpp : 。
//
#include "stdafx.h"
#include "../3dparty/include/gdal_priv.h"
#include "../3dparty/include/gdal_alg.h"
#include "../3dparty/include/gdalwarper.h"
#include
#ifdef _DEBUG
#pragma comment(lib, "../3dparty/lib/gdal_i.lib")
#else
#pragma comment(lib, "../3dparty/lib/gdal_i.lib")
#endif
#define SIZEOFARRAY(array) sizeof(array)/sizeof(array[0])
static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
static double dfXRes=0.0, dfYRes=0.0;
static int bTargetAlignedPixels = FALSE;
static int nForcePixels=0, nForceLines=0, bQuiet = FALSE;
static int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
static int bVRT = FALSE;
static char *pszDstFileName = "../data/out.tif";
static char *apszSrcFiles[] =
{
/*"../data/H50G038013DOM.tif",
"../data/H50G038014DOM.tif",
"../data/H50G038015DOM.tif",
"../data/H50G039012DOM.tif",
"../data/H50G039013DOM.tif",
"../data/H50G039014DOM.tif",*/
"../data/H50G039015DOM.tif",
"../data/H50G039016DOM.tif"
};
GDALDataType eWorkingType = GDT_Byte;
static GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
static GDALDatasetH GDALWarpCreateOutput( char **papszSrcFiles,
int nFilesCount,
const char *pszFilename,
const char *pszFormat,
GDALDataType eDT );
static void CreateMask(char* const* paFilePath, const int nFilesCount);
int _tmain(int argc, _TCHAR* argv[])
{
GDALAllRegister();
GDALDatasetH hSrcDS = NULL;
GDALDatasetH hDstDS = NULL;
int nFilesCount = SIZEOFARRAY(apszSrcFiles);
hDstDS = GDALWarpCreateOutput(apszSrcFiles, nFilesCount,
pszDstFileName, "GTiff", GDT_Byte);
char **papszSrcFiles = apszSrcFiles;
CreateMask(papszSrcFiles, nFilesCount);
for( int iSrc = 0; iSrc < nFilesCount; iSrc++ )
{
GDALDatasetH hSrcDS;
hSrcDS = GDALOpen( apszSrcFiles[iSrc], GA_ReadOnly );
GDALDataset *pSrcDS = (GDALDataset *)hSrcDS;
assert(pSrcDS->GetRasterBand(1)->GetMaskBand() != NULL);
if( hSrcDS == NULL )
exit( 2 );
if ( GDALGetRasterCount(hSrcDS) == 0 )
{
fprintf(stderr, "Input file %s has no raster bands.
", apszSrcFiles[iSrc] );
exit( 1 );
}
if( !bQuiet )
printf( "Processing input file %s.
", apszSrcFiles[iSrc] );
if ( eResampleAlg != GRA_NearestNeighbour &&
GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
{
if( !bQuiet )
fprintf( stderr, "Warning: Input file %s has a color table, which will likely lead to "
"bad results when using a resampling method other than "
"nearest neighbour. Converting the dataset prior to 24/32 bit "
"is advised.
", apszSrcFiles[iSrc] );
}
if( GDALGetRasterColorInterpretation(
GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) )
== GCI_AlphaBand
&& !bEnableSrcAlpha )
{
bEnableSrcAlpha = TRUE;
if( !bQuiet )
printf( "Using band %d of source image as alpha.
",
GDALGetRasterCount(hSrcDS) );
}
void *hTransformArg = NULL;
hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, NULL );
if( hTransformArg == NULL )
exit( 1 );
GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform;
GDALWarpOptions *psWO = GDALCreateWarpOptions();
psWO->papszWarpOptions = NULL;
psWO->eWorkingDataType = eWorkingType;
psWO->eResampleAlg = eResampleAlg;
psWO->hSrcDS = hSrcDS;
psWO->hDstDS = hDstDS;
psWO->pfnTransformer = pfnTransformer;
psWO->pTransformerArg = hTransformArg;
if( !bQuiet )
psWO->pfnProgress = GDALTermProgress;
if( bEnableSrcAlpha )
psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
else
psWO->nBandCount = GDALGetRasterCount(hSrcDS);
psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
for( int nBandIndex = 0; nBandIndex < psWO->nBandCount; nBandIndex++ )
{
psWO->panSrcBands[nBandIndex] = nBandIndex+1;
psWO->panDstBands[nBandIndex] = nBandIndex+1;
}
if( bEnableSrcAlpha )
psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);
if( !bEnableDstAlpha
&& GDALGetRasterCount(hDstDS) == psWO->nBandCount+1
&& GDALGetRasterColorInterpretation(
GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS)))
== GCI_AlphaBand )
{
if( !bQuiet )
printf( "Using band %d of destination image as alpha.
",
GDALGetRasterCount(hDstDS) );
bEnableDstAlpha = TRUE;
}
GDALWarpOperation oWO;
if( oWO.Initialize( psWO ) == CE_None )
{
CPLErr eErr;
eErr = oWO.ChunkAndWarpImage( 0, 0,
GDALGetRasterXSize( hDstDS ),
GDALGetRasterYSize( hDstDS ) );
}
GDALDestroyGenImgProjTransformer( hTransformArg );
GDALDestroyWarpOptions( psWO );
GDALClose( hSrcDS );
}
GDALClose( hDstDS );
return 0;
}
static GDALDatasetH
GDALWarpCreateOutput( char **papszSrcFiles, int nFilesCount,
const char *pszFilename, const char *pszFormat,
GDALDataType eDT )
{
GDALDriverH hDriver;
GDALDatasetH hDstDS;
void *hTransformArg;
GDALColorTableH hCT = NULL;
double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0;
double dfWrkResX=0, dfWrkResY=0;
int nDstBandCount = 0;
hDriver = GDALGetDriverByName( pszFormat );
if( hDriver == NULL
|| GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
{
int iDr;
printf( "Output driver `%s' not recognised or does not support
",
pszFormat );
printf( "direct output file creation. The following format drivers are configured
"
"and support direct output:
" );
for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
{
GDALDriverH hDriver = GDALGetDriver(iDr);
if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
{
printf( " %s: %s
",
GDALGetDriverShortName( hDriver ),
GDALGetDriverLongName( hDriver ) );
}
}
printf( "
" );
exit( 1 );
}
int iSrc = 0;
for( iSrc = 0; iSrc < nFilesCount; iSrc++ )
{
GDALDatasetH hSrcDS;
hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
if( hSrcDS == NULL )
exit( 1 );
if ( GDALGetRasterCount(hSrcDS) == 0 )
{
fprintf(stderr, "Input file %s has no raster bands.
", papszSrcFiles[iSrc] );
exit( 1 );
}
if( eDT == GDT_Unknown )
eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
if( iSrc == 0 )
{
nDstBandCount = GDALGetRasterCount(hSrcDS);
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
if( hCT != NULL )
{
hCT = GDALCloneColorTable( hCT );
}
}
hTransformArg =
GDALCreateGenImgProjTransformer2( hSrcDS, NULL, NULL);
if( hTransformArg == NULL )
{
GDALClose( hSrcDS );
return NULL;
}
double adfThisGeoTransform[6];
double adfExtent[4];
int nThisPixels, nThisLines;
if( GDALSuggestedWarpOutput2( hSrcDS,
GDALGenImgProjTransform, hTransformArg,
adfThisGeoTransform,
&nThisPixels, &nThisLines,
adfExtent, 0 ) != CE_None )
{
GDALClose( hSrcDS );
return NULL;
}
if (CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", NULL ) == NULL)
{
double MinX = adfExtent[0];
double MaxX = adfExtent[2];
double MaxY = adfExtent[3];
double MinY = adfExtent[1];
int bSuccess = TRUE;
#define N_STEPS 20
int i = 0;
int j = 0;
for(i=0; i <= N_STEPS && bSuccess; i++)
{
for(j=0; j <= N_STEPS && bSuccess; j++)
{
double dfRatioI = i * 1.0 / N_STEPS;
double dfRatioJ = j * 1.0 / N_STEPS;
double expected_x = (1 - dfRatioI) * MinX + dfRatioI * MaxX;
double expected_y = (1 - dfRatioJ) * MinY + dfRatioJ * MaxY;
double x = expected_x;
double y = expected_y;
double z = 0;
if (!GDALGenImgProjTransform(hTransformArg, TRUE, 1, &x, &y, &z, &bSuccess)
|| !bSuccess)
bSuccess = FALSE;
if (!GDALGenImgProjTransform(hTransformArg, FALSE, 1, &x, &y, &z, &bSuccess)
|| !bSuccess)
bSuccess = FALSE;
if (fabs(x - expected_x) > (MaxX - MinX) / nThisPixels ||
fabs(y - expected_y) > (MaxY - MinY) / nThisLines)
bSuccess = FALSE;
}
}
if (!bSuccess)
{
CPLSetConfigOption( "CHECK_WITH_INVERT_PROJ", "TRUE" );
CPLDebug("WARP", "Recompute out extent with CHECK_WITH_INVERT_PROJ=TRUE");
GDALDestroyGenImgProjTransformer(hTransformArg);
hTransformArg =
GDALCreateGenImgProjTransformer2( hSrcDS, NULL, NULL );
if( GDALSuggestedWarpOutput2( hSrcDS,
GDALGenImgProjTransform, hTransformArg,
adfThisGeoTransform,
&nThisPixels, &nThisLines,
adfExtent, 0 ) != CE_None )
{
GDALClose( hSrcDS );
return NULL;
}
}
}
if( dfWrkMaxX == 0.0 && dfWrkMinX == 0.0 )
{
dfWrkMinX = adfExtent[0];
dfWrkMaxX = adfExtent[2];
dfWrkMaxY = adfExtent[3];
dfWrkMinY = adfExtent[1];
dfWrkResX = adfThisGeoTransform[1];
dfWrkResY = ABS(adfThisGeoTransform[5]);
}
else
{
dfWrkMinX = MIN(dfWrkMinX,adfExtent[0]);
dfWrkMaxX = MAX(dfWrkMaxX,adfExtent[2]);
dfWrkMaxY = MAX(dfWrkMaxY,adfExtent[3]);
dfWrkMinY = MIN(dfWrkMinY,adfExtent[1]);
dfWrkResX = MIN(dfWrkResX,adfThisGeoTransform[1]);
dfWrkResY = MIN(dfWrkResY,ABS(adfThisGeoTransform[5]));
}
GDALDestroyGenImgProjTransformer( hTransformArg );
GDALClose( hSrcDS );
}
if( nDstBandCount == 0 )
{
CPLError( CE_Failure, CPLE_AppDefined,
"No usable source images." );
return NULL;
}
double adfDstGeoTransform[6];
int nPixels, nLines;
adfDstGeoTransform[0] = dfWrkMinX;
adfDstGeoTransform[1] = dfWrkResX;
adfDstGeoTransform[2] = 0.0;
adfDstGeoTransform[3] = dfWrkMaxY;
adfDstGeoTransform[4] = 0.0;
adfDstGeoTransform[5] = -1 * dfWrkResY;
nPixels = (int) ((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5);
nLines = (int) ((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5);
if( dfXRes != 0.0 && dfYRes != 0.0 )
{
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
{
dfMinX = adfDstGeoTransform[0];
dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
dfMaxY = adfDstGeoTransform[3];
dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
}
if ( bTargetAlignedPixels )
{
dfMinX = floor(dfMinX / dfXRes) * dfXRes;
dfMaxX = ceil(dfMaxX / dfXRes) * dfXRes;
dfMinY = floor(dfMinY / dfYRes) * dfYRes;
dfMaxY = ceil(dfMaxY / dfYRes) * dfYRes;
}
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
adfDstGeoTransform[1] = dfXRes;
adfDstGeoTransform[5] = -dfYRes;
}
else if( nForcePixels != 0 && nForceLines != 0 )
{
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
{
dfMinX = dfWrkMinX;
dfMaxX = dfWrkMaxX;
dfMaxY = dfWrkMaxY;
dfMinY = dfWrkMinY;
}
dfXRes = (dfMaxX - dfMinX) / nForcePixels;
dfYRes = (dfMaxY - dfMinY) / nForceLines;
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
adfDstGeoTransform[1] = dfXRes;
adfDstGeoTransform[5] = -dfYRes;
nPixels = nForcePixels;
nLines = nForceLines;
}
else if( nForcePixels != 0 )
{
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
{
dfMinX = dfWrkMinX;
dfMaxX = dfWrkMaxX;
dfMaxY = dfWrkMaxY;
dfMinY = dfWrkMinY;
}
dfXRes = (dfMaxX - dfMinX) / nForcePixels;
dfYRes = dfXRes;
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
adfDstGeoTransform[1] = dfXRes;
adfDstGeoTransform[5] = -dfYRes;
nPixels = nForcePixels;
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
}
else if( nForceLines != 0 )
{
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
{
dfMinX = dfWrkMinX;
dfMaxX = dfWrkMaxX;
dfMaxY = dfWrkMaxY;
dfMinY = dfWrkMinY;
}
dfYRes = (dfMaxY - dfMinY) / nForceLines;
dfXRes = dfYRes;
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
adfDstGeoTransform[1] = dfXRes;
adfDstGeoTransform[5] = -dfYRes;
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
nLines = nForceLines;
}
else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
{
dfXRes = adfDstGeoTransform[1];
dfYRes = fabs(adfDstGeoTransform[5]);
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
dfXRes = (dfMaxX - dfMinX) / nPixels;
dfYRes = (dfMaxY - dfMinY) / nLines;
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
adfDstGeoTransform[1] = dfXRes;
adfDstGeoTransform[5] = -dfYRes;
}
if( bEnableSrcAlpha )
nDstBandCount--;
if( bEnableDstAlpha )
nDstBandCount++;
if( !bQuiet )
printf( "Creating output file that is %dP x %dL.
", nPixels, nLines );
hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
nDstBandCount, eDT, NULL );
if( hDstDS == NULL )
{
return NULL;
}
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
if( bEnableDstAlpha )
{
GDALSetRasterColorInterpretation(
GDALGetRasterBand( hDstDS, nDstBandCount ),
GCI_AlphaBand );
}
if( hCT != NULL )
{
GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
GDALDestroyColorTable( hCT );
}
return hDstDS;
}
void CreateMask(char* const* paFilePath, const int nFilesCount)
{
for ( int nSrcIndex = 0; nSrcIndex < nFilesCount; nSrcIndex++ )
{
GDALDatasetH hSrcDS;
hSrcDS = GDALOpen( apszSrcFiles[nSrcIndex], GA_Update );
GDALDataset *pSrcDS = (GDALDataset *)hSrcDS;
pSrcDS->CreateMaskBand(GMF_ALL_VALID);
for (int nBandIndex = 1; nBandIndex <= pSrcDS->GetRasterCount(); ++nBandIndex)
{
GDALRasterBand *pSrcBand = pSrcDS->GetRasterBand(nBandIndex);
GDALRasterBand *pMaskBand = pSrcBand->GetMaskBand();
double dfNoDataValue = pMaskBand->GetNoDataValue();
dfNoDataValue = pSrcBand->GetNoDataValue();
int nMaskXSize = pMaskBand->GetXSize();
int nMaskYSize = pMaskBand->GetYSize();
unsigned char *pMaskData = new unsigned char[nMaskXSize * nMaskYSize];
unsigned int *pMaskDataInt = (unsigned int *)pMaskData;
pMaskBand->RasterIO(GF_Read,0, 0, nMaskXSize, nMaskYSize,
pMaskData, nMaskXSize, nMaskYSize, GDT_Byte, 0, 0);
unsigned char *pRasterData = new unsigned char[pSrcBand->GetXSize() *
pSrcBand->GetYSize()];
pSrcBand->RasterIO(GF_Read,0, 0, pSrcBand->GetXSize(), pSrcBand->GetYSize(),
pRasterData, pSrcBand->GetXSize(), pSrcBand->GetYSize(), GDT_Byte, 0, 0);
for ( int nMaskPos = 0; nMaskPos < nMaskXSize * nMaskYSize; ++nMaskPos )
{
if (pRasterData[nMaskPos] == 0)
{
pMaskData[nMaskPos] = 0x00;
}
else
{
pMaskData[nMaskPos] = 0xFF;
}
}
pMaskBand->RasterIO(GF_Write, 0, 0, pMaskBand->GetXSize(),
pMaskBand->GetYSize(), pMaskData, nMaskXSize,
nMaskYSize, GDT_Byte, 0, 0);
delete []pRasterData;
pRasterData = NULL;
delete []pMaskData;
pMaskData = NULL;
}
GDALClose(hSrcDS);
}
}
マスクを自分で作る.