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

  Program:   Visualization Toolkit
  Module:    vtkSelectEnclosedPoints.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkSelectEnclosedPoints.h"

#include "vtkDataSet.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkPolyData.h"
#include "vtkUnsignedCharArray.h"
#include "vtkExecutive.h"
#include "vtkFeatureEdges.h"
#include "vtkCellLocator.h"
#include "vtkGenericCell.h"
#include "vtkMath.h"
#include "vtkGarbageCollector.h"

vtkStandardNewMacro(vtkSelectEnclosedPoints);

//----------------------------------------------------------------------------
// Construct object.
vtkSelectEnclosedPoints::vtkSelectEnclosedPoints()
{
  this->SetNumberOfInputPorts(2);

  this->CheckSurface = 0;
  this->InsideOut = 0;
  this->Tolerance = 0.001;

  this->InsideOutsideArray = NULL;
  
  this->CellLocator = vtkCellLocator::New();
  this->CellIds = vtkIdList::New();
  this->Cell = vtkGenericCell::New();
}

//----------------------------------------------------------------------------
vtkSelectEnclosedPoints::~vtkSelectEnclosedPoints()
{
  if ( this->InsideOutsideArray )
    {
    this->InsideOutsideArray->Delete();
    }

  if ( this->CellLocator )
    {
    vtkCellLocator *loc = this->CellLocator;
    this->CellLocator = NULL;
    loc->Delete();
    }

  this->CellIds->Delete();
  this->Cell->Delete();
}

//----------------------------------------------------------------------------
int vtkSelectEnclosedPoints::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *in2Info = inputVector[1]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  // get the two inputs and output
  vtkDataSet *input = vtkDataSet::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *surface = vtkPolyData::SafeDownCast(
    in2Info->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataSet *output = vtkDataSet::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkDebugMacro("Selecting enclosed points");

  // If requested, check that the surface is closed
  if ( this->CheckSurface && ! this->IsSurfaceClosed(surface) )
    {
    return 0;
    }

  // Initiailize search structures
  this->Initialize(surface);

  // Create array to mark inside/outside
  if ( this->InsideOutsideArray )
    {
    this->InsideOutsideArray->Delete();
    }
  this->InsideOutsideArray = vtkUnsignedCharArray::New();
  vtkUnsignedCharArray *marks = this->InsideOutsideArray;
  marks->SetName("SelectedPointsArray");

  // Loop over all input points determining inside/outside
  vtkIdType numPts = input->GetNumberOfPoints();
  marks->SetNumberOfValues(numPts);
  vtkIdType ptId;
  double x[3];
  
  int abort=0;
  vtkIdType progressInterval=numPts/20+1;
  for ( ptId=0; ptId < numPts && !abort; ptId++ )
    {
    if ( ! (ptId % progressInterval) ) //manage progress / early abort
      {
      this->UpdateProgress ((double)ptId / numPts);
      abort = this->GetAbortExecute();
      }

    input->GetPoint(ptId,x);

    if ( this->IsInsideSurface(x) )
      {
      marks->SetValue(ptId,(this->InsideOut?0:1));
      }
    else
      {
      marks->SetValue(ptId,(this->InsideOut?1:0));
      }
    }
  
  // Copy all the input geometry and data to the output.
  output->CopyStructure(input);
  output->GetPointData()->PassData(input->GetPointData());
  output->GetCellData()->PassData(input->GetCellData());

  // Add the new scalars array to the output.
  marks->SetName("SelectedPoints");
  output->GetPointData()->SetScalars(marks);

  // release memory
  this->Complete();

  return 1;
}

//----------------------------------------------------------------------------
int vtkSelectEnclosedPoints::IsSurfaceClosed(vtkPolyData *surface)
{
  vtkPolyData *checker = vtkPolyData::New();
  checker->CopyStructure(surface);
  
  vtkFeatureEdges *features = vtkFeatureEdges::New();
  features->SetInput(checker);
  features->BoundaryEdgesOn();
  features->NonManifoldEdgesOn();
  features->ManifoldEdgesOff();
  features->FeatureEdgesOff();
  features->Update();
  
  vtkIdType numCells = features->GetOutput()->GetNumberOfCells();
  features->Delete();
  checker->Delete();
  
  if ( numCells > 0 )
    {
    return 0;
    }
  else
    {
    return 1;
    }
}

//----------------------------------------------------------------------------
void vtkSelectEnclosedPoints::Initialize(vtkPolyData *surface)
{
  if ( ! this->CellLocator )
    {
    this->CellLocator = vtkCellLocator::New();
    }

  this->Surface = surface;
  surface->GetBounds(this->Bounds);
  this->Length = surface->GetLength();

  // Set up structures for acceleration ray casting
  this->CellLocator->SetDataSet(surface);
  this->CellLocator->BuildLocator();
}

//----------------------------------------------------------------------------
int vtkSelectEnclosedPoints::IsInside(vtkIdType inputPtId)
{
  if ( !this->InsideOutsideArray ||
       this->InsideOutsideArray->GetValue(inputPtId) == 0 )
    {
    return 0;
    }
  else 
    {
    return 1;
    }
}

//----------------------------------------------------------------------------
int vtkSelectEnclosedPoints::IsInsideSurface(double x, double y, double z)
{
  double xyz[3];
  xyz[0] = x;
  xyz[1] = y;
  xyz[2] = z;
  return this->IsInsideSurface(xyz);
}


#define VTK_CERTAIN 1
#define VTK_UNCERTAIN 0
#define VTK_MAX_ITER 10    //Maximum iterations for ray-firing
#define VTK_VOTE_THRESHOLD 3
//----------------------------------------------------------------------------
int vtkSelectEnclosedPoints::IsInsideSurface(double x[3])
{
  // do a quick bounds check
  if ( x[0] < this->Bounds[0] || x[0] > this->Bounds[1] ||
       x[1] < this->Bounds[2] || x[1] > this->Bounds[3] ||
       x[2] < this->Bounds[4] || x[2] > this->Bounds[5])
    {
    return 0;
    }
  
  //  Perform in/out by shooting random rays. Multiple rays are fired
  //  to improve accuracy of the result.
  //
  //  The variable iterNumber counts the number of rays fired and is
  //  limited by the defined variable VTK_MAX_ITER.
  //
  //  The variable deltaVotes keeps track of the number of votes for
  //  "in" versus "out" of the surface.  When deltaVotes > 0, more votes
  //  have counted for "in" than "out".  When deltaVotes < 0, more votes
  //  have counted for "out" than "in".  When the delta_vote exceeds or
  //  equals the defined variable VTK_VOTE_THRESHOLD, then the
  //  appropriate "in" or "out" status is returned.
  //
  double rayMag, ray[3], xray[3], t, pcoords[3], xint[3];
  int i, numInts, iterNumber, deltaVotes, subId;
  vtkIdType idx, numCells;
  double tol = this->Tolerance*this->Length;

  for (deltaVotes = 0, iterNumber = 1;
       (iterNumber < VTK_MAX_ITER) && (abs(deltaVotes) < VTK_VOTE_THRESHOLD);
       iterNumber++) 
    {
    //  Define a random ray to fire.
    rayMag = 0.0;
    while (rayMag == 0.0 )
      {
      for (i=0; i<3; i++)
        {
        ray[i] = vtkMath::Random(-1.0,1.0);
        }
      rayMag = vtkMath::Norm(ray);
      }

    // The ray must be appropriately sized wrt the bounding box. (It has to go
    // all the way through the bounding box.)
    for (i=0; i<3; i++)
      {
      xray[i] = x[i] + (this->Length/rayMag)*ray[i];
      }

    // Retrieve the candidate cells from the locator
    this->CellLocator->FindCellsAlongLine(x,xray,tol,this->CellIds);

    // Intersect the line with each of the candidate cells
    numInts = 0;
    numCells = this->CellIds->GetNumberOfIds();
    for ( idx=0; idx < numCells; idx++ )
      {
      this->Surface->GetCell(this->CellIds->GetId(idx), this->Cell);
      if ( this->Cell->IntersectWithLine(x, xray, tol, t, xint, pcoords, subId) )
        {
        numInts++;
        }
      } //for all candidate cells
    
    // Count the result
    if ( (numInts % 2) == 0)
      {
      --deltaVotes;
      }
    else
      {
      ++deltaVotes;
      }
    } //try another ray

  //   If the number of votes is positive, the point is inside
  //
  return ( deltaVotes < 0 ? 0 : 1 );
}

#undef VTK_CERTAIN
#undef VTK_UNCERTAIN
#undef VTK_MAX_ITER
#undef VTK_VOTE_THRESHOLD


//----------------------------------------------------------------------------
// Specify a source object at a specified table location.
void vtkSelectEnclosedPoints::SetSurfaceConnection(vtkAlgorithmOutput* algOutput)
{
  this->SetInputConnection(1, algOutput);
}

//----------------------------------------------------------------------------
// Specify a source object at a specified table location.
void vtkSelectEnclosedPoints::SetSurface(vtkPolyData *pd)
{
  this->SetInput(1, pd);
}

//----------------------------------------------------------------------------
// Get a pointer to a source object at a specified table location.
vtkPolyData *vtkSelectEnclosedPoints::GetSurface()
{
  return vtkPolyData::SafeDownCast(
    this->GetExecutive()->GetInputData(1, 0));
}

//----------------------------------------------------------------------------
vtkPolyData* vtkSelectEnclosedPoints::GetSurface(vtkInformationVector *sourceInfo)
{
  vtkInformation *info = sourceInfo->GetInformationObject(1);
  if (!info)
    {
    return NULL;
    }
  return vtkPolyData::SafeDownCast(info->Get(vtkDataObject::DATA_OBJECT()));
}

//----------------------------------------------------------------------------
int vtkSelectEnclosedPoints::FillInputPortInformation(int port, vtkInformation *info)
{
  if (port == 0)
    {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
    }
  else if (port == 1)
    {
    info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 0);
    info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 0);
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
    }

  return 1;
}

//----------------------------------------------------------------------------
void vtkSelectEnclosedPoints::Complete()
{
  this->CellLocator->FreeSearchStructure();
}

//----------------------------------------------------------------------------
void vtkSelectEnclosedPoints::ReportReferences(vtkGarbageCollector* collector)
{
  this->Superclass::ReportReferences(collector);
  // These filters share our input and are therefore involved in a
  // reference loop.
  vtkGarbageCollectorReport(collector, this->CellLocator, "CellLocator");
}

//----------------------------------------------------------------------------
void vtkSelectEnclosedPoints::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  os << indent << "Check Surface: "
     << (this->CheckSurface ? "On\n" : "Off\n");

  os << indent << "Inside Out: "
     << (this->InsideOut ? "On\n" : "Off\n");
  
  os << indent << "Tolerance: " << this->Tolerance << "\n";
}