/*=========================================================================

This script is written to visualize the point cloud data (PCD) generated from LiDAR system,
and provide the result of voxelization process for a given PCD by defining the dimension characteristics of each voxel.  It will serve as a basis of point cloud slicing algorithm (PCS) 
in order to retrieve the forest biophysical information such as (forest inventory parameters, leaf area index, etc.)

Guang Zheng
guangz@u.washington.edu

Remote Sensing and Geospatial Analysis Laboratory (RSGAL)
School of Forest Resources, College of Environment
University of Washington

=========================================================================*/





#include <iostream>
#include <fstream>
#include <iomanip>
#include <algorithm>
#include <vector>

#include "vtkPoints.h"
#include "vtkPointData.h"
#include "vtkDoubleArray.h"
#include "vtkPolyData.h"
#include "vtkPolyVertex.h"
#include "vtkThresholdPoints.h"
#include "vtkTransform.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkPolyDatamapper.h"
#include "vtkPolyDataNormals.h"
#include "vtkLODActor.h"
#include "vtkScalarBarActor.h"
#include "vtkOutlineFilter.h"
#include "vtkActor.h"
#include "vtkTextProperty.h"
#include "vtkCubeAxesActor2D.h"
#include "vtkCamera.h"
#include "vtkPlaneSource.h"
#include "vtkFollower.h"
#include "vtkTextActor.h"
#include "vtkImageData.h"
#include "vtkExtractEdges.h"
#include "vtkTubeFilter.h"
#include "vtkDataSetMapper.h"
#include "vtkSphereSource.h"
#include "vtkThresholdPoints.h"
#include "vtkGlyph3D.h"
#include "vtkAnnotatedCubeActor.h"
#include "vtkAxesActor.h"
#include "vtkPropAssembly.h"
#include "vtkOrientationMarkerWidget.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkInteractorStyleTrackballCamera.h"
#include "vtkProperty.h"

using namespace std;




int main()
{

	

	cout <<"//---------------------------------------------------------//"<<endl;
	cout <<endl<<"This script is written to visualize the point cloud data (PCD) "<<endl;
	cout <<"generated from LiDAR system,and provide the result of voxelization"<<endl;
	cout <<"process for a given PCD by defining the dimension,characteristics"<<endl;
	cout <<"of each voxel.  It will serve as a basis of point cloud slicing "<<endl;
	cout <<" algorithm (PCS)in order to retrieve the forest biophysical"<<endl;
	cout <<"information such as (forest inventory parameters, leaf area index, etc.)"<<endl;
	cout <<endl;
	cout <<"Guang Zheng guangz@u.washington.edu"<<endl;
	cout <<"Remote Sensing and Geospatial Analysis Laboratory (RSGAL)"<<endl;
	cout <<"School of Forest Resources, College of Environment, University of Washington"<<endl;
	cout <<endl<<"//--------------------------------------------------------//"<<endl;


/***************** Section one    *.txt file reading***************/



cout <<endl<<endl<<"Please specify the input path and file: "<<endl;
	char file[100];

	cin>>file;


//read the point cloud data from a *.txt file
//ifstream file_to_read("twotree.txt");
ifstream file_to_read(file);

//defining the maximum number of char in a line and the header line number
int max_num_of_char_in_a_line=1000,num_of_header_lines=0;

//skip the header characterline
for (int i=0; i<num_of_header_lines; ++i)
    file_to_read.ignore(max_num_of_char_in_a_line,'\n');
	
//create the separate vector data structure for each point
vector<double> xcoord_array, ycoord_array, zcoord_array;

//report the error message when reading file is wrong
if (!file_to_read)
	cerr<<"File could not be opened!"<<endl;

cout <<endl<<endl<<"reading file..."<<endl;
	
//read the x,y,z coordinates into three individual vector structure
double x,y,z;
while(!file_to_read.eof())
{
file_to_read>>x;
file_to_read.ignore(1);
file_to_read>>y;
file_to_read.ignore(1);
file_to_read>>z;
file_to_read.ignore(1);


xcoord_array.push_back(x);
ycoord_array.push_back(y);
zcoord_array.push_back(z);
};

/***************** Section one    *.txt file reading ends here***************/


cout <<endl<<endl<<"Visualizing points..."<<endl;

/***************** Section two    prepare the points for VTK pipeline start here***************/

vtkPoints *pt=vtkPoints::New();
    for (int i=0; i<xcoord_array.size(); i++)
	pt->InsertPoint(i,xcoord_array[i],ycoord_array[i],zcoord_array[i]);

int num=pt->GetNumberOfPoints();

	
vtkPolyVertex *vertex=vtkPolyVertex::New();
	vertex->GetPointIds()->SetNumberOfIds(num);
	for(int i=0; i<xcoord_array.size(); ++i)
	vertex->GetPointIds()->SetId(i,i);
	
vtkDoubleArray *scalars=vtkDoubleArray::New();
	for (int i=0; i<zcoord_array.size(); ++i)
	scalars->InsertTuple1(i,zcoord_array[i]);
	
//vtkPolyData *polydata=vtkPolyData::New();
//	polydata->Allocate(1,1);
//	polydata->InsertNextCell(vertex->GetCellType(), vertex->GetPointIds());
//	polydata->SetPoints(pt);
//	polydata->GetPointData()->SetScalars(scalars);
//	polydata->Update();
	vtkPolyData* polydata = vtkPolyData::New();
	polydata->Allocate(1, 1);
	polydata->InsertNextCell(vertex->GetCellType(), vertex->GetPointIds());
	polydata->SetPoints(pt);
	polydata->GetPointData()->SetScalars(scalars);

	// Removed the Update() call

	
vtkThresholdPoints *threshold=vtkThresholdPoints::New();
	threshold->SetInputData(polydata);
	threshold->ThresholdByUpper(-100);
	
	
vtkTransform *transform =vtkTransform::New();
	transform->Identity();
	transform->RotateY(0);
	transform->Scale(1,1,1);

vtkTransformPolyDataFilter *transformpoly=vtkTransformPolyDataFilter::New();
	transformpoly->SetTransform(transform);
	transformpoly->SetInputConnection(threshold->GetOutputPort());

vtkPolyDataNormals *normals=vtkPolyDataNormals::New();
	normals->SetInputConnection(transformpoly->GetOutputPort());
	
vtkPolyDataMapper *pointsmapper=vtkPolyDataMapper::New();
	pointsmapper->SetInputConnection(normals->GetOutputPort());
	pointsmapper->SetScalarRange(scalars->GetRange());
	
vtkLODActor *pointsactor=vtkLODActor::New();
	pointsactor->SetMapper(pointsmapper);



	
//create the scalarbar
vtkScalarBarActor *scalarbar=vtkScalarBarActor::New();
	scalarbar->SetLookupTable(pointsmapper->GetLookupTable());
	scalarbar->SetTitle("Tree Height(m)");
	scalarbar->GetPositionCoordinate()->SetCoordinateSystemToNormalizedViewport();
	scalarbar->GetActualPositionCoordinate()->SetValue(0.04,0.2);
	scalarbar->SetOrientationToVertical();
	scalarbar->SetWidth(0.1);
	scalarbar->SetHeight(0.8);
	
	
/***************** Section two    prepare the points for VTK pipeline ends here***************/
	

	
/***************** Section three    outline with scalar of point cloud data starts here***************/

vtkOutlineFilter *outline =vtkOutlineFilter::New();
	outline->SetInputConnection(normals->GetOutputPort());
	
vtkPolyDataMapper *mapOutline=vtkPolyDataMapper::New();
	mapOutline->SetInputConnection(outline->GetOutputPort());

vtkActor *outlineactor=vtkActor::New();
	outlineactor->SetMapper(mapOutline);
	outlineactor->GetProperty()->SetColor(0,0,0);

vtkTextProperty *tprop=vtkTextProperty::New();
	tprop->SetColor(1,1,1);
	tprop->ShadowOn();


	

/***************** Section three    outline with scalar of point cloud dataset ends here***************/

	
	
/***************** Section four    bottom plane for scene starts here***************/

vtkPlaneSource *plane=vtkPlaneSource::New();
	plane->SetXResolution(800);
	plane->SetYResolution(800);
	plane->SetPoint1(-1000,1000,0);
	plane->SetPoint2(1000,1000,0);
	plane->SetCenter(0,0,0);
	
vtkPolyDataMapper *planemapper=vtkPolyDataMapper::New();
	planemapper->SetInputConnection(plane->GetOutputPort());
	
vtkFollower *planeActor=vtkFollower::New();
	planeActor->SetMapper(planemapper);
	
vtkActor *planeactor=vtkActor::New();
	planeactor->SetMapper(planemapper);
	planeactor->GetProperty()->SetColor(0,1,0);
	planeactor->GetProperty()->SetOpacity(0.8);
	planeactor->GetProperty()->SetRepresentationToWireframe();
	
/***************** Section four    bottome plane for scene starts here***************/	



/***************** Section five    output the total number on the screen***************/

vtkTextActor *textactor=vtkTextActor::New();
	textactor->SetDisplayPosition(1,1);
	textactor->SetInput("The total number of points is:  ");
	
/***************** Section five    output the total number on the screen***************/
	


/***************** Section six    sort the order of the point dataset***************/

double len=xcoord_array.size();

double* x_first=&(xcoord_array[0]);
double* x_last=&(xcoord_array[len-1]);
double* x_max=std::max_element(x_first, x_last);
double* x_min=std::min_element(x_first, x_last);

double* y_first=&(ycoord_array[0]);
double* y_last=&(ycoord_array[len-1]);
double* y_max=std::max_element(y_first, y_last);
double* y_min=std::min_element(y_first, y_last);

double* z_first=&(zcoord_array[0]);
double* z_last=&(zcoord_array[len-1]);
double* z_max=std::max_element(z_first, z_last);
double* z_min=std::min_element(z_first, z_last);

cout <<"The total number of points in this file is:  " <<len<<endl;

cout <<"Minimum value of X is " <<*x_min<<endl;
cout <<"Maximum value of X is " <<*x_max<<endl;
cout <<"Minimum value of Y is " <<*y_min<<endl;
cout <<"Maximum value of Y is " <<*y_max<<endl;
cout <<"Minimum value of Z is " <<*z_min<<endl;
cout <<"Maximum value of Z is " <<*z_max<<endl;

/***************** Section six    sort the order of the point dataset***************/



/***************** Section seven  voxelization process starts******************************/


double x_width=10;
double y_length=10;
double z_height=10;

//int dims[3]={((*x_max)-(*x_min))/x_width+2,((*y_max)-(*y_min))/y_length+2,((*z_max)-(*z_min))/z_height+2};

int dims[3]={x_width+1,y_length+1,z_height+1};

vtkDoubleArray *voxelscalar=vtkDoubleArray::New();
	for (int i=0; i<dims[0]*dims[1]*dims[2]; ++i)
		voxelscalar->InsertNextValue(1.0);
	

vtkDoubleArray *Scalars=vtkDoubleArray::New();
	for (int i=0; i<dims[0]*dims[1]*dims[2]; ++i)
		Scalars->InsertNextValue(1.0);
	
	
vtkImageData *boundingvolume =vtkImageData::New();
	boundingvolume->SetDimensions(dims);
	boundingvolume->SetOrigin((*x_min), (*y_min), (*z_min));
	double sp_x=((*x_max)-(*x_min))/10;
	double sp_y=((*y_max)-(*y_min))/10;
	double sp_z=((*z_max)-(*z_min))/10;
	boundingvolume->SetSpacing(sp_x,sp_y,sp_z);
	boundingvolume->GetPointData()->SetScalars(Scalars);

	
vtkExtractEdges *voxeledge=vtkExtractEdges::New();
	voxeledge->SetInputData(boundingvolume);

vtkTubeFilter *voxeledgetube=vtkTubeFilter::New();
	voxeledgetube->SetInputConnection(voxeledge->GetOutputPort());
	voxeledgetube->SetRadius(0.01);
	voxeledgetube->SetNumberOfSides(10);
	voxeledgetube->SetDefaultNormal(0.577,0.577,0.577);
	
vtkPolyDataMapper *voxeledgemapper=vtkPolyDataMapper::New();
	voxeledgemapper->SetInputConnection(voxeledgetube->GetOutputPort());
	
vtkActor *voxeledgeactor=vtkActor::New();                            //actor for voxel edge
	voxeledgeactor->SetMapper(voxeledgemapper);
	voxeledgeactor->GetProperty()->SetColor(1,1,1);
	
vtkDataSetMapper *boundingvolumemapper=vtkDataSetMapper::New();
	boundingvolumemapper->SetInputData(boundingvolume);
	
vtkActor *boundingvolumeactor=vtkActor::New();                        // actor for bounding volume
	boundingvolumeactor->SetMapper(boundingvolumemapper);
	boundingvolumeactor->GetProperty()->SetRepresentationToWireframe();
	boundingvolumeactor->GetProperty()->SetOpacity(0.5);
	
/***************** Section seven  voxelization process ends******************************/	
	

	
/***************** Section eight  create 3D Glyph for voxel starts******************************/

	
vtkSphereSource *Sphere=vtkSphereSource::New();
	Sphere->SetRadius((((*x_max)-(*x_min))/x_width+1)/50);
	Sphere->SetPhiResolution(20);
	Sphere->SetThetaResolution(20);

vtkThresholdPoints *ThresholdIn2 =vtkThresholdPoints::New();
	ThresholdIn2->SetInputData(boundingvolume);
	ThresholdIn2->ThresholdByUpper(0.0);

vtkGlyph3D *Vertices =vtkGlyph3D::New();
	Vertices->SetInputConnection(ThresholdIn2->GetOutputPort());
	Vertices->SetSourceData(Sphere->GetOutput());
	Vertices->SetScaleModeToScaleByScalar();
    Vertices->SetScaleFactor(2.0);

vtkPolyDataMapper *SphereMapper=vtkPolyDataMapper::New();
	SphereMapper->SetInputConnection(Vertices->GetOutputPort());
	SphereMapper->ScalarVisibilityOff();

vtkActor *sphereActor=vtkActor::New();
	sphereActor->SetMapper(SphereMapper);
	sphereActor->GetProperty()->SetDiffuseColor(1,0,0);

/***************** Section eight  create 3D Glyph for voxel ends******************************/


/***************** Section nine create the origin point for coordinate system******************************/

vtkAnnotatedCubeActor *cube=vtkAnnotatedCubeActor::New();
	cube->SetScale(3,3,3);
	cube->SetCubeVisibility(1);

vtkAxesActor *axes1=vtkAxesActor::New();
	axes1->SetShaftTypeToCylinder();
	axes1->SetXAxisLabelText("x");
	axes1->SetYAxisLabelText("Y");
	axes1->SetZAxisLabelText("Z");
	axes1->SetScale(1,1,1);
	axes1->SetTotalLength(3,3,3);

vtkPropAssembly *assembly =vtkPropAssembly::New();
	assembly->AddPart(axes1);
	assembly->AddPart(cube);

vtkOrientationMarkerWidget *marker=vtkOrientationMarkerWidget::New();
	marker->SetOutlineColor(0,0,0);
	marker->SetOrientationMarker(assembly);
	marker->SetViewport(0.0,0.0,0.15,0.3);	
	
/***************** Section nine create the origin point for coordinate system******************************/



/***************** Section ten  create the vtk environment stuff******************************/

vtkRenderer *ren = vtkRenderer::New();
	ren->SetBackground(0,0,0);
	
	ren->AddViewProp(scalarbar);
	//ren->AddViewProp(sphereActor);
	//ren->AddViewProp(boundingvolumeactor);
	//ren->AddViewProp(voxeledgeactor);
	ren->AddViewProp(textactor);
	//ren->AddViewProp(planeactor);
	ren->AddViewProp(pointsactor);
	//ren->AddViewProp(outlineactor);
	
	
	
	//ren->AddViewProp(cube);
	
	//ren->AddViewProp(assembly);


/////////////////////////////////////////////////////////////////////

vtkCubeAxesActor2D *axes=vtkCubeAxesActor2D::New();
    axes->SetInputData(normals->GetOutput());
	axes->SetCamera(ren->GetActiveCamera());
	axes->SetLabelFormat("%6.4g");
	axes->SetFontFactor(0.8);
	axes->SetFlyModeToOuterEdges();
	axes->SetAxisTitleTextProperty (tprop);
	axes->SetAxisLabelTextProperty(tprop);
	axes->ScalingOn();

vtkCamera *cam1 = vtkCamera::New();
	cam1->SetPosition((*x_min)+((*x_max)-(*x_min))/2,(*y_min)+((*y_max)-(*y_min))/2,(*z_min)+2*((*z_max)-(*z_min)));
	cam1->SetFocalPoint(normals->GetOutput()->GetCenter());
	//cam1->Azimuth(45);
	//cam1->Elevation(0);
	cam1->OrthogonalizeViewUp();
	

	cam1->SetViewUp(0, 1,0 );


	ren->AddViewProp(axes);
	ren->SetActiveCamera(cam1);


/////////////////////////////////////////////////////////////////////

	
vtkRenderWindow* renWin = vtkRenderWindow::New();
	renWin->AddRenderer(ren);
	renWin->SetSize(1024,768);
	renWin->SetStereoTypeToAnaglyph();
	renWin->Render();
	renWin->SetWindowName("RSGAL_Guang Zheng_2010");
	renWin->Modified();
	
vtkRenderWindowInteractor* iren = vtkRenderWindowInteractor::New();
	iren->SetRenderWindow(renWin);

vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New();
	iren->SetInteractorStyle(style);
	iren->Initialize();
	iren->Start();

/***************** Section ten  create the vtk environment stuff******************************/
	
	
	

pointsactor->Delete();
scalarbar->Delete();
outlineactor->Delete();
axes->Delete();
planeactor->Delete();
textactor->Delete();
voxeledgeactor->Delete();
boundingvolumeactor->Delete();
sphereActor->Delete();
cube->Delete();
assembly->Delete();


return 0;



}