#include "gdalwarper.h"
#include "cpl_conv.h"

int main(int argc, char **argv) {
	if (argc != 3) {
		printf("%d\n", argc);
		printf("%s source destination\n", argv[0]);
		exit(1);
	}
	GDALDatasetH  hSrcDS, hDstDS;

	// Open input and output files.

	GDALAllRegister();

	hSrcDS = GDALOpen( argv[1], GA_ReadOnly );
	hDstDS = GDALOpen( argv[2], GA_Update );

	// Setup warp options.


	GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();

	psWarpOptions->eResampleAlg = GRA_CubicSpline; 

	psWarpOptions->hSrcDS = hSrcDS;
	psWarpOptions->hDstDS = hDstDS;
	
	psWarpOptions->nBandCount = 3;
	psWarpOptions->panSrcBands =
		(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
	psWarpOptions->panSrcBands[0] = 1;
	psWarpOptions->panSrcBands[1] = 2;
	psWarpOptions->panSrcBands[2] = 3;
	psWarpOptions->panDstBands =
		(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
	psWarpOptions->panDstBands[0] = 1;
	psWarpOptions->panDstBands[1] = 2;
	psWarpOptions->panDstBands[2] = 3;

	psWarpOptions->pfnProgress = GDALTermProgress;

	// Establish reprojection transformer.

	psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(
		hSrcDS,
		GDALGetProjectionRef(hSrcDS),
		hDstDS,
		GDALGetProjectionRef(hDstDS),
		FALSE, 0.0, 1
	);

	psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

	// Initialize and execute the warp operation.

	GDALWarpOperation oOperation;

	oOperation.Initialize( psWarpOptions );
	oOperation.ChunkAndWarpImage( 
		0, 0, 
		GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS )
	);

	GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
	GDALDestroyWarpOptions( psWarpOptions );



	psWarpOptions = GDALCreateWarpOptions();

	psWarpOptions->eResampleAlg = GRA_NearestNeighbour;

	psWarpOptions->hSrcDS = hSrcDS;
	psWarpOptions->hDstDS = hDstDS;
	
	psWarpOptions->nBandCount = 1;
	psWarpOptions->panSrcBands =
		(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
	psWarpOptions->panSrcBands[0] = 4;
	psWarpOptions->panDstBands =
		(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
	psWarpOptions->panDstBands[0] = 4;

	psWarpOptions->pfnProgress = GDALTermProgress;

	// Establish reprojection transformer.

	psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(
		hSrcDS,
		GDALGetProjectionRef(hSrcDS),
		hDstDS,
		GDALGetProjectionRef(hDstDS),
		FALSE, 0.0, 1
	);

	psWarpOptions->pfnTransformer = GDALGenImgProjTransform;


	// Initialize and execute the warp operation.

	oOperation.Initialize( psWarpOptions );
	oOperation.ChunkAndWarpImage( 
		0, 0, 
		GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS )
	);


	
	GDALDestroyWarpOptions( psWarpOptions );
	GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );

	GDALClose( hDstDS );
	GDALClose( hSrcDS );

	return 0;
}