//=========================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt 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 "vtkEnclosingModelEntityOperator.h"

#include "vtkBoundingBox.h"
#include "vtkCellData.h"
#include "vtkCellLocator.h"
#include "vtkDiscreteModel.h"
#include "vtkDiscreteModelFace.h"
#include "vtkDiscreteModelWrapper.h"
#include "vtkModelUserName.h"
#include "vtkDoubleArray.h"
#include "vtkGenericCell.h"
#include "vtkIdTypeArray.h"
#include "vtkMath.h"

#include "vtkModelRegion.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"
#include "vtkPolygon.h"
#include "vtkNew.h"


vtkStandardNewMacro(vtkEnclosingModelEntityOperator);

//----------------------------------------------------------------------------
vtkEnclosingModelEntityOperator::vtkEnclosingModelEntityOperator()
{
  this->OperateSucceeded = 0;
  this->CellLocator = 0;
  this->EnclosingEntity = 0;
}

//----------------------------------------------------------------------------
vtkEnclosingModelEntityOperator::~vtkEnclosingModelEntityOperator()
{
  if (this->CellLocator)
    {
    this->CellLocator->Delete();
    }
}

//----------------------------------------------------------------------------
void vtkEnclosingModelEntityOperator::BuildLinks(vtkDiscreteModelWrapper* modelWrapper)
{
  modelWrapper->GetModel()->GetMesh().BuildLinks();
}

//----------------------------------------------------------------------------
void vtkEnclosingModelEntityOperator::Operate(vtkDiscreteModelWrapper* modelWrapper,
                                              double *pt)
{
  vtkDebugMacro("Operating on a model.");
  this->OperateSucceeded = 0;

  const DiscreteMesh& mesh = modelWrapper->GetModel()->GetMesh();

  if (!this->CellLocator)
    {
    vtkSmartPointer<vtkPolyData> mainPD = mesh.ShallowCopyFaceData();
    this->CellLocator = vtkCellLocator::New();
    this->CellLocator->SetDataSet( mainPD );
    this->CellLocator->BuildLocator();
    }

  // reset, per chance we don't handle condition specific to the point and closest cell
  this->EnclosingEntity = 0;

  vtkIdType closestCellId;
  int subId;
  double dist2, closestPt[3];
  vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
  this->CellLocator->FindClosestPoint(pt, closestPt, cell, closestCellId, subId, dist2);

  // if closest isn't a poly, then start asking for cells within a bounds, and check polys
  // that are returned for which is closest?
//    virtual void FindCellsWithinBounds(double *bbox, vtkIdList *cells);

  if (cell->GetCellType() != VTK_TRIANGLE && cell->GetCellType() != VTK_POLYGON)
    {
    vtkErrorMacro("Check point location not supported if closest cell isn't a polygon!");
    return;
    }

  // determine whether we are on the positive side or negative side of the
  // polygon... note, need to handle cases where closest to edge of the polygon

  double normal[3];
  mesh.ComputeCellNormal(closestCellId,normal);


  vtkDiscreteModel::ClassificationType& classified =
      modelWrapper->GetModel()->GetMeshClassification();

  vtkModelEntity* entity = classified.GetEntity(closestCellId)->GetThisModelEntity();

  vtkDiscreteModelFace *faceEntity = vtkDiscreteModelFace::SafeDownCast(entity);

  // don't want failure due to bounds test
  vtkBoundingBox bbox(cell->GetBounds());
  bbox.AddPoint(pt);
  double bounds[6];
  bbox.GetBounds(bounds);

  double *cellPts =
    vtkDoubleArray::SafeDownCast(cell->GetPoints()->GetData())->GetPointer(0);
  int inside = vtkPolygon::PointInPolygon(pt, cell->GetNumberOfPoints(),
    cellPts, bounds, normal);

  // if "inside", can do simple normal test... more interesting if "outside", which
  // will happen plenty often
  if (inside != 1 )
    {
    // check the edges... which one are we closest to;  then check polys
    // that share that edge;  if dot of vector from centroid of poly1 ("cell")
    // to centroid of poly2 ("neighborCell") with normal vector of poly1 > 0,
    // then the point is "inside" region 0 (of face); else if < 0, then point
    // is outside region 0 (inside region 1) of the face
    int closestEdge = -1;
    double pCoords[3], weights[2], bestDist2 = -1,  edgeDist2;
    for (int i = 0; i < cell->GetNumberOfEdges(); i++)
      {
      vtkCell *edge = cell->GetEdge(i);
      if (edge->EvaluatePosition(pt, closestPt, subId, pCoords, edgeDist2, weights) == -1)
        {
        vtkErrorMacro("numerical error while anlysinz edge");
        continue;
        }
      if (closestEdge == -1 || edgeDist2 < bestDist2)
        {
        closestEdge = i;
        bestDist2 = edgeDist2;
        }
      }

    vtkIdType ptIds[2] = {cell->GetPointId(closestEdge),
      cell->GetPointId((closestEdge+1) % cell->GetNumberOfPoints())};

    vtkSmartPointer<vtkIdList> neighborIds  = vtkSmartPointer<vtkIdList>::New();
    mesh.GetCellEdgeNeighbors(closestCellId, ptIds[0], ptIds[1], neighborIds);

    if (neighborIds->GetNumberOfIds() == 1)
      {
      // centroid of neighber cell
      double neighborCentroid[3];
      mesh.ComputeCellCentroid(neighborIds->GetId(0), neighborCentroid);

      double centroid[3];
      mesh.ComputeCellCentroid(closestCellId, centroid);

      // vector from centroid of cell to centroid of neighborCell
      double centroidVector[3] = {neighborCentroid[0] - centroid[0],
        neighborCentroid[1] - centroid[1], neighborCentroid[2] - centroid[2]};

      if (vtkMath::Dot(centroidVector, normal) > 0)
        {
        this->EnclosingEntity = faceEntity->GetModelRegion(0);
        }
      else
        {
        this->EnclosingEntity = faceEntity->GetModelRegion(1);
        }
      this->OperateSucceeded = 1;
      return;
      }
    else
      {
      vtkErrorMacro("Condition not handled!");
      return;
      }
    }

  // if inside, check to see if "on" the polygon

  // check to see if outside or inside the region specified
  double testVector[3] = {pt[0] - cellPts[0], pt[1] - cellPts[1],
    pt[2] - cellPts[2]};

  if (vtkMath::Dot(testVector, normal) > 0)
    {
    this->EnclosingEntity = faceEntity->GetModelRegion(1);
    }
  else
    {
    this->EnclosingEntity = faceEntity->GetModelRegion(0);
    }

  vtkDebugMacro("Finished operating on a model.");
  this->OperateSucceeded = 1;
  return;
}

//----------------------------------------------------------------------------
void vtkEnclosingModelEntityOperator::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  os << indent << "OperateSucceeded: " << this->OperateSucceeded << endl;
}
