/* file:        GDALRasterizeBurnZTests.h
 * author:        Chris Hanson http://xenon.arcticus.com/ xenon@alphapixel.com 2009-10-30
 * copyright:        (C) 2009 Chris Hanson
 * license:        
*/

#include <iostream>
#include <assert.h>
#include <gdal_priv.h> // for GDALDataset
#include "gdal_alg.h" // for GDALRasterizeGeometries


class GDALRasterizeBurnZTestCase /* : public TestCase */ {
public:
	GDALRasterizeBurnZTestCase(std::string name) /* : TestCase(name) */ {}

	// method to GDAL MEM dataset construction, fill and internal buffer access
	void testGDALRasterizeBurnZ(void)
	{
		std::cout << std::endl << "testGDALRasterizeBurnZ" << std::endl << std::endl;
		char *textPtrRollinsville = "LINESTRING(-105.6278770976 39.9035832722 2783.114, -105.6292072836 39.9033961358 2785.117, -105.6307144906 39.9031837106 2786.884, -105.6313922280 39.9030876135 2787, -105.6318196071 39.9030522093 2787.464, -105.6322444574 39.9030572671 2788, -105.6325681529 39.9030876135 2788.2, -105.6330435806 39.9031533642 2789, -105.6334734886 39.9032444035 2790, -105.6341866301 39.9034113090 2791, -105.6345760762 39.9034972906 2793.729, -105.6348795407 39.9035326948 2794, -105.6353094488 39.9035529258 2794.108, -105.6357899342 39.9035276370 2796.598, -105.6361642071 39.9034871751 2797.162, -105.6379622342 39.9032696922 2798.4, -105.6391305725 39.9031280755 2799.362, -105.6419300324 39.9027487449 2806, -105.6438115122 39.9024908000 2807.401, -105.6448660513 39.9023643565 2810.035, -105.6459787545 39.9022025088 2814.337, -105.6459787545 39.9022025088 2814.337)";
		OGRLineString lineStringOne;
		GDALDataset *_workingElevationGrid;
		OGRFeature *featureOne = NULL;

		// initialize GDAL (only needs to happen once, but apparently no harm in doing it extra times)
		GDALAllRegister();

		// Rollinsville Test code 
		if(OGRFeatureDefn *featureDefn = new OGRFeatureDefn("SimplePolygon"))
		{
			featureDefn->SetGeomType(wkbLineString);
			if(featureOne = OGRFeature::CreateFeature(featureDefn))
			{
				char *textPtrRollinsvilleTemp = textPtrRollinsville; // GDAL rewrites the textPtrRollinsvilleTemp pointer after we pass it, so we need a temp copy
				osg::notify(osg::DEBUG_INFO) << "ReaderWriterMODIFYTERRAIN WKT: " << textPtrRollinsvilleTemp << std::endl;
				OGRErr status = lineStringOne.importFromWkt(&textPtrRollinsvilleTemp);
				assert(status == OGRERR_NONE);
				featureOne->SetGeometry( &lineStringOne ); 

			} // if
		} // if

		// create sub-grid using GDAL
		// inspired by example code found here:
		// http://article.gmane.org/gmane.comp.gis.gdal.devel/14908
		// Acquire access to in-memory GDAL dataset driver. We are apparently
		// not required to release this resource explicitly, nor that of the 
		// GDALDriverManager singleton
		if(GDALDriver* memDriver = GetGDALDriverManager()->GetDriverByName("MEM"))
		{
			// use MEM driver to create a new dataset to our specifications
			if(_workingElevationGrid = memDriver->Create("testGDALMEMaccess_dataset", 2825, 354, 1 /* nBands */, GDT_Float32, NULL /* paramlist */)) // dimensions found empirically
			{
				double padfTransform[6];
				padfTransform[2] = padfTransform[4] = 0.0; // clear unused terms
				padfTransform[0] = -105.64667287591845; // The upper left corner of the upper left pixel (X)
				padfTransform[3] = 39.904474122037868; // The upper left corner of the upper left pixel (Y)
				padfTransform[1] = 6.8991686501776537e-006; // padfTransform[1] is the pixel width,
				padfTransform[5] = -8.9932101263546030e-006; // padfTransform[5] is the pixel height (usually negative)
				_workingElevationGrid->SetGeoTransform(padfTransform); // set geospatial-to-raster affine transformation coefficients
				// set up projection as WGS84 geographic. Not 100% sure if this is necessary for our usage.
				_workingElevationGrid->SetProjection(WKTGeographicWGS84.c_str());

				// initialize dataset to known NODATA value
				GDALRasterBand *primaryBand = NULL;
				if(primaryBand = _workingElevationGrid->GetRasterBand(1)) // GDAL counts bands Pascal-style from 1... It does not appear you need to release the GDALRasterBand object
				{
					primaryBand->SetNoDataValue(0.0);
					primaryBand->Fill(0.0);

					std::vector<int> anBandList;
					std::vector<double> adfFullBurnValues;
					std::vector<OGRGeometryH> ahGeometries;
					char **papszRasterizeOptions = NULL;

					papszRasterizeOptions = CSLSetNameValue( papszRasterizeOptions, "BURN_VALUE_FROM", "Z" );
					
					anBandList.push_back(1); // operate on band 1. GDAL counts bands Pascal-style, from 1
					adfFullBurnValues.resize(1, 10.0); // setup elevation offset value (10.0) to burn in for this geometry
					
					// insert the Feature's Geometry into Geometries container
					ahGeometries.push_back(featureOne->GetGeometryRef());

					AbstractMapping::notify(osg::DEBUG_INFO) << "DeformerPolyline GDALRasterizeGeometries" << std::endl;
					AbstractMapping::notify(osg::DEBUG_INFO) << "_workingElevationGrid size " << _workingElevationGrid->GetRasterXSize() << ", " << _workingElevationGrid->GetRasterXSize() << std::endl;

					// call GDALRasterizeGeometries to do the work of rasterizing the polygon
					if(GDALRasterizeGeometries(_workingElevationGrid, anBandList.size(), &(anBandList[0]),
						ahGeometries.size(), &(ahGeometries[0]), NULL /* transform */, NULL /* transform data */,
						&(adfFullBurnValues[0]), papszRasterizeOptions, GDALDummyProgress, NULL /* progress arg */) == CE_None)
					{
						AbstractMapping::notify(osg::DEBUG_INFO) << "GDALRasterizeGeometries done, " << ahGeometries.size() << " geometries, " << anBandList.size() << " bands." << std::endl;

							
						// read and validate some values
						float result32 = 0.0;
						int panBandMap = 1; // operate on band 1. GDAL counts bands Pascal-style, from 1

						// read data from left end (101, 253 in GDAL's top-left 0,0 notation)
						_workingElevationGrid->RasterIO(GF_Read, 101, 253, 1 /* pixels */, 1 /* pixels */, &result32,
							1 /* destw */, 1 /* desth */, GDT_Float32, 1 /* bands */, &panBandMap ,
							0 /* nPixelSpace */, 0 /* nLineSpace */, 0 /* nBandSpace */);
						assert(fabs(result32 - 2824.337) < .005); // 2824.337 = 2814.337 + 10m offset

						// read data from right end (2724, 99 in GDAL's top-left 0,0 notation)
						_workingElevationGrid->RasterIO(GF_Read, 2724, 99, 1 /* pixels */, 1 /* pixels */, &result32,
							1 /* destw */, 1 /* desth */, GDT_Float32, 1 /* bands */, &panBandMap ,
							0 /* nPixelSpace */, 0 /* nLineSpace */, 0 /* nBandSpace */);
						assert(fabs(result32 - 2793.114) < .005); // 2793.114 = 2783.114 + 10m offset

						// read data from middle (1636, 102 in GDAL's top-left 0,0 notation)
						_workingElevationGrid->RasterIO(GF_Read, 1636, 102, 1 /* pixels */, 1 /* pixels */, &result32,
							1 /* destw */, 1 /* desth */, GDT_Float32, 1 /* bands */, &panBandMap ,
							0 /* nPixelSpace */, 0 /* nLineSpace */, 0 /* nBandSpace */);
						assert(fabs(result32 - 2804.499) < .005); // 2804.499 = 2794.499 + 10m offset


					} // if


				} // if
				assert(primaryBand);

			} // if Dataset created ok
			assert(_workingElevationGrid);

			// clean up
			delete _workingElevationGrid; _workingElevationGrid = NULL;
		} // if MEM driver available


	} // testGDALRasterizeBurnZ


	// method to create a suite of tests
	/* static Test *suite (); */
}; // GDALRasterizeBurnZTestCase

/*
// method to create a suite of tests
Test *DeformTestCase::suite () {
  TestSuite *testSuite = new TestSuite ("GDALRasterizeBurnZTestCase");
  
  // add the tests
  testSuite->addTest (new TestCaller <GDALRasterizeBurnZTestCase> ("testGDALRasterizeBurnZ", &GDALRasterizeBurnZTestCase::testGDALRasterizeBurnZ));
  return testSuite;
}
*/