// This example creates an unstructured grid and attempts to choose
// one of the cells of the grid with a vtkCellPicker.  The chosen
// point is inside the grid, but it ought to return a point on the
// surface.

// With the given geometry, when vtkCellLocator::IntersectWithLine is
// examining the tetrahedron with cellID=25,
// vtkTetra::IntersectWithLine returns the wrong location.  This
// happens because when it's examining the triangle with corners at
// (0, 5, 10), (5, 5, 10), and (5, 10, 10), it computes an
// intersection point at z = 10 + 1.78e-15, and concludes that it's
// outside of the locator's octant, which is limited by z<=10. 
// vtkTetra::IntersectWithLine therefore returns an intersection point
// on a different face, which is farther away from the camera, and not
// on the exterior surface of the grid at z=10.

// Changing the camera or grid parameters very slightly from the
// values used here can make the effect go away, although there are
// other values for which the effect remains.

// This was tested with vtk 5.10.1 on the following 64-bit systems:
//   OS X 10.6.8, g++ 4.2.2
//   OS X 10.7.5, clang++ 4.2
//   Debian Linux 6.0, g++ 4.4 
// No special compilation flags.  Link with -lvtkCommon -lvtkFiltering
// -lvtkRendering -lvtkGraphics -lvtkHybrid.


#include <iostream>
#include <math.h>

#include <vtkActor.h>
#include <vtkAxesActor.h>
#include <vtkCamera.h>
#include <vtkCaptionActor2D.h>
#include <vtkCellLocator.h>
#include <vtkCellPicker.h>
#include <vtkDataSetMapper.h>
#include <vtkPoints.h>
#include <vtkProp3DCollection.h>
#include <vtkProperty.h>
#include <vtkProperty2D.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkSmartPointer.h>
#include <vtkTetra.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>
#include <vtkTransform.h>
#include <vtkUnstructuredGrid.h>
#include <vtkXOpenGLRenderWindow.h>

int main(int, char**) {

  // Divide a cube into subcubes, each split into 5 tetrahedra.
  double gridsize = 10.;	// length of cube side
  int n = 2;			// number of subcubes per side
  int dx = gridsize/n;

  double center[] = {gridsize/2, gridsize/2, gridsize/2};

  // Make the points
  vtkSmartPointer<vtkPoints> points =
    vtkSmartPointer<vtkPoints>::New();
  int nTotal = (n+1)*(n+1)*(n+1);
  points->Allocate(nTotal, nTotal);
  for(int i=0; i<=n; i++) {     // y
    for(int j=0; j<=n; j++) {	// x
      for(int k=0; k<=n; k++) { // z
	double x[3];
	if(j == n)
	  x[0] = gridsize;
	else
	  x[0] = j*dx;
	if(i == n)
	  x[1] = gridsize;
	else
	  x[1] = i*dx;
	if(k==n)
	  x[2] = gridsize;
	else
	  x[2] = k*dx;
	points->InsertNextPoint(x);
      }
    }
  }

  vtkSmartPointer<vtkUnstructuredGrid> grid =
    vtkSmartPointer<vtkUnstructuredGrid>::New();
  grid->SetPoints(points);

  // Make the tetrahedra
  bool flip = false;
  grid->Allocate(5*n*n*n, 5*n*n*n);
  for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) {
      for(int k=0; k<n; k++) {
	int ulf = (i+1)*(n+1)*(n+1)+j*(n+1)+k+1;      // upper left front
	int urf = (i+1)*(n+1)*(n+1)+(j+1)*(n+1)+k+1;  // upper right front
	int lrf = i*(n+1)*(n+1)+(j+1)*(n+1)+k+1;      // lower right front
	int llf = i*(n+1)*(n+1)+j*(n+1)+k+1;          // lower left front 
	int ulb = (i+1)*(n+1)*(n+1)+j*(n+1)+k;        // upper left back
	int urb = (i+1)*(n+1)*(n+1)+(j+1)*(n+1)+k;    // upper right back 
	int lrb = i*(n+1)*(n+1)+(j+1)*(n+1)+k;        // lower right back
	int llb = i*(n+1)*(n+1)+j*(n+1)+k;            // lower left back

	if(!flip) {
	  vtkIdType ids1[4] = {llf, urf, lrf, lrb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids1);
	  vtkIdType ids2[4] = {llf, ulf, urf, ulb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids2);
	  vtkIdType ids3[4] = {lrb, urf, urb, ulb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids3);
	  vtkIdType ids4[4] = {llf, lrb, llb, ulb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids4);
	  vtkIdType ids5[4] = {llf, ulb, urf, lrb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids5);
	}
	else {
	  vtkIdType ids1[4] = {llf, ulf, lrf, llb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids1);
	  vtkIdType ids2[4] = {ulf, urf, lrf, urb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids2);
	  vtkIdType ids3[4] = {ulf, ulb, urb, llb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids3);
	  vtkIdType ids4[4] = {lrf, urb, lrb, llb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids4);
	  vtkIdType ids5[4] = {ulf, urb, lrf, llb};
	  grid->InsertNextCell(VTK_TETRA, 4, ids5);
	}
	flip = !flip;
      }
      flip = !flip;
    }
    flip = !flip;
  }

  vtkSmartPointer<vtkXOpenGLRenderWindow> render_window =
    vtkSmartPointer<vtkXOpenGLRenderWindow>::New();
  vtkSmartPointer<vtkRenderer> renderer =
    vtkSmartPointer<vtkRenderer>::New();
  render_window->AddRenderer(renderer);

  // Draw the grid.
  vtkSmartPointer<vtkActor> faceActor = 
    vtkSmartPointer<vtkActor>::New();
  vtkSmartPointer<vtkDataSetMapper> faceMapper = 
    vtkSmartPointer<vtkDataSetMapper>::New();
  faceMapper->SetInputConnection(grid->GetProducerPort());
  faceActor->SetMapper(faceMapper);
  faceActor->GetProperty()->SetRepresentationToSurface();
  faceActor->PickableOn();
  renderer->AddViewProp(faceActor);

  // Draw axes.
  vtkSmartPointer<vtkAxesActor> axes = vtkSmartPointer<vtkAxesActor>::New();
  renderer->AddActor(axes);
  axes->SetTotalLength(11., 11., 11.);
  vtkSmartPointer<vtkTransform> transf = vtkSmartPointer<vtkTransform>::New();
  transf->Translate(-0.5, -0.5, -0.5);
  axes->SetUserTransform(transf);

  vtkSmartPointer<vtkCellLocator> locator =
    vtkSmartPointer<vtkCellLocator>::New();
  locator->Initialize();
  locator->LazyEvaluationOn();
  grid->Update();
  locator->SetDataSet(grid);

  // Compute initial view parameters

  // depth is the distance from the center of the grid to its edge, in
  // the direction of the camera.
  double depth = center[2];
  // width is the same thing, but perpendicular to the camera direction.
  double width = center[0];
  double angle = 30.;
  // dist is the distance from the front of the grid to the camera, so
  // that the whole grid is visible, including a fudge factor.
  double dist = 1.3*width/tan(M_PI/180.*angle/2.);

  double cameraZ = center[2] + depth + dist;
  
  // The exact value of the camera position here is crucial.  Using
  // "34.2583302491977" instead of cameraZ moves the
  // intersection point across the octant boundary (by a distance of
  // 2*1.77635683940025e-15)

  double cameraPos0[] = {center[0], center[1], cameraZ};
  double cameraUp0[] = {0., 1., 0.};
  vtkSmartPointer<vtkCamera> camera = renderer->GetActiveCamera();
  camera->SetPosition(cameraPos0);
  camera->SetFocalPoint(center);
  //camera->SetViewUp(cameraUp0);
  renderer->ResetCameraClippingRange();

  // Magic parameters here come from the script driving the original
  // program, and were generated there by recording a GUI session.
  render_window->SetSize(691, 652);
  double cameraPos1[] = {-7.32641, 7.86392, 31.3801};
  double cameraUp1[] = {0.000351236, 0.994176, -0.107767};
  camera->SetPosition(cameraPos1);
  camera->SetFocalPoint(center);
  camera->SetViewUp(cameraUp1);

  vtkSmartPointer<vtkCellPicker> picker =
    vtkSmartPointer<vtkCellPicker>::New();
  picker->AddLocator(locator);

  // x and y are the position of a mouse click.  Changing y to 353
  // also puts the picked point inside the cube.  352 puts it on the
  // surface.
  double x = 448.;
  double y = 351.;
  
  if(picker->Pick(x, y, 0.0, renderer)) {
    vtkSmartPointer<vtkPoints> pickedPoints = picker->GetPickedPositions();
    double pickedPoint[3];
    pickedPoints->GetPoint(0, pickedPoint);
    std::cout 
      << "The picked point should be on the surface of the cube, at z=10." 
      << std::endl;
    std::cout << "Picked point:"
	      << " x=" << pickedPoint[0] 
	      << ", y = " << pickedPoint[1]
	      << ", z = " << pickedPoint[2]
	      << std::endl;
  }
  else 
    std::cout << "Nothing picked" << std::endl;

  return 0;
}