#include "vtkSphereSource.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkCellLocator.h"
#include "vtkPolyData.h"
#include "vtkPointData.h"
#include "vtkGenericCell.h"

int main( int argc, char *argv[] )
{

  vtkDataArray *normalsSphere1;
  double actualPoint[3], normalVectorSphere1[3], endPoint[3];
  double t, intersect[3], pcoords[3];
  int subId, cellId;
  vtkGenericCell *cell = vtkGenericCell::New();

  int numberOfFoundIntersections = 0;
  int numberOfNotFoundIntersections = 0;

  vtkSphereSource *sphere1 = vtkSphereSource::New();
		sphere1->SetThetaResolution(100);
    sphere1->SetPhiResolution(100);
		sphere1->SetRadius(1);
    sphere1->Update();

  vtkSphereSource *sphere2 = vtkSphereSource::New();
		sphere2->SetThetaResolution(100);
    sphere2->SetPhiResolution(100);
		sphere2->SetRadius(0.8);
    sphere2->Update();

  // Normals of sphere2
	normalsSphere1 = sphere1->GetOutput()->GetPointData()->GetNormals();

  vtkCellLocator *locator = vtkCellLocator::New();
    locator->SetDataSet((vtkDataSet*)sphere2->GetOutput());
    locator->CacheCellBoundsOn();
    locator->AutomaticOn();
		locator->BuildLocator();

  // this for loop changes the search depth
  for (double searchDepth = 0.05; searchDepth < 0.75; searchDepth = searchDepth + 0.05){
    //reset counting variables
    numberOfFoundIntersections = 0;
    numberOfNotFoundIntersections = 0;

    // this loop traverses all points on sphere1 and looks for an intersection on sphere2
    for (int i=0; i < sphere1->GetOutput()->GetNumberOfPoints(); i++){
      sphere1->GetOutput()->GetPoint(i, actualPoint);
		  normalsSphere1->GetTuple(i, normalVectorSphere1);

      // search in negative direction for an intersection with sphere1
		  endPoint[0] = actualPoint[0] - searchDepth * normalVectorSphere1[0];
		  endPoint[1] = actualPoint[1] - searchDepth * normalVectorSphere1[1];
		  endPoint[2] = actualPoint[2] - searchDepth * normalVectorSphere1[2];

      if(locator->IntersectWithLine(actualPoint, endPoint, 0.001, t, intersect, pcoords, subId, cellId, cell)){
        numberOfFoundIntersections++;
      }
      else{
        numberOfNotFoundIntersections++;
      }
    }

    cout << "actual search depth                              : " << searchDepth << endl;
    cout << "number of points where an intersection was found : " << numberOfFoundIntersections << endl;
    cout << "number of points with no intersection            : " << numberOfNotFoundIntersections << endl;
    cout << endl;
  }

  cout << endl;
  cout << "difference of the two radii is 0.2" << endl;

  return 0;
}