From b3abcf73c367d574bb2e1ce9a0da3460578e83c3 Mon Sep 17 00:00:00 2001 From: Riccardo Elitropi Date: Tue, 4 Nov 2025 16:02:04 +0100 Subject: [PATCH] EgtGeomKernel 2.7k2 : - aggiunta funzione per il calcolo della curvatura in un vertice per una superficie TriMesh. --- EgtGeomKernel.rc | Bin 11718 -> 11718 bytes SurfTriMesh.h | 4 ++ SurfTriMeshFaceting.cpp | 151 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 155 insertions(+) diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 97503c860207c7daf927c438ef4cfa4edb1f79e9..f3c849c992d94942feee72fb8d1101a54f5f4aed 100644 GIT binary patch delta 265 zcmX>WeJpyzI%aiGhCGIJh8%`ehGHOAV9;O)W+-7u1ma8}FBB-731Uyas3kpFPd#gM z9Lqc=Mx)L5`OY&hX`jGK#;zcG^Fwu6)J>CBt_NlipBAM6KDupe@OZchZce6ybVG-gbx%@?)4G64W1 CWk|>X delta 260 zcmX>WeJpyzI%aKGhCGIJh8%`WhGK>c1_cHUhESkLCPN8BDnkm89}HwAf_cT0FKS6| z{=m$^%xJjzHs5)s$u(R$lRrqRF&a+ZC?LLBKq-!K^8{`$h_o$`0t*`mtzZFbKoS8e z?2(2l#9~eesm?k$`JT?a$* #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 ; +}