//----------------------------------------------------------------------------
//  EgalTech 2018-2018
//----------------------------------------------------------------------------
// File : IntersTriaTria.cpp       Data : 27.08.18         Versione : 1.9h3
// Contenuto : Implementazione della intersezione triangolo/triangolo.
//
//
//
// Modifiche : 27.08.18 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "ProjPlane.h"
#include "IntersLineTria.h"
#include "CurveComposite.h"
#include "SurfFlatRegion.h"
#include "Triangulate.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EGkIntersTriaTria.h"
#include "/EgtDev/Include/EGkIntersPlanePlane.h"
#include <array>

using namespace std ;

//----------------------------------------------------------------------------
int nTriaTriaIntersCases[6][6] = { { ITTTS_NO,  ITTTS_NO,           ITTTS_NO,             ITTTS_NO,         ITTTS_NO,             ITTTS_NO }, 
                                   { ITTTS_NO,  ITTTS_INT_INT_SEG,  ITTTS_INT_EDGE,       ITTTS_NO,         ITTTS_NO,             ITTTS_NO },
                                   { ITTTS_NO,  ITTTS_EDGE_INT,     ITTTS_EDGE_EDGE_SEG,  ITTTS_NO,         ITTTS_NO,             ITTTS_NO },
                                   { ITTTS_NO,  ITTTS_NO,           ITTTS_NO,             ITTTS_VERT_VERT,  ITTTS_VERT_EDGE,      ITTTS_VERT_INT },
                                   { ITTTS_NO,  ITTTS_NO,           ITTTS_NO,             ITTTS_EDGE_VERT,  ITTTS_EDGE_EDGE_PNT,  ITTTS_NO },
                                   { ITTTS_NO,  ITTTS_NO,           ITTTS_NO,             ITTTS_INT_VERT,   ITTTS_NO,             ITTTS_NO } } ;

//----------------------------------------------------------------------------
static int IntersCoplanarTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, TRIA3DVECTOR& vTria) ;

//----------------------------------------------------------------------------
int
IntersTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, Point3d& ptInt, Point3d& ptInt2, TRIA3DVECTOR& vTria)
{
  // piano del secondo triangolo
   Plane3d plTria2 ;
   plTria2.Set( trTria2.GetCentroid(), trTria2.GetN()) ;
  // calcolo le distanze dei vertici del primo triangolo dal piano del secondo
   array< double, 3> vDist1 ;
   for ( int i = 0 ; i < 3 ; ++i)
      vDist1[i] = DistPointPlane( trTria1.GetP( i), plTria2) ;
  // verifico posizione del primo triangolo rispetto al piano del secondo
   int nVertPos1 = 0 ; int nVertNeg1 = 0 ;
   for ( const auto& dDist : vDist1) {
      if ( dDist > EPS_SMALL)
         ++ nVertPos1 ;
      else if ( dDist < -EPS_SMALL)
         ++ nVertNeg1 ;
   }
  // se il triangolo giace tutto da una parte del piano, nessuna intersezione
   if ( nVertPos1 == 3  ||  nVertNeg1 == 3)
      return ITTT_NO ;

  // piano del primo triangolo
   Plane3d plTria1 ;
   plTria1.Set( trTria1.GetCentroid(), trTria1.GetN()) ;
  // calcolo le distanze dei vertici del secondo triangolo dal piano del primo
   array< double, 3> vDist2 ;
   for ( int i = 0 ; i < 3 ; ++i)
      vDist2[i] = DistPointPlane( trTria2.GetP( i), plTria1) ;
  // verifico posizione del secondo triangolo rispetto al piano del primo
   int nVertPos2 = 0 ; int nVertNeg2 = 0 ;
   for ( const auto& dDist : vDist2) {
      if ( dDist > EPS_SMALL)
         ++ nVertPos2 ;
      else if ( dDist < -EPS_SMALL)
         ++ nVertNeg2 ;
   }
  // se il triangolo giace tutto da una parte del piano, nessuna intersezione
   if ( nVertPos2 == 3  ||  nVertNeg2 == 3)
      return ITTT_NO ;

  // intersezione tra i piani dei due triangoli
   Point3d ptL1 ; Vector3d vtL1 ;
   int nResPP = IntersPlanePlane( plTria1, plTria2, ptL1, vtL1) ;
   if ( nResPP == IPPT_NO)
      return ITTT_NO ;

  // se i triangoli sono complanari
   if ( nResPP == IPPT_OVERLAPS  ||
        ((( nVertPos1 == 0  &&  nVertNeg1 == 0)  ||  ( nVertPos2 == 0  &&  nVertNeg2 == 0)) &&
         ( trTria1.GetN() ^ trTria2.GetN()).SqLen() < 25 * SIN_EPS_ANG_SMALL * SIN_EPS_ANG_SMALL))
      return IntersCoplanarTriaTria( trTria1, trTria2, vTria) ;

  // limito la linea di intersezione con il primo triangolo
   int nRes1 = IntersCoplanarLineTria( ptL1, vtL1, 100.0, trTria1, ptInt, ptInt2, false) ;
   if ( nRes1 == ILTT_NO)
      return ITTT_NO ;

  // il segmento calcolato va limitato col secondo triangolo
   Point3d ptL2 = ptInt ;
   Vector3d vtL2 = vtL1 ;
   double dLen = 0 ;
   if ( nRes1 != ILTT_VERT) {
      vtL2 =  ptInt2 - ptInt ;
      dLen = vtL2.Len() ;
      vtL2 /= dLen ;
   }  
   Point3d ptI1 = ptInt ;
   Point3d ptI12 = ptInt2 ;
   int nRes2 = IntersCoplanarLineTria( ptL2, vtL2, dLen, trTria2, ptInt, ptInt2, true) ;
  
  // gestione casi speciali apparentemente incongruenti
   if ( ( nRes1 == ILTT_SEGM  &&  nRes2 == ILTT_EDGE)  ||  ( nRes1 == ILTT_EDGE  &&  nRes2 == ILTT_SEGM)) {
      if ( nRes1 == ILTT_SEGM) {
         if ( PointInTria( ptI1, trTria2) != PTT_OUT  &&  PointInTria( ptI12, trTria2) != PTT_OUT) {
            ptInt = ( ptI1 + ptI12) / 2 ;
            ptInt2 = ( ptI1 + ptI12) / 2 ;
            return ITTTS_EDGE_EDGE_PNT ;
         }
         else
            return ITTTS_NO ;
      }
      else {
         Point3d ptI3, ptI32 ;
         IntersCoplanarLineTria( ptL1, vtL1, 100.0, trTria2, ptI3, ptI32, false) ;
         if ( PointInTria( ptI3, trTria1) != PTT_OUT  &&  PointInTria( ptI32, trTria1) != PTT_OUT) {
            ptInt = ( ptI3 + ptI32) / 2 ;
            ptInt2 = ( ptI3 + ptI32) / 2 ;
            return ITTTS_EDGE_EDGE_PNT ;
         }
         else
            return ITTTS_NO ;
      }      
   }
   return nTriaTriaIntersCases[nRes1][nRes2] ;
}

//----------------------------------------------------------------------------
int
IntersCoplanarTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, TRIA3DVECTOR& vTria)
{
  // Se i tre vertici del secondo triangolo stanno tutti a destra di almeno un lato del primo, non c'è intersezione
   for ( int i = 0 ; i < 3 ; ++ i) {
      Vector3d vtSide = trTria1.GetP( ( i + 1) % 3) - trTria1.GetP( i) ;
      Vector3d vtSN = vtSide ^ trTria1.GetN() ;
      vtSN.Normalize() ;
      if ( ( trTria2.GetP( 0) - trTria1.GetP( i)) * vtSN > - EPS_TRIA_H  &&
           ( trTria2.GetP( 1) - trTria1.GetP( i)) * vtSN > - EPS_TRIA_H  &&
           ( trTria2.GetP( 2) - trTria1.GetP( i)) * vtSN > - EPS_TRIA_H)
         return ITTT_NO ;
   }
  // Se i tre vertici del primo triangolo stanno tutti a destra di almeno un lato del secondo, non c'è intersezione
   for ( int i = 0 ; i < 3 ; ++ i) {
      Vector3d vtSide = trTria2.GetP( ( i + 1) % 3) - trTria2.GetP( i) ;
      Vector3d vtSN = vtSide ^ trTria2.GetN() ;
      vtSN.Normalize() ;
      if ( ( trTria1.GetP( 0) - trTria2.GetP( i)) * vtSN > - EPS_TRIA_H  &&
           ( trTria1.GetP( 1) - trTria2.GetP( i)) * vtSN > - EPS_TRIA_H  &&
           ( trTria1.GetP( 2) - trTria2.GetP( i)) * vtSN > - EPS_TRIA_H)
         return ITTT_NO ;
   }

  // determino il piano medio dei due triangoli
   Plane3d plMed ;
   plMed.Set( ( trTria2.GetCentroid() + trTria2.GetCentroid()) / 2, trTria1.GetN() + trTria2.GetN()) ;

  // creo la regione equivalente al primo triangolo
   SurfFlatRegion sfrTria1 ;
   PtrOwner<CurveComposite> pCcTria1( CreateBasicCurveComposite()) ;
   if ( IsNull( pCcTria1))
      return ITTT_NO ;
   pCcTria1->AddPoint( ProjectPointOnPlane( trTria1.GetP( 0), plMed)) ;
   pCcTria1->AddLine( ProjectPointOnPlane( trTria1.GetP( 1), plMed)) ;
   pCcTria1->AddLine( ProjectPointOnPlane( trTria1.GetP( 2), plMed)) ;
   pCcTria1->Close() ;
   if ( ! sfrTria1.AddExtLoop( Release( pCcTria1)))
      return ITTT_NO ;

  // creo la regione equivalente al secondo triangolo
   SurfFlatRegion sfrTria2 ;
   PtrOwner<CurveComposite> pCcTria2( CreateBasicCurveComposite()) ;
   if ( IsNull( pCcTria2))
      return ITTT_NO ;
   pCcTria2->AddPoint( ProjectPointOnPlane( trTria2.GetP( 0), plMed)) ;
   pCcTria2->AddLine( ProjectPointOnPlane( trTria2.GetP( 1), plMed)) ;
   pCcTria2->AddLine( ProjectPointOnPlane( trTria2.GetP( 2), plMed)) ;
   pCcTria2->Close() ;
   if ( ! sfrTria2.AddExtLoop( Release( pCcTria2)))
      return ITTT_NO ;
   if ( sfrTria1.GetNormVersor() * sfrTria2.GetNormVersor() < 0)
      sfrTria2.Invert() ;

  // calcolo l'intersezione tra le due regioni
   if ( ! sfrTria1.Intersect( sfrTria2)  ||  ! sfrTria1.IsValid())
      return ITTT_NO ;

  // recupero il contorno esterno del risultato come polilinea
   PolyLine PL ;
   PtrOwner<ICurve> pLoop( sfrTria1.GetLoop( 0, 0)) ;
   if ( IsNull( pLoop)  ||  pLoop->GetType() != CRV_COMPO)
      return ITTT_NO ;
   ICurveComposite* pCoLoop = GetBasicCurveComposite( pLoop) ;
   for ( int i = 0 ; i < pCoLoop->GetCurveCount() ; ++ i) {
      const ICurve* pCrv = pCoLoop->GetCurve( i) ;
      if ( i == 0) {
         Point3d ptS ; pCrv->GetStartPoint( ptS) ;
         PL.AddUPoint( i, ptS) ;
      }
      Point3d ptE ; pCrv->GetEndPoint( ptE) ;
      PL.AddUPoint( i+1, ptE) ;
   }

  // eseguo una triangolazione del contorno chiuso
   PNTVECTOR vPnt ;
   INTVECTOR vTrVert ;
   Triangulate Tri ;
   if ( ! Tri.Make( PL, vPnt, vTrVert))
      return ITTT_NO ;
   int nTrVert = int( vTrVert.size()) / 3 ;
   for ( int i = 0 ; i < nTrVert ; ++i) {
      Triangle3d Tria ;
      Tria.Set( vPnt[vTrVert[3*i]], vPnt[vTrVert[3*i+1]], vPnt[vTrVert[3*i+2]]) ;
      if ( Tria.Validate( true))
         vTria.emplace_back( Tria) ;
   }

   return ITTT_OVERLAPS ;
}