EgtGeomKernel :

- correzione calcolo PCA per punti (aumentata precisione).

PointsPCA.cpp

@@ -14,36 +14,28 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "PointsPCA.h"
+#include "/EgtDev/Extern/Eigen/Dense"
 
 
 //----------------------------------------------------------------------------
 PointsPCA::PointsPCA( void)
 {
-  // azzero numero punti
-   m_dTotW = 0 ;
-  // azzero il baricentro
-   m_ptCen.Set( 0, 0, 0) ;
-  // azzero matrice di covarianza
-   m_CovMat.setZero() ;
   // inizializzo il rank ad un valore assurdo
    m_nRank = - 1 ;
+  // inizializzo baricentro
+   m_ptCen.Set( 0, 0, 0) ;
+  // inizializzo il peso totale
+   m_dTotW = 0 ;
+  // riservo memoria per la matrice di punti e pesi
+   m_vPntW.reserve( 128) ;
 }
 
 //----------------------------------------------------------------------------
 void
 PointsPCA::AddPoint( const Point3d& ptP, double dW)
 {
-  // incremento numero punti
-   m_dTotW += dW ;
-  // aggiorno il baricentro
-   m_ptCen += dW * ptP ;
-  // aggiorno la matrice di covarianza (solo triangolo superiore perchè simmetrica)
-   m_CovMat(0,0) += dW * ( ptP.x * ptP.x) ;
-   m_CovMat(1,1) += dW * ( ptP.y * ptP.y) ;
-   m_CovMat(2,2) += dW * ( ptP.z * ptP.z) ;
-   m_CovMat(0,1) += dW * ( ptP.x * ptP.y) ;
-   m_CovMat(0,2) += dW * ( ptP.x * ptP.z) ;
-   m_CovMat(1,2) += dW * ( ptP.y * ptP.z) ;
+  // salvo i dati
+   m_vPntW.emplace_back( ptP, dW) ;
 }
 
 //----------------------------------------------------------------------------
@@ -58,29 +50,47 @@ PointsPCA::Finalize( void)
    m_nRank = 0 ;
 
   // se non sono stati assegnati punti, esco
-   if ( m_dTotW < EPS_ZERO)
+   if ( m_vPntW.empty())
       return false ;
 
-  // fattore di scala per numero di punti
+  // calcolo del peso totale
+   for ( const auto& PntW : m_vPntW)
+      m_dTotW += PntW.second ;
+   if ( m_dTotW < EPS_ZERO)
+      return false ;
+  // fattore di scala
    double dScale = 1 / m_dTotW ;
 
   // calcolo del baricentro
+   for ( const auto& PntW : m_vPntW)
+      m_ptCen += PntW.second * PntW.first ;
    m_ptCen *= dScale ;
 
-  // completo la matrice di covarianza
-   m_CovMat(0,0) = m_CovMat(0,0) * dScale - m_ptCen.x * m_ptCen.x ;
-   m_CovMat(1,1) = m_CovMat(1,1) * dScale - m_ptCen.y * m_ptCen.y ;
-   m_CovMat(2,2) = m_CovMat(2,2) * dScale - m_ptCen.z * m_ptCen.z ;
-   m_CovMat(0,1) = m_CovMat(0,1) * dScale - m_ptCen.x * m_ptCen.y ;
-   m_CovMat(0,2) = m_CovMat(0,2) * dScale - m_ptCen.x * m_ptCen.z ;
-   m_CovMat(1,2) = m_CovMat(1,2) * dScale - m_ptCen.y * m_ptCen.z ;
-   m_CovMat(1,0) = m_CovMat(0,1) ;
-   m_CovMat(2,0) = m_CovMat(0,2) ;
-   m_CovMat(2,1) = m_CovMat(1,2) ;
+  // matrice di covarianza
+   Eigen::Matrix3d CovMat ;
+   CovMat.setZero() ;
+   for ( const auto& PntW : m_vPntW) {
+      Point3d ptP( PntW.first.x - m_ptCen.x, PntW.first.y - m_ptCen.y, PntW.first.z - m_ptCen.z) ;
+      CovMat(0,0) += PntW.second * ( ptP.x * ptP.x) ;
+      CovMat(1,1) += PntW.second * ( ptP.y * ptP.y) ;
+      CovMat(2,2) += PntW.second * ( ptP.z * ptP.z) ;
+      CovMat(0,1) += PntW.second * ( ptP.x * ptP.y) ;
+      CovMat(0,2) += PntW.second * ( ptP.x * ptP.z) ;
+      CovMat(1,2) += PntW.second * ( ptP.y * ptP.z) ;
+   }
+   CovMat(0,0) = CovMat(0,0) * dScale ;
+   CovMat(1,1) = CovMat(1,1) * dScale ;
+   CovMat(2,2) = CovMat(2,2) * dScale ;
+   CovMat(0,1) = CovMat(0,1) * dScale ;
+   CovMat(0,2) = CovMat(0,2) * dScale ;
+   CovMat(1,2) = CovMat(1,2) * dScale ;
+   CovMat(1,0) = CovMat(0,1) ;
+   CovMat(2,0) = CovMat(0,2) ;
+   CovMat(2,1) = CovMat(1,2) ;
 
-  // calcolo gli autovalori e autovettori (essendo matrice 3x3 uso metodo diretto)
+  // calcolo gli autovalori e autovettori
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> es ;
-   es.compute( m_CovMat) ;        // non usare computeDirect : errore nell'ordine degli autovettori
+   es.compute( CovMat) ;        // non usare computeDirect : errore nell'ordine degli autovettori
    if ( es.info() == Eigen::NoConvergence)
       return false ;
 

PointsPCA.h

@@ -13,8 +13,7 @@
 
 #pragma once
 
-#include "/EgtDev/Include/EGkPoint3d.h"
-#include "/EgtDev/Extern/Eigen/Dense"
+#include "/EgtDev/Include/EGkGeoCollection.h"
 
 //----------------------------------------------------------------------------
 class PointsPCA
@@ -33,9 +32,9 @@ class PointsPCA
       static const int MAX_RANK = 3 ;
 
    private :
-      double          m_dTotW ;            // peso totale (se pesi tutti unitari, allora è numero punti)
+      PNTUVECTOR      m_vPntW ;            // vettore dei punti con i loro pesi
       Point3d         m_ptCen ;            // baricentro
-      Eigen::Matrix3d m_CovMat ;           // matrice di covarianza
+      double          m_dTotW ;            // peso totale
       int             m_nRank ;            // numero delle componenti principali (MAX_RANK = 3)
       Vector3d        m_vtPC[MAX_RANK] ;   // direzioni delle componenti principali
 } ;