#include "vtkSystemIncludes.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkSTLReader.h"
#include "vtkPolyData.h"
#include "vtkPolygon.h"
#include "vtkPoints.h"
#include "vtkNormals.h"
#include "vtkTCoords.h"
#include "vtkCell.h"
#include "vtkPolyDataNormals.h"
#include "vtkStripper.h"
#include "vtkPolyDataMapper.h"
#include "vtkTexture.h"
#include "vtkPNMReader.h"
#include "vtkDecimatePro.h"
#include "vtkInteractorStyleFlight.h"

#include "paramFile.h"
#include "vector3d.h"
using osl::vector3d;


/*
Tile a rectangular grid of the given dimensions with
points with the given heights.  Write grid to given output.
*/
class grid {
	int nx,ny;
	float pixSize;
	const float *ht;
	int htRow; //Distance between rows of height data
	
	//Clipping height access routine
	float at(int x,int y) {
		if (x<0) x=0; if (x>=nx) x=nx-1;
		if (y<0) y=0; if (y>=ny) y=ny-1;
		return ht[y*htRow+x];
	}
	//Estimate surface normal at point (uses parabolic interpolation)
	vector3d normal(int x,int y) {
		float             u=at(x,y+1);
		float l=at(x-1,y),            r=at(x+1,y);
		float             d=at(x,y-1);
		vector3d gx(2*pixSize,0,r-l);
		vector3d gy(0,2*pixSize,u-d);
		return gx.cross(gy).dir();
	}
	
public:
	grid(int nx_,int ny_,float pixSize_,
		const float *ht_,int htRow_) 
		:nx(nx_), ny(ny_), pixSize(pixSize_),
		ht(ht_), htRow(htRow_) {}
	
	void build(vtkPolyData *output,bool withNormals,bool asQuads);

};

void grid::build(vtkPolyData *output,bool withNormals,bool asQuads)
{
  register int x,y;
  int numPts;
  int numPolys;
//Allocate
  numPts = nx * ny;
  numPolys = (nx-1)*(ny-1);
  vtkPoints *pts=vtkPoints::New();
  pts->Allocate(numPts); 
  vtkNormals *norm=NULL;
  if (withNormals) {
    norm=vtkNormals::New();
    norm->Allocate(numPts);
  }
  vtkTCoords *tc=vtkTCoords::New();
  tc->Allocate(numPts,2);
  vtkCellArray *poly=vtkCellArray::New();
  if (asQuads) 
	poly->Allocate(poly->EstimateSize(numPolys,4));
  else
	poly->Allocate(poly->EstimateSize(2*numPolys,3));
  
//Place points and texture coordinates
  double tDx=1.0/(nx);
  double tDy=1.0/(ny);
  float t[2];
  t[0]=t[1]=0.0;
  for (y=0;y<ny;y++) {
	const float *hRow=&ht[y*htRow];
	float pt[3]; pt[0]=0; pt[1]=y*pixSize;
//	t[0]=0; t[1]=y*tDy;
	for (x=0;x<nx;x++) {
		t[0]=x*tDx;
		t[1]=y*tDy;
		pt[0]+=pixSize;
		pt[2]=hRow[x];
		pts->InsertNextPoint(pt);
		tc->InsertNextTCoord(t);
		if (withNormals) {
			vector3d n=normal(x,y);
			float na[3];
			na[0]=n.x; na[1]=n.y; na[2]=n.z;
			norm->InsertNextNormal(na);
		}
	}
  }

//Connect triangles
  for (y=0;y<ny-1;y++) {
	for (x=0;x<nx-1;x++) {
		int tl=(y  )*nx+(x  );
		int tr=(y  )*nx+(x+1);
		int bl=(y+1)*nx+(x  );
		int br=(y+1)*nx+(x+1);
		if (asQuads) {//Output quadrilaterals
			int conn[4];
			conn[0] = tl;
			conn[1] = tr;
			conn[2] = br;
			conn[3] = bl;
			poly->InsertNextCell(4,conn);
		} else { //Output triangles
			int t[3], b[3];
			t[0]=tl; t[1]=tr;
			t[2]=bl; 
			         b[0]=tr;
		    b[2]=bl; b[1]=br;
			poly->InsertNextCell(3,t);
			poly->InsertNextCell(3,b);
		}
	}
  }


//Output
  output->SetPoints(pts);
  pts->Delete();

  if (withNormals) {
    output->GetPointData()->SetNormals(norm);
    norm->Delete();
  }

  output->GetPointData()->SetTCoords(tc);
  tc->Delete();

  output->SetPolys(poly);
  poly->Delete();
}

/*Read a just-the-bytes-ma'm little-endian short integers
DEM file of the given size from the given file.*/
const char *readFlatDEM(const char *fName,float *dest,int nx,int ny) 
{
	FILE *f=fopen(fName,"rb");
	if (f==NULL) return "Can't open input file";
	int rowBytes=2*nx;
	unsigned char *row=new unsigned char[rowBytes];
	int y;
	for (y=ny-1;y>=0;y--) //Read DEM bottom-up
	{
		if (rowBytes!=fread((void *)row,1,rowBytes,f)) break;
		for (int x=0;x<nx;x++)
		{
			int ht=row[2*x+0]+(row[2*x+1]<<8);
			dest[y*nx+x]=(float)ht;
		}
	}
	fclose(f);
	delete[] row;
	if (y!=-1) return "Unexpected EOF";
	else return NULL;
}

int main (int argc,char *argv[])
{	
  vtkRenderer *aren    = vtkRenderer::New();
  vtkRenderWindow *renWin  = vtkRenderWindow::New();
  renWin->AddRenderer(aren);
  vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
  iren->SetRenderWindow(renWin);
  aren->SetBackground(0.4,0.6,1.0); //A nice sky blue

  vtkPolyDataMapper *mapper = vtkPolyDataMapper::New();
  vtkPolyData *pd=vtkPolyData::New();
  
  const char *paramFileName=argv[1];
  bool forWrite=false;
  if (paramFileName==NULL) {paramFileName="demparam.txt"; /*forWrite=true;*/}
  parameterFile param(paramFileName,forWrite);

  int origSize=param("dem.origSize","Original size of DEM (pixels)",4096);
  double origSpacing=param("dem.origPix","Original pixel spacing (meters/pixel)",40.0);
  int divide=param("dem.downsize","Downsampling factor to apply to DEM",16);
  const char *demPattern=param("dem.name","Filename pattern for raw DEM file","final_%d.dem");
  bool useTexture=param("textureOn","Use texture mapping",true);
  const char *textureName=NULL;
  if (useTexture)
	textureName=param("textureImage","PPM Texture map name","final_4.ppm");
  bool doDecimate=param("doDecimate","Apply mesh decimation filter",true);

  float demSpacing=origSpacing*divide;
  int demSize=origSize/divide;
  char demNameBuf[1024];
  sprintf(demNameBuf,demPattern,divide);

  int nx=demSize;
  int ny=nx;
  float *ht=new float[nx*ny];
  const char *err;
  if (NULL!=(err=readFlatDEM(demNameBuf,ht,nx,ny))) {
	  fprintf(stderr,"ERROR> %s while reading file '%s'!\n",err,demNameBuf);
	  exit(1);
  }
  else
	  printf("Read %dx%d DEM from %s\n",nx,ny,demNameBuf);
  
  grid g(nx,ny,demSpacing,ht,nx);
  bool addNormals=true;
  bool asQuads=false;
  g.build(pd,addNormals,asQuads);

  int lookLocX=demSize*0.61, lookLocY=demSize*0.39;
  double lookPos[3]={lookLocX*demSpacing,lookLocY*demSpacing,ht[lookLocY*nx+lookLocX]};
  double altitude=1.0e3; //Initial height of user in meters above ground
  double camPos[3]={lookPos[0],lookPos[1]+20.0e3,lookPos[2]+altitude};
  double upVec[3]={0,0,1};
  
  if (false) { /*Fly around instead of rotating-- buggy!*/
    vtkInteractorStyleFlight *iflight=vtkInteractorStyleFlight::New();
    iflight->JumpTo(camPos,lookPos);
	iflight->SetFixUpVector(1);
	iflight->SetFixedUpVector(upVec);
    iren->SetInteractorStyle(iflight);
	iflight->Delete();
  } else {
	vtkCamera *cam=aren->GetActiveCamera();
	cam->SetFocalPoint(lookPos);
	cam->SetPosition(camPos);
	cam->SetViewUp(upVec);
	cam->SetClippingRange(10.0,1000.0e3);
  }

  delete ht;

  if (doDecimate) { /*Decimate model*/
	printf("Applying decimate filter...\n");
	vtkDecimatePro *dec=vtkDecimatePro::New();
	dec->SetTargetReduction(0.999);
	dec->SetPreserveTopology(1);
	dec->SetBoundaryVertexDeletion(1);
	dec->SetMaximumError(1.0e-8);
	dec->SetInput(pd);
	
	if (true) {/*Rebuild normals and strip triangles*/
		vtkPolyDataNormals *norm=vtkPolyDataNormals::New();
		norm->SetInput(dec->GetOutput());
		vtkStripper *strip=vtkStripper::New();
		strip->SetInput(norm->GetOutput());
		mapper->SetInput(strip->GetOutput());
		norm->Delete();
		strip->Delete();
	}
	dec->Delete();
  }
  else
	mapper->SetInput(pd);
  
  vtkActor *actor1 = vtkActor::New();
  actor1->SetMapper(mapper);
  aren->AddActor(actor1);

  if (useTexture) {
    printf("Reading texture image...\n");
    vtkTexture *atext=vtkTexture::New();
    vtkPNMReader *pnmReader=vtkPNMReader::New();
    pnmReader->SetFileName(textureName);
    atext->SetInput(pnmReader->GetOutput());
    atext->InterpolateOn();
    actor1->SetTexture(atext);
    actor1->GetProperty()->SetDiffuseColor(1,1,1);
    actor1->GetProperty()->SetDiffuse(0.3);
    actor1->GetProperty()->SetSpecular(0.0);
    actor1->GetProperty()->SetAmbientColor(1,1,1);
    actor1->GetProperty()->SetAmbient(0.7);
	atext->Delete();
  }

  printf("About to interact\n");
  renWin->Render();
  iren->Start();
  //renWin->Render();
  // interact with data

  // Clean up
  aren->Delete();
  renWin->Delete();
  iren->Delete();
  pd->Delete();
  mapper->Delete();
  actor1->Delete();
  return 0;
}
