diff --git a/PointsPCA.cpp b/PointsPCA.cpp index ec7874c..7ca4bfd 100644 --- a/PointsPCA.cpp +++ b/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 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 ; diff --git a/PointsPCA.h b/PointsPCA.h index f628a87..bd33ab7 100644 --- a/PointsPCA.h +++ b/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 } ; \ No newline at end of file