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

  Program:   Visualization Toolkit
  Module:    vtkPolyDataSilhouette.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.

=========================================================================*/
// .SECTION Thanks
// Contribution by Thierry Carrard <br>
// CEA/DIF - Commissariat a l'Energie Atomique, Centre DAM Ile-De-France <br>
// BP12, F-91297 Arpajon, France. <br>

#include "vtkPolyDataSilhouette.h"

#include "vtkCamera.h"
#include "vtkCellData.h"
#include "vtkGenericCell.h"
#include "vtkMath.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkProp3D.h"
#include "vtkTransform.h"
#include "vtkUnsignedIntArray.h"
#include "vtkIdTypeArray.h"
#include "vtkPoints.h"
#include "vtkCellArray.h"
#include "vtkPolygon.h"

#include <map>

vtkStandardNewMacro(vtkPolyDataSilhouette);

vtkCxxSetObjectMacro(vtkPolyDataSilhouette,Camera,vtkCamera);

struct vtkOrderedEdge
{
  inline vtkOrderedEdge(vtkIdType a, vtkIdType b)
    {
    if(a<=b) { p1=a; p2=b; }
    else     { p1=b; p2=a; }
    }
  inline bool operator < (const vtkOrderedEdge& oe) const
    {
    return (p1<oe.p1) || ( (p1==oe.p1) && (p2<oe.p2) ) ;
    }
  vtkIdType p1,p2;
};

struct vtkTwoNormals
{
  double leftNormal[3]; // normal of the left polygon
  double rightNormal[3]; // normal of the right polygon
  inline vtkTwoNormals()
    {
    leftNormal[0] = 0.0;
    leftNormal[1] = 0.0;
    leftNormal[2] = 0.0;
    rightNormal[0] = 0.0;
    rightNormal[1] = 0.0;
    rightNormal[2] = 0.0;
    }
};

class vtkPolyDataEdges
{
public:
  vtkTimeStamp mtime;
  double vec[3];
  std::map<vtkOrderedEdge,vtkTwoNormals> edges;
  bool * edgeFlag;
  vtkCellArray* lines;
  inline vtkPolyDataEdges() : edgeFlag(0), lines(0) { vec[0]=vec[1]=vec[2]=0.0; }
};

vtkPolyDataSilhouette::vtkPolyDataSilhouette()
{
  this->Camera = NULL;
  this->Prop3D = NULL;
  this->Direction = VTK_DIRECTION_CAMERA_ORIGIN;
  this->Vector[0] = this->Vector[1] = this->Vector[2] = 0.0;
  this->Origin[0] = this->Origin[1] = this->Origin[2] = 0.0;
  this->Transform = vtkTransform::New();
  this->EnableFeatureAngle = 1;
  this->FeatureAngle = 60;
  this->BorderEdges = 0;
  this->PieceInvariant = 1;
  this->PreComp = new vtkPolyDataEdges();
}

vtkPolyDataSilhouette::~vtkPolyDataSilhouette()
{
  this->Transform->Delete();
  
  if ( this->Camera )
    {
    this->Camera->Delete();
    }

  if (this->PreComp->edgeFlag)
    {
    delete [] this->PreComp->edgeFlag;
    }
  if (this->PreComp->lines)
    {
    this->PreComp->lines->Delete();
    }
  delete this->PreComp;

  //Note: vtkProp3D is not deleted to avoid reference count cycle
}

// Don't reference count to avoid nasty cycle
void vtkPolyDataSilhouette::SetProp3D(vtkProp3D *prop3d)
{
  if ( this->Prop3D != prop3d )
    {
    this->Prop3D = prop3d;
    this->Modified();
    }
}

vtkProp3D *vtkPolyDataSilhouette::GetProp3D()
{
  return this->Prop3D;
}

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

  // get the input and ouptut
  vtkPolyData *input = vtkPolyData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  if(input==0 || output==0)
    {
    vtkErrorMacro(<<"Need correct connections");
    return 0;
    }

  vtkDebugMacro(<<"RequestData\n");

  const double featureAngleCos = cos( vtkMath::RadiansFromDegrees( this->FeatureAngle ) );

  bool vectorMode = true;
  double vector[3];
  double origin[3];

  // Compute the sort vector
  switch( this->Direction )
    {
  case VTK_DIRECTION_SPECIFIED_VECTOR :
    vector[0] = this->Vector[0];
    vector[1] = this->Vector[1];
    vector[2] = this->Vector[2];
    break;

  case VTK_DIRECTION_SPECIFIED_ORIGIN :
    origin[0] = this->Origin[0];
    origin[1] = this->Origin[1];
    origin[2] = this->Origin[2];
    vectorMode = false;
    break;

  case VTK_DIRECTION_CAMERA_ORIGIN :
    vectorMode = false;

  case VTK_DIRECTION_CAMERA_VECTOR :
    if ( this->Camera == NULL)
      {
      vtkErrorMacro(<<"Need a camera when direction is set to VTK_DIRECTION_CAMERA_*");
      return 0;
      }
    this->ComputeProjectionVector(vector, origin);
    break;
    }

  vtkIdType nPolys = input->GetNumberOfPolys();
  vtkIdTypeArray* polysArray = input->GetPolys()->GetData();
  vtkPoints* inPoints = input->GetPoints();

  if( input->GetMTime() > this->PreComp->mtime.GetMTime() )
    {
    vtkDebugMacro(<<"Compute edge-face connectivity and face normals\n");

    this->PreComp->mtime.Modified();
    this->PreComp->edges.clear();

    vtkIdType* polys = polysArray->GetPointer(0);

    for(vtkIdType i=0;i<nPolys;i++)
      {
      int np = *(polys); ++polys;
      double normal[3];
      vtkPolygon::ComputeNormal(inPoints, np, polys, normal);

      for(int j=0;j<np;j++)
        {
        vtkIdType p1=j, p2=(j+1)%np;
        vtkOrderedEdge oe( polys[p1], polys[p2] );
        vtkTwoNormals& tn = this->PreComp->edges[oe];
        if( polys[p1] < polys[p2] )
          {
#ifdef DEBUG
          if( vtkMath::Dot(tn.leftNormal,tn.leftNormal) > 0 )
            {
            vtkDebugMacro(<<"Warning: vtkPolyDataSilhouette: non-manifold mesh: edge-L ("<<polys[p1]<<","<<polys[p2]<<") counted more than once\n");
            }
#endif
          tn.leftNormal[0] = normal[0];
          tn.leftNormal[1] = normal[1];
          tn.leftNormal[2] = normal[2];
          }
        else
          {
#ifdef DEBUG
          if( vtkMath::Dot(tn.rightNormal,tn.rightNormal) > 0 )
            {
            vtkDebugMacro(<<"Warning: vtkPolyDataSilhouette: non-manifold mesh: edge-R ("<<polys[p1]<<","<<polys[p2]<<") counted more than once\n");
            }
#endif
          tn.rightNormal[0] = normal[0];
          tn.rightNormal[1] = normal[1];
          tn.rightNormal[2] = normal[2];
          }
        }
      polys += np;
      }

    if( this->PreComp->edgeFlag != 0 ) delete [] this->PreComp->edgeFlag;
    this->PreComp->edgeFlag = new bool[ this->PreComp->edges.size() ];
    }

  bool vecChanged = false;
  for(int d=0;d<3;d++)
    {
    vecChanged = vecChanged || this->PreComp->vec[d]!=vector[d] ;
    }

  if( ( this->PreComp->mtime.GetMTime() > output->GetMTime() ) ||
      ( this->Camera && this->Camera->GetMTime() > output->GetMTime() ) ||
      ( this->Prop3D && this->Prop3D->GetMTime() > output->GetMTime() ) ||
      vecChanged )
    {
    vtkDebugMacro(<<"Extract edges\n");

    vtkIdType i=0, silhouetteEdges=0;

    for(std::map<vtkOrderedEdge,vtkTwoNormals>::iterator it=this->PreComp->edges.begin(); it!=this->PreComp->edges.end(); ++it)
      {
      double d1,d2;

      // does this edge have two co-faces ?
      bool winged =  vtkMath::Norm(it->second.leftNormal)>0.5 && vtkMath::Norm(it->second.rightNormal)>0.5 ;

      // cosine of feature angle, to be compared with scalar product of two co-faces normals
      double edgeAngleCos = vtkMath::Dot( it->second.leftNormal, it->second.rightNormal );

      if( vectorMode ) // uniform direction
        {
        d1 = vtkMath::Dot( vector, it->second.leftNormal );
        d2 = vtkMath::Dot( vector, it->second.rightNormal );
        }
      else // origin to edge's center direction
        {
        double p1[3];
        double p2[3];
        double vec[3];
        inPoints->GetPoint( it->first.p1 , p1 );
        inPoints->GetPoint( it->first.p2 , p2 );
        vec[0] = origin[0] - ( (p1[0]+p2[0])*0.5 );
        vec[1] = origin[1] - ( (p1[1]+p2[1])*0.5 );
        vec[2] = origin[2] - ( (p1[2]+p2[2])*0.5 );
        d1 = vtkMath::Dot( vec, it->second.leftNormal );
        d2 = vtkMath::Dot( vec, it->second.rightNormal );     
        }

      // shall we output this edge ?
      bool outputEdge =
        ( winged && (d1*d2)<0 )
        || ( this->EnableFeatureAngle && edgeAngleCos<featureAngleCos )
        || ( this->BorderEdges && !winged );

      if( outputEdge ) // add this edge
        {
        this->PreComp->edgeFlag[i] = true ;
        ++silhouetteEdges;
        }
      else // skip this edge
        {
        this->PreComp->edgeFlag[i] = false ;
        }
      ++i;
      }


    // build output data set (lines)
    vtkIdTypeArray* la = vtkIdTypeArray::New();
    la->SetNumberOfValues( 3*silhouetteEdges );
    vtkIdType* laPtr = la->WritePointer(0,3*silhouetteEdges);

    i=0;
    silhouetteEdges=0;
    for(std::map<vtkOrderedEdge,vtkTwoNormals>::iterator it=this->PreComp->edges.begin(); it!=this->PreComp->edges.end(); ++it)
      {
      if( this->PreComp->edgeFlag[i] )
        {
        laPtr[ silhouetteEdges*3+0 ] = 2 ;
        laPtr[ silhouetteEdges*3+1 ] = it->first.p1 ;
        laPtr[ silhouetteEdges*3+2 ] = it->first.p2 ;
        ++silhouetteEdges;
        }
      ++i;
      }

    if( this->PreComp->lines == 0 )
      {
      this->PreComp->lines = vtkCellArray::New();
      }
    this->PreComp->lines->SetCells( silhouetteEdges, la );
    la->Delete();
    }

  output->Initialize();
  output->SetPoints(inPoints);
  output->SetLines(this->PreComp->lines);

  return 1;
}

void vtkPolyDataSilhouette::ComputeProjectionVector(double vector[3], 
                                                   double origin[3])
{
  double *focalPoint = this->Camera->GetFocalPoint();
  double *position = this->Camera->GetPosition();
 
  // If a camera is present, use it
  if ( !this->Prop3D )
    {
    for(int i=0; i<3; i++)
      { 
      vector[i] = focalPoint[i] - position[i];
      origin[i] = position[i];
      }
    }

  else  //Otherwise, use Prop3D
    {
    double focalPt[4], pos[4];
    int i;

    this->Transform->SetMatrix(this->Prop3D->GetMatrix());
    this->Transform->Push();
    this->Transform->Inverse();

    for(i=0; i<4; i++)
      {
      focalPt[i] = focalPoint[i];
      pos[i] = position[i];
      }

    this->Transform->TransformPoint(focalPt,focalPt);
    this->Transform->TransformPoint(pos,pos);

    for (i=0; i<3; i++) 
      {
      vector[i] = focalPt[i] - pos[i];
      origin[i] = pos[i];
      }
    this->Transform->Pop();
  }
}

unsigned long int vtkPolyDataSilhouette::GetMTime()
{
  unsigned long mTime=this->Superclass::GetMTime();
 
  if ( this->Direction != VTK_DIRECTION_SPECIFIED_VECTOR )
    {
    unsigned long time;
    if ( this->Camera != NULL )
      {
      time = this->Camera->GetMTime();
      mTime = ( time > mTime ? time : mTime );
      }

    if ( this->Prop3D != NULL )
      {
      time = this->Prop3D->GetMTime();
      mTime = ( time > mTime ? time : mTime );
      }
    }

  return mTime;
}

void vtkPolyDataSilhouette::PrintSelf(ostream& os, vtkIndent indent)
{

  this->Superclass::PrintSelf(os,indent);

  if ( this->Camera )
    {
    os << indent << "Camera:\n";
    this->Camera->PrintSelf(os,indent.GetNextIndent());
    }
  else
    {
    os << indent << "Camera: (none)\n";
    }

  if ( this->Prop3D )
    {
    os << indent << "Prop3D:\n";
    this->Prop3D->PrintSelf(os,indent.GetNextIndent());
    }
  else
    {
    os << indent << "Prop3D: (none)\n";
    }

  os << indent << "Direction: ";
#define DIRECTION_CASE(name) case VTK_DIRECTION_##name :os << "VTK_DIRECTION_" << #name <<"\n"; break
  switch(this->Direction)
    {
    DIRECTION_CASE( SPECIFIED_ORIGIN );
    DIRECTION_CASE( SPECIFIED_VECTOR );
    DIRECTION_CASE( CAMERA_ORIGIN );
    DIRECTION_CASE( CAMERA_VECTOR );
    }
#undef DIRECTION_CASE

  if( this->Direction == VTK_DIRECTION_SPECIFIED_VECTOR ) 
    {
    os << "Specified Vector: (" << this->Vector[0] << ", " << this->Vector[1] << ", " << this->Vector[2] << ")\n";
    }
  if( this->Direction == VTK_DIRECTION_SPECIFIED_ORIGIN )
    {
    os << "Specified Origin: (" << this->Origin[0] << ", " << this->Origin[1] << ", " << this->Origin[2] << ")\n";
    }

  os << indent << "PieceInvariant: " << this->PieceInvariant << "\n";
  os << indent << "FeatureAngle: " << this->FeatureAngle << "\n";
  os << indent << "EnableFeatureAngle: " << this->EnableFeatureAngle << "\n";
  os << indent << "BorderEdges: " << this->BorderEdges << "\n";
}