//----------------------------------------------------------------------------
//  EgalTech 2014-2015
//----------------------------------------------------------------------------
// File : PolygonPlane.cpp          Data : 23.02.15          Versione : 1.6b7
// Contenuto : Implementazione della classe PolygonPlane.
//             Calcolo della normale e delle aree con il metodo di Newell.
//
//
// Modifiche : 12.08.14 DS Creazione modulo.
//             23.02.15 DS Aggiunta gestione centro geometrico (centroid).
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "PolygonPlane.h"

//----------------------------------------------------------------------------
void
PolygonPlane::AddPoint( const Point3d& ptP)
{
  // se è il primo punto (parto da -1 perchè verrà contato alla chiusura)
   if ( m_nPntNbr == -1) {
     // inizializzazioni
      m_dLenN = -1 ;
      m_vtN = V_NULL ;
      m_ptMid = ORIG ;
      m_dSXy = 0 ; m_dSXz = 0 ;
      m_dSYz = 0 ; m_dSYx = 0 ;
      m_dSZx = 0 ; m_dSZy = 0 ;
     // salvo il punto come primo (per verificare chiusura alla fine)
      m_ptFirst = ptP ;
   }
  // dal secondo punto posso cominciare ad eseguire le somme parziali
   else {
     // Compute normal as being proportional to projected areas of polygon onto the yz,
     // xz, and xy planes. Also compute centroid as representative point on the plane
      m_vtN.x += ( m_ptLast.y - ptP.y) * ( m_ptLast.z + ptP.z) ;   // projection on yz
      m_vtN.y += ( m_ptLast.z - ptP.z) * ( m_ptLast.x + ptP.x) ;   // projection on xz
      m_vtN.z += ( m_ptLast.x - ptP.x) * ( m_ptLast.y + ptP.y) ;   // projection on xy
     // Very approximate center
      m_ptMid += ptP ;
     // First moment of Area (momento statico)
      double dTmpX = ( m_ptLast.y * ptP.z - ptP.y * m_ptLast.z) ;
      m_dSXy += ( m_ptLast.y + ptP.y) * dTmpX ;
      m_dSXz += ( m_ptLast.z + ptP.z) * dTmpX ;
      double dTmpY = ( m_ptLast.z * ptP.x - ptP.z * m_ptLast.x) ;
      m_dSYz += ( m_ptLast.z + ptP.z) * dTmpY ;
      m_dSYx += ( m_ptLast.x + ptP.x) * dTmpY ;
      double dTmpZ = ( m_ptLast.x * ptP.y - ptP.x * m_ptLast.y) ;
      m_dSZx += ( m_ptLast.x + ptP.x) * dTmpZ ;
      m_dSZy += ( m_ptLast.y + ptP.y) * dTmpZ ;
   }
  // salvo punto e incremento numero di punti
   m_ptLast = ptP ;
   ++ m_nPntNbr ;
}

//----------------------------------------------------------------------------
bool
PolygonPlane::Finalize( void)
{
  // almeno 3 punti (il triangolo è il poligono con minimo numero di lati)
   if ( m_nPntNbr + 1 < 3)
      return false ;
  // se il poligono non è stato chiuso, aggiungo calcolo per ultimo lato
   if ( ! AreSamePointExact( m_ptFirst, m_ptLast)) {
     // aggiungo il primo punto per far eseguire i conti sul lato di chiusura
      AddPoint( m_ptFirst) ;
   }
  // se non effettuato, eseguo il calcolo finale
   if ( m_dLenN < 0) {
     // lunghezza della normale (doppio dell'area del poligono)
      m_dLenN = m_vtN.Len() ;
      if ( m_dLenN < SQ_EPS_SMALL)
         return false ;
     // normalizzo
      m_vtN /= m_dLenN ;
     // sistemo baricentro
      m_ptMid /= m_nPntNbr ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
PolygonPlane::GetCentroid( Point3d& ptCen)
{
  // verifico completamento conti
   if ( ! Finalize())
      return false ;
  // assegno i dati del centro geometrico
   if ( abs( m_vtN.z) > abs( m_vtN.x)  &&
        abs( m_vtN.z) > abs( m_vtN.y)) {
     // calcoli nel piano xy perpendicolare a Z
      ptCen.x = m_dSZx / ( 3 * m_vtN.z * m_dLenN) ;
      ptCen.y = m_dSZy / ( 3 * m_vtN.z * m_dLenN) ;
      ptCen.z = (( m_ptMid.x - ptCen.x) * m_vtN.x + ( m_ptMid.y - ptCen.y) * m_vtN.y + m_ptMid.z * m_vtN.z) / m_vtN.z ;
   }
   else if ( abs( m_vtN.x) > abs( m_vtN.y)) {
     // calcoli nel piano yz perpendicolare a X
      ptCen.y = m_dSXy / ( 3 * m_vtN.x * m_dLenN) ;
      ptCen.z = m_dSXz / ( 3 * m_vtN.x * m_dLenN) ;
      ptCen.x = ( m_ptMid.x * m_vtN.x + ( m_ptMid.y - ptCen.y) * m_vtN.y + ( m_ptMid.z - ptCen.z) * m_vtN.z) / m_vtN.x ;
   }
   else {
     // calcoli nel piano zx perpendicolare a Y
      ptCen.z = m_dSYz / ( 3 * m_vtN.y * m_dLenN) ;
      ptCen.x = m_dSYx / ( 3 * m_vtN.y * m_dLenN) ;
      ptCen.y = (( m_ptMid.x - ptCen.x) * m_vtN.x + m_ptMid.y * m_vtN.y + ( m_ptMid.z - ptCen.z) * m_vtN.z) / m_vtN.y ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
PolygonPlane::GetNormal( Vector3d& vtN)
{
  // verifico completamento conti
   if ( ! Finalize())
      return false ;
  // assegno i dati del centro
   vtN = m_vtN ;
   return true ;
}

//----------------------------------------------------------------------------
bool
PolygonPlane::GetPlane( Plane3d& plPlane)
{
  // verifico completamento conti
   if ( ! Finalize())
      return false ;
  // assegno i dati al piano
   return plPlane.Set( m_ptMid, m_vtN) ;
}

//----------------------------------------------------------------------------
bool
PolygonPlane::GetArea( double& dArea)
{
  // verifico completamento conti
   if ( ! Finalize())
      return false ;
  // assegno l'area
   dArea = 0.5 * m_dLenN ;
   return true ;
}