diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 97503c8..f3c849c 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/SurfTriMesh.h b/SurfTriMesh.h index 8c5a650..3e4345e 100644 --- a/SurfTriMesh.h +++ b/SurfTriMesh.h @@ -322,6 +322,10 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW bool GetEdge( int nInd, int& nV1, int& nV2, int& nFl, int& nFr, double& dAng) const override ; bool GetEdge( int nInd, Point3d& ptP1, Point3d& ptP2, double& dAng) const override ; bool GetEdges( ICURVEPOVECTOR& vpCurve) const override ; + bool GetCurvature( int nV, + double& dMinK, Vector3d& vtMinK, double& dMaxK, Vector3d& vtMaxK, bool& bPlanar, Vector3d& vtNorm, + double* dGaussK = nullptr, double* dMeanK = nullptr, + double* dMaxNormK = nullptr, double* dMinNormK = nullptr) const override ; bool Cut( const Plane3d& plPlane, bool bSaveOnEq) override ; bool GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) override ; bool Add( const ISurfTriMesh& Other) override ; diff --git a/SurfTriMeshFaceting.cpp b/SurfTriMeshFaceting.cpp index bbb1529..2fc7347 100644 --- a/SurfTriMeshFaceting.cpp +++ b/SurfTriMeshFaceting.cpp @@ -21,6 +21,8 @@ #include #include #include +#define EIGEN_NO_IO +#include "/EgtDev/Extern/Eigen/Dense" using namespace std ; @@ -1015,3 +1017,152 @@ SurfTriMesh::GetEdges( ICURVEPOVECTOR& vpCurve) const } return true ; } + +//----------------------------------------------------------------------------- +// Funzione per il calcolo della curvatura massima e minima in un vertice della superficie TriMesh +// vtMinK : versore direzione curvatura minima ( dMinK è il suo modulo con segno) +// vtMaxK : versore direzione curvatura massima ( dMaxK è il suo modulo con segno) +// bPlanar : Flag per indicare le superficie localmente piana +// vtNorm : versore normale del piano tangente alla supericie nel vertice +// dGaussK : valore della curvatura gaussiana con segno ( da associarsi a vtNorm) +// dMeanK : valore della curvatura media con segno ( da associarsi a vtNorm) +// dMaxNormK : valore della curvatura normale massima con segno ( da associarsi a vtNorm) +// dMinNormK : valore della curvatura normale minima con segno ( da associarsi a vtNorm) +bool +SurfTriMesh::GetCurvature( int nV, + double& dMinK, Vector3d& vtMinK, double& dMaxK, Vector3d& vtMaxK, bool& bPlanar, Vector3d& vtNorm, + double* dGaussK, double* dMeanK, double* dMaxNormK, double* dMinNormK) const +{ + // Controllo la validità della TriMesh e del Vertice + if ( ! IsValid()) + return false ; + Point3d ptCurr ; + if ( ! GetVertex( nV, ptCurr)) + return false ; + bPlanar = false ; + + // Recupero tutti i vertici attorno al vertice corrente ( se presenza di lato libero, allora errore) + INTVECTOR vT ; + bool bCirc ; + if ( ! GetAllTriaAroundVertex( nV, vT, bCirc) || ! bCirc) + return false ; + + // Calcolo la normale del vertice pesata mediante angolo sotteso + // [Meyer et al. (2003) – Discrete Differential-Geometry Operators for Triangulated 2-Manifolds + // Section 5.4.2. Tip-Angle Weights.] + // [https://www.cs.cmu.edu/~kmcrane/Projects/DDG/paper.pdf] + INTSET setIndTriaNeightbors ; + for ( const int& nT : vT) { + // Recupero le direzioni delle semirette per l'angolo al vertice corrente + Triangle3d Tria ; + int nIdVert[3] ; + if ( ! GetTriangle( nT, Tria) || ! GetTriangle( nT, nIdVert)) + return false ; + Vector3d vtDir0 = V_NULL, vtDir1 = V_NULL ; + for ( int i = 0 ; i < 3 ; ++ i) { + if ( AreSamePointApprox( Tria.GetP( i), ptCurr)) { + vtDir0 = Tria.GetP( ( i + 1) % 3) - Tria.GetP( i) ; + vtDir0.Normalize() ; + vtDir1 = Tria.GetP( ( i + 2) % 3) - Tria.GetP( i) ; + vtDir1.Normalize() ; + break ; + } + } + for ( int i = 0 ; i < 3 ; ++ i) { + if ( nV == nIdVert[i]) { + setIndTriaNeightbors.insert( nIdVert[( i + 1) % 3]) ; + setIndTriaNeightbors.insert( nIdVert[( i + 2) % 3]) ; + break ; + } + } + // Calcolo l'angolo sotteso + double dCosTheta = max( -1., min( 1., ( vtDir0 * vtDir1))) ; + double dTheta = acos( dCosTheta) ; + // Aggiorno il contributo della normale al vertice + vtNorm += ( dTheta * Tria.GetN()) ; + } + vtNorm.Normalize() ; + + // [ESTIMATING CURVATURE ON TRIANGULAR MESHES, cap. 2.1. Fitting Methods, + // par. 2.1.2. Quadric Fitting] + // [https://people.eecs.berkeley.edu/~jrs/meshpapers/GatzkeGrimm.pdf] + + // Definisco il piano tangente al vertice corrente + Plane3d plTan ; + if ( ! plTan.Set( ptCurr, vtNorm)) + return false ; + + // Proietto i punti Neightbors(1) sul piano tangete ( senza ripeterli) + BIPNTVECTOR vPtPtProj ; vPtPtProj.reserve( setIndTriaNeightbors.size()) ; + for ( auto nIter = setIndTriaNeightbors.begin() ; nIter != setIndTriaNeightbors.end() ; ++ nIter) { + Point3d ptNeightbors ; + if ( ! GetVertex( *nIter, ptNeightbors)) + return false ; + vPtPtProj.emplace_back( make_pair( ptNeightbors, ProjectPointOnPlane( ptNeightbors, plTan))) ; + } + + // Recupero due versori perpendicolari nel piano definito + Vector3d vtTan1 = V_NULL ; + if ( abs( vtNorm.x) < 1./64. && abs( vtNorm.y) < 1./64.) + vtTan1 = Y_AX ^ vtNorm ; + else + vtTan1 = Z_AX ^ vtNorm ; + vtTan1.Normalize() ; + Vector3d vtTan2 = vtNorm ^ vtTan1 ; + vtTan2.Normalize() ; + + // Sistema da risolvere mediante minimi quadrati : z(u,v) = Au^2 + Buv + Cv^2 + // Questo sistema è molto più semplice e robusto rispetto a z(u,v) = Au^2 + Buv + Cv^2 + Du + Ev + F + // per il fatto che ora le coordinate sono in locale al piano tangente alla superficie ( mettendo + // quindi a 0 i coefficienti del primo ordine e il termine noto F) + // Definizione della matrice A(Nx3) + const int nPts = int( vPtPtProj.size()) ; + Eigen::MatrixXd mat_A( nPts, 3) ; + Eigen::VectorXd vec_b( nPts) ; + for ( int i = 0 ; i < nPts ; ++ i) { + Vector3d vtEdgeProj = vPtPtProj[i].second - ptCurr ; + Vector3d vtEdge = ( vPtPtProj[i].first - vPtPtProj[i].second) ; + double dU = vtEdgeProj * vtTan1 ; + double dV = vtEdgeProj * vtTan2 ; + double dW = vtEdge * vtNorm ; + mat_A( i, 0) = dU * dU ; + mat_A( i, 1) = dU * dV ; + mat_A( i, 2) = dV * dV ; + vec_b( i) = dW ; + } + // Risoluzione del sistema + Eigen::VectorXd vec_x = mat_A.colPivHouseholderQr().solve( vec_b) ; + + // Costruzione della matrice hessiana + Eigen::Matrix2d mat_H ; + mat_H << 2. * vec_x( 0), vec_x( 1), + vec_x( 1), 2. * vec_x( 2) ; + + // Calcolo gli autovalori ( quindi le curvature principali) e restituisco i risultati + Eigen::SelfAdjointEigenSolver Solver( mat_H) ; + dMinK = Solver.eigenvalues()( 0) ; + dMaxK = Solver.eigenvalues()( 1) ; + Eigen::Vector2d dir_Min = Solver.eigenvectors().col( 0) ; + Eigen::Vector2d dir_Max = Solver.eigenvectors().col( 1) ; + vtMinK = dir_Min( 0) * vtTan1 + dir_Min( 1) * vtTan2 ; + vtMaxK = dir_Max( 0) * vtTan1 + dir_Max( 1) * vtTan2 ; + bPlanar = ( abs( dMinK) < EPS_SMALL && abs( dMaxK) < EPS_SMALL) ; + + // Se rischiesto calcolo della curvatura Gaussiana + if ( dGaussK != nullptr) + *dGaussK = ( dMinK * dMaxK) ; + + // Se richiesto calcolo della curvatura Media + if ( dMeanK != nullptr) + *dMeanK = ( dMinK + dMaxK) / 2. ; + + // Se richiesto calcolo della curvatura normale lungo la direzione massima + if ( dMaxNormK != nullptr) + *dMaxNormK = dir_Max.transpose() * mat_H * dir_Max ; + + // Se richiesto calcolo della curvatura normale lungo la direzione minima + if ( dMinNormK != nullptr) + *dMinNormK = dir_Min.transpose() * mat_H * dir_Min ; + + return true ; +}