From 87e3800d1a3f8f492f02eed58fc6e834b0842124 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Sun, 12 Jan 2014 22:48:09 +0000 Subject: [PATCH] =?UTF-8?q?EgtGeomKernel=201.5a2=20:=20Aggiunto=20calcolo?= =?UTF-8?q?=20singolarit=C3=A0=20in=20Curve=20Bezier.=20Migliorato=20calco?= =?UTF-8?q?lo=20distanza=20punto-curva=20di=20Bezier.=20Aggiunta=20classe?= =?UTF-8?q?=20PolynomialPoint3d.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CurveBezier.cpp | 222 ++++++++++++++++++++++- CurveBezier.h | 9 +- DistPointCrvBezier.cpp | 321 +++++++++++++++++++++++++--------- EgtGeomKernel.rc | Bin 7666 -> 7650 bytes EgtGeomKernel.vcxproj | 10 +- EgtGeomKernel.vcxproj.filters | 18 +- GdbExecutor.cpp | 45 +++-- GdbExecutor.h | 2 +- OutScl.cpp | 96 +++++++++- OutScl.h | 3 +- PolynomialPoint3d.cpp | 225 ++++++++++++++++++++++++ PolynomialPoint3d.h | 193 ++++++++++++++++++++ stdafx.h | 1 + 13 files changed, 1035 insertions(+), 110 deletions(-) create mode 100644 PolynomialPoint3d.cpp create mode 100644 PolynomialPoint3d.h diff --git a/CurveBezier.cpp b/CurveBezier.cpp index 03b998a..e6cece4 100644 --- a/CurveBezier.cpp +++ b/CurveBezier.cpp @@ -1,7 +1,7 @@ //---------------------------------------------------------------------------- -// EgalTech 2013-2013 +// EgalTech 2013-2014 //---------------------------------------------------------------------------- -// File : CurveBezier.cpp Data : 22.11.13 Versione : 1.3a1 +// File : CurveBezier.cpp Data : 12.01.14 Versione : 1.5a2 // Contenuto : Implementazione della classe Curva di Bezier. // // @@ -17,8 +17,10 @@ #include "CurveBezier.h" #include "DistPointLine.h" #include "GeoObjFactory.h" -#include "\EgtDev\Include\EGnStringUtils.h" +#include "PolynomialPoint3d.h" #include "\EgtDev\Include\EGkCurveArc.h" +#include "\EgtDev\Include\ENkPolynomial.h" +#include "\EgtDev\Include\EGnStringUtils.h" #include using namespace std ; @@ -28,7 +30,8 @@ GEOOBJ_REGISTER( CRV_BEZ, "C_BEZ", CurveBezier) ; //---------------------------------------------------------------------------- CurveBezier::CurveBezier( void) - : m_nStatus( TO_VERIFY), m_nDeg( 0), m_bRat( false), m_nDimArr( 0), m_aPtCtrl( nullptr), m_aWeCtrl( nullptr) + : m_nStatus( TO_VERIFY), m_nDeg( 0), m_bRat( false), m_dParSing( -2), + m_nDimArr( 0), m_aPtCtrl( nullptr), m_aWeCtrl( nullptr) { } @@ -89,6 +92,9 @@ CurveBezier::Init( int nDeg, bool bIsRational) // salvo flag di razionale m_bRat = bIsRational ; + // dichiaro non si è ancora fatta analisi presenza singolarità + m_dParSing = - 2 ; + // pulisco, assegnando 0 ai punti e 1 ai pesi memset( m_aPtCtrl, 0, sizeof( Point3d) * ( m_nDeg + 1)) ; if ( m_bRat) { @@ -115,6 +121,9 @@ CurveBezier::SetControlPoint( int nInd, const Point3d& ptCtrl) if ( m_bRat) m_aWeCtrl[nInd] = 1 ; + // annullo analisi presenza singolarità + m_dParSing = - 2 ; + return true ; } @@ -134,6 +143,9 @@ CurveBezier::SetControlPoint( int nInd, const Point3d& ptCtrl, double dW) m_aPtCtrl[nInd] = ptCtrl ; m_aWeCtrl[nInd] = dW ; + // annullo analisi presenza singolarità + m_dParSing = - 2 ; + return true ; } @@ -162,7 +174,7 @@ CurveBezier::SetFromArc( const ICurveArc& crArc) vtDir = ptMed - crArc.GetCenter() ; ptNew = crArc.GetCenter() + vtDir / ( dCosAhalf * dCosAhalf) ; - Init( 2, true ) ; + Init( 2, true) ; SetControlPoint( 0, ptStart, 1) ; SetControlPoint( 1, ptNew, dCosAhalf) ; SetControlPoint( 2, ptEnd, 1) ; @@ -487,6 +499,206 @@ CurveBezier::IncreaseBernsteinOneDegree( double dU, int nDeg, double dBern[]) co } +//---------------------------------------------------------------------------- +int +CurveBezier::GetSingularParam( double& dPar) const +{ + // se da calcolare, lo faccio + if ( m_dParSing < -1.5) + CalcSingularParam() ; + + // se esiste singolarità + if ( m_dParSing > 0 && m_dParSing < 1) { + dPar = m_dParSing ; + return 1 ; + } + // altrimenti, non esiste + else + return 0 ; +} + +//---------------------------------------------------------------------------- +bool +CurveBezier::CalcSingularParam( void) const +{ + // non si considerano singolarità agli estremi ( non creano problemi) + + // le curve di grado 1 non possono avere singolarità + if ( m_nDeg <= 1) { + m_dParSing = -1 ; + return true ; + } + + // le curve di grado 2 possono avere singolarità solo se i punti sono allineati con ordine P0 P2 P1 + if ( m_nDeg == 2) { + Vector3d vtV1 = m_aPtCtrl[1] - m_aPtCtrl[0] ; + Vector3d vtV2 = m_aPtCtrl[2] - m_aPtCtrl[1] ; + if ( ( vtV1 * vtV2) > - EPS_ZERO) { + m_dParSing = -1 ; + return true ; + } + } + + // calcolo della derivata o del suo numeratore se curva razionale + PolynomialPoint3d pol3P ; + + // se polinomiale + if ( ! m_bRat) { + // derivata come curva di Bezier delle differenze tra i punti di controllo + CurveBezier cbDeriv ; + cbDeriv.Init( m_nDeg - 1, false) ; + for ( int i = 0 ; i < m_nDeg ; ++ i) + cbDeriv.SetControlPoint( i, m_nDeg * ( m_aPtCtrl[i+1] + ( - m_aPtCtrl[i]))) ; + // trasformo nella forma monomia + cbDeriv.ToPowerBase( pol3P) ; + } + // altrimenti razionale + else { + // numeratore e denominatore in forma monomia + PolynomialPoint3d pol3Num ; + Polynomial polDen ; + ToPowerBase( pol3Num, polDen) ; + // derivata del numeratore + PolynomialPoint3d pol3NumD ; + pol3NumD.Derive( pol3Num) ; + // derivata del denominatore + Polynomial polDenD ; + polDenD.Derive( polDen) ; + // numeratore della derivata complessiva ( NumD * Den - Num * DenD) + pol3P = pol3NumD * polDen - pol3Num * polDenD ; + pol3P.AdjustDegree() ; + } + + // calcolo gli zeri della componente più importante della derivata + DBLVECTOR vdRoot ; + int nZ = pol3P.FindMainComponentRoots( vdRoot) ; + // scarto gli zeri fuori dall'intervallo aperto 0-1 e quelli ripetuti + nZ = FilterMultipleAndOutOfRangeRoots( vdRoot, 0 + EPS_ZERO, 1 - EPS_ZERO, EPS_ZERO) ; + + // verifico se in corrispondenza di qualche zero si annulla la derivata + Point3d ptPos ; + Vector3d vtDer1 ; + for ( int i = 0 ; i < nZ ; ++ i) { + GetPointD1D2( vdRoot[i], FROM_MINUS, ptPos, &vtDer1) ; + if ( vtDer1.IsZero()) { + m_dParSing = vdRoot[i] ; + LOG_DBG_INFO( GetEGkLogger(), "INFO : Found Singularity in CurveBezier") + return true ; + } + } + + // non ci sono singolarità + m_dParSing = -1 ; + return true ; +} + +//---------------------------------------------------------------------------- +// i coefficienti sono ordinati dalla potenza più bassa alla più alta +// se razionale, ritorna la sola forma monomiale del numeratore +//---------------------------------------------------------------------------- +bool +CurveBezier::ToPowerBase( PolynomialPoint3d& pol3P) const +{ + // creo un vettore di punti ausiliario (pulito e di dimensioni adeguate) + PNTVECTOR vPow ; + vPow.clear() ; + vPow.reserve( m_nDeg + 1) ; + + // algoritmo di Sederberg (BYU) CAGD cap. 3 + // copio i punti di controllo + for ( int i = 0 ; i <= m_nDeg ; ++ i) + vPow.push_back( m_aPtCtrl[i]) ; + // se razionale, li moltiplico per i relativi pesi + if ( m_bRat) { + for ( int i = 0 ; i <= m_nDeg ; ++ i) + vPow[i] *= m_aWeCtrl[i] ; + } + // applico differenze successive, lasciando sul posto il risultato + for ( int i = 1 ; i <= m_nDeg ; ++ i) { + for ( int j = m_nDeg ; j >= i ; -- j) + vPow[j] += - vPow[j-1] ; + } + // calcolo i coefficienti binomiali + double dBinom[ MAXDEG+1] ; + dBinom[0] = 1 ; + for ( int i = 1 ; i <= m_nDeg ; ++ i) + dBinom[i] = 0 ; + for ( int i = 1 ; i <= m_nDeg ; ++ i) { + for ( int j = i ; j > 0 ; -- j) + dBinom[j] += dBinom[j-1] ; + } + // ottengo i coefficienti delle potenze + for ( int i = 1 ; i <= m_nDeg ; ++ i) + vPow[i] *= dBinom[i] ; + + // assegno il risultato al parametro di ritorno + pol3P.Set( m_nDeg, vPow) ; + + return true ; +} + +//---------------------------------------------------------------------------- +// i coefficienti sono ordinati dalla potenza più bassa alla più alta +// se polinomiale, la forma monomiale del denominatore è un 1 +//---------------------------------------------------------------------------- +bool +CurveBezier::ToPowerBase( PolynomialPoint3d& pol3Num, Polynomial& polDen) const +{ + // creo un vettore di punti ausiliario (pulito e di dimensioni adeguate) + PNTVECTOR vPow ; + vPow.reserve( m_nDeg + 1) ; + // creo un vettore di double ausiliario (pulito e di dimensioni adeguate) + DBLVECTOR dPow ; + if ( m_bRat) + dPow.reserve( m_nDeg + 1) ; + + // algoritmo di Sederberg (BYU) CAGD cap. 3 + // copio i punti di controllo + for ( int i = 0 ; i <= m_nDeg ; ++ i) + vPow.push_back( m_aPtCtrl[i]) ; + // se razionale + if ( m_bRat) { + // moltiplico i punti per i relativi pesi + for ( int i = 0 ; i <= m_nDeg ; ++ i) + vPow[i] *= m_aWeCtrl[i] ; + // copio i pesi + for ( int i = 0 ; i <= m_nDeg ; ++ i) + dPow.push_back( m_aWeCtrl[i]) ; + } + // applico differenze successive, lasciando sul posto il risultato + for ( int i = 1 ; i <= m_nDeg ; ++ i) { + for ( int j = m_nDeg ; j >= i ; -- j) { + vPow[j] += - vPow[j-1] ; + if ( m_bRat) + dPow[j] += - dPow[j-1] ; + } + } + // calcolo i coefficienti binomiali + double dBinom[ MAXDEG+1] ; + dBinom[0] = 1 ; + for ( int i = 1 ; i <= m_nDeg ; ++ i) + dBinom[i] = 0 ; + for ( int i = 1 ; i <= m_nDeg ; ++ i) { + for ( int j = i ; j > 0 ; -- j) + dBinom[j] += dBinom[j-1] ; + } + // ottengo i coefficienti delle potenze + for ( int i = 1 ; i <= m_nDeg ; ++ i) { + vPow[i] *= dBinom[i] ; + if ( m_bRat) + dPow[i] *= dBinom[i] ; + } + + // assegno i risultati ai parametri di ritorno + pol3Num.Set( m_nDeg, vPow) ; + if ( m_bRat) + polDen.Set( m_nDeg, dPow) ; + else + polDen.SetToConstant( 1) ; + + return true ; +} + //---------------------------------------------------------------------------- bool CurveBezier::GetControlPolygonLength( double& dLen) const diff --git a/CurveBezier.h b/CurveBezier.h index e37c378..94532c4 100644 --- a/CurveBezier.h +++ b/CurveBezier.h @@ -13,8 +13,10 @@ #pragma once -#include "/EgtDev/Include/EGkCurveBezier.h" #include "CurveAux.h" +#include "PolynomialPoint3d.h" +#include "/EgtDev/Include/EGkCurveBezier.h" +#include "/EgtDev/Include/ENkNumCollection.h" //---------------------------------------------------------------------------- class CurveBezier : public ICurveBezier @@ -73,6 +75,7 @@ class CurveBezier : public ICurveBezier virtual const Point3d& GetControlPoint( int nInd, bool* pbOk = NULL) const ; virtual double GetControlWeight( int nInd, bool* pbOk = NULL) const ; virtual bool GetControlPolygonLength( double& dLen) const ; + virtual int GetSingularParam( double& dPar) const ; public : CurveBezier( void) ; @@ -89,12 +92,15 @@ class CurveBezier : public ICurveBezier private : bool Validate( void) ; void IncreaseBernsteinOneDegree( double dU, int nDeg, double dBern[]) const ; + bool CalcSingularParam( void) const ; double GetSegmentLength( int nLev, double dU0, double dU1, double dU2, const Point3d& ptP0, const Point3d& ptP1, const Point3d& ptP2) const ; bool GetSegmentParam( double dLen, double& dCurrLen, double& dSegLen, double& dUIni, double& dUFin) const ; bool FlatOrSplit( int nLev, const CurveBezier& crvBez, double dParStart, double dParEnd, double dLinTol, double dAngTolDeg, PolyLine& PL) const ; + bool ToPowerBase( PolynomialPoint3d& pol3P) const ; + bool ToPowerBase( PolynomialPoint3d& pol3Num, Polynomial& polDen) const ; private : enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ; @@ -105,6 +111,7 @@ class CurveBezier : public ICurveBezier Status m_nStatus ; // stato int m_nDeg ; // grado bool m_bRat ; // flag di razionale/polinomiale + mutable double m_dParSing ; // eventuale parametro della singolarità (-1=no, -2=da calcolare) int m_nDimArr ; // dimensione dell'array dinamico di punti Point3d* m_aPtCtrl ; // array dei punti di controllo double* m_aWeCtrl ; // array dei pesi di controllo diff --git a/DistPointCrvBezier.cpp b/DistPointCrvBezier.cpp index f0583f4..0450611 100644 --- a/DistPointCrvBezier.cpp +++ b/DistPointCrvBezier.cpp @@ -13,88 +13,73 @@ //--------------------------- Include ---------------------------------------- #include "stdafx.h" +#include "DllMain.h" #include "DistPointCrvBezier.h" #include "DistPointLine.h" -#include "\EgtDev\Include\EgtTrace.h" - //---------------------------------------------------------------------------- -DistPointCrvBezier::DistPointCrvBezier( const Point3d& ptP, const ICurveBezier& CrvBez) +static const double LIN_TOL_APPROX = 1 ; +static const double ANG_TOL_APPROX_DEG = 45 ; + +//---------------------------------------------------------------------------- +struct MinDistCalc { + double dDist ; + double dPar ; + double dParMin ; + double dParMax ; + bool bParMinSing ; + bool bParMaxSing ; + Point3d ptQ ; + MinDistCalc( void) + : dDist( 0), dPar( 0), dParMin( 0), dParMax( 0), + bParMinSing( false), bParMaxSing( false), ptQ( 0, 0, 0) {} + MinDistCalc( double dD, double dP, double dPm, double dPM, Point3d pT) + : dDist( dD), dPar( dP), dParMin( dPm), dParMax( dPM), + bParMinSing( false), bParMaxSing( false), ptQ( pT) {} + MinDistCalc( double dD, double dP, double dPm, double dPM, bool bMinS, bool bMaxS, Point3d pT) + : dDist( dD), dPar( dP), dParMin( dPm), dParMax( dPM), + bParMinSing( bMinS), bParMaxSing( bMaxS), ptQ( pT) {} +} ; +typedef std::vector MDCVECTOR ; // vettore di MinDistCalc + +//---------------------------------------------------------------------------- +enum MdiType { MDI_NORMAL, MDI_START_CONT, MDI_END_CONT} ; +struct MinDistInfo { + MdiType nFlag ; + double dPar ; + Point3d ptQ ; + MinDistInfo( void) + : nFlag( MDI_NORMAL), dPar( 0), ptQ( 0, 0, 0) {} + MinDistInfo( MdiType nF, double dP, Point3d pT) + : nFlag( nF), dPar( dP), ptQ( pT) {} +} ; +typedef std::vector MDIVECTOR ; // vettore di MinDistInfo + +//---------------------------------------------------------------------------- +bool +PolishMinDistPointCurve( const Point3d& ptP, const ICurve& cCurve, + const MinDistCalc& approxMin, double& dPrevPar, Point3d& ptQ) { - const double LIN_TOL_APPROX = 1 ; - const double ANG_TOL_APPROX_DEG = 45 ; - - - // distanza non calcolata - m_dDist = - 1 ; - - if ( ! CrvBez.IsValid()) - return ; - - // creo una polilinea di approssimazione - PolyLine PL ; - if ( ! CrvBez.ApproxWithLines( LIN_TOL_APPROX, ANG_TOL_APPROX_DEG, PL)) - return ; - - // cerco la minima distanza per la polilinea - double dSqDist ; - double dMinPar ; - double dMinParIni ; - double dMinParFin ; + const int MAX_COUNT = 16 ; + ICurve::Side nSide ; double dPar ; - double dUIni ; - double dUFin ; - Point3d ptMinDist ; - Point3d ptIni ; - Point3d ptFin ; - bool bFound = false ; - bool bOnEnd = false ; - double dSqMinDist = INFINITO * INFINITO ; - for ( bool bLine = PL.GetFirstULine( &dUIni, &ptIni, &dUFin, &ptFin) ; - bLine ; - bLine = PL.GetNextULine( &dUIni, &ptIni, &dUFin, &ptFin)) { - // calcolo la distanza del punto dal segmento - DistPointLine dstPtLn( ptP, ptIni, ptFin) ; - if ( dstPtLn.GetSqDist( dSqDist) && dSqDist < dSqMinDist) { - bFound = true ; - dSqMinDist = dSqDist ; - dstPtLn.GetPointMinDist( ptMinDist) ; - dstPtLn.GetParamAtPointMinDist( dPar) ; - dMinPar = ( 1 - dPar) * dUIni + dPar * dUFin ; - dMinParIni = dUIni ; - dMinParFin = dUFin ; - bOnEnd = ( AreSamePointNear( ptMinDist, ptFin)) ; - } - else if ( bOnEnd) { - bOnEnd = false ; - dMinParFin = dUFin ; - } - } - if ( ! bFound) - return ; - - // raffino i punti trovati (con algoritmo tipo Newton da TheNurbsBook pag 230) - const int MAX_COUNT = 10 ; - double dMinDist = sqrt( dSqMinDist) ; - double dPrevPar ; double dTemp ; double dSqCosA ; - Point3d ptQ ; Vector3d vtDer1 ; Vector3d vtDer2 ; Vector3d vtDiff ; + + + // raffino i punti trovati (con algoritmo tipo Newton da TheNurbsBook pag 230) int nCount = 0 ; - dPar = dMinPar ; + bool bClampedFromSing = false ; + dPar = approxMin.dPar ; do { // contatore iterazioni nCount ++ ; // calcolo P, D1 e D2 - CrvBez.GetPointD1D2( dPar, ICurve::FROM_MINUS, ptQ, &vtDer1, &vtDer2) ; - // se D1 nulla e distanza sicuramente non minima, riprovo sul centro dell'intervallo - if ( vtDer1.IsZero() && Dist( ptQ, ptP) > __min( 1.2 * dMinDist, dMinDist + LIN_TOL_APPROX)) { - dPar = 0.5 * ( dMinParIni + dMinParFin) ; - CrvBez.GetPointD1D2( dPar, ICurve::FROM_MINUS, ptQ, &vtDer1, &vtDer2) ; - } + nSide = ( fabs( dPar - approxMin.dParMin) < EPS_ZERO) ? ICurve::FROM_PLUS : ICurve::FROM_MINUS ; + cCurve.GetPointD1D2( dPar, nSide, ptQ, &vtDer1, &vtDer2) ; // vettore dal punto al piede sulla curva vtDiff = ptQ - ptP ; // angolo tra vettore e tangente @@ -108,24 +93,204 @@ DistPointCrvBezier::DistPointCrvBezier( const Point3d& ptP, const ICurveBezier& dTemp = vtDer2 * vtDiff + vtDer1.SqLen() ; if ( fabs( dTemp) > EPS_ZERO) dPar = dPrevPar - ( vtDer1 * vtDiff) / dTemp ; - // clipping parametro in 0-1 - if ( dPar < dMinParIni) - dPar = dMinParIni ; - else if ( dPar > dMinParFin) - dPar = dMinParFin ; + // clipping parametro + if ( dPar < approxMin.dParMin) { + if ( approxMin.bParMinSing && ! bClampedFromSing) { + dPar = approxMin.dParMax ; + bClampedFromSing = true ; + } + else + dPar = approxMin.dParMin ; + } + else if ( dPar > approxMin.dParMax) { + if ( approxMin.bParMaxSing && ! bClampedFromSing) { + dPar = approxMin.dParMin ; + bClampedFromSing = true ; + } + else + dPar = approxMin.dParMax ; + } } while ( nCount < MAX_COUNT && fabs( dPar - dPrevPar) > EPS_ZERO && fabs( dSqCosA) > COS_ORTO_ANG_ZERO * COS_ORTO_ANG_ZERO) ; -#if defined( _DEBUG) if ( nCount == MAX_COUNT) - EGT_TRACE( "Exceeded MAX_COUNT in PolishingDistPointCrv.") ; -#endif + LOG_DBG_ERR( GetEGkLogger(), "ERROR : Exceeded recursions") ; + + return true ; +} + + +//---------------------------------------------------------------------------- +DistPointCrvBezier::DistPointCrvBezier( const Point3d& ptP, const ICurveBezier& CrvBez) +{ + // distanza non calcolata + m_dDist = - 1 ; + + if ( ! CrvBez.IsValid()) + return ; + + // creo una polilinea di approssimazione + PolyLine PL ; + if ( ! CrvBez.ApproxWithLines( LIN_TOL_APPROX, ANG_TOL_APPROX_DEG, PL)) + return ; + + // cerco la minima distanza per la polilinea + MinDistCalc approxMin ; + MDCVECTOR vApproxMin ; + MDCVECTOR::iterator Iter ; + vApproxMin.reserve( 4) ; + double dSqDist ; + double dPar ; + double dUIni ; + double dUFin ; + Point3d ptIni ; + Point3d ptFin ; + bool bFound = false ; + bool bOnEnd = false ; + double dMinDist ; + double dSqMinDist ; + for ( bool bLine = PL.GetFirstULine( &dUIni, &ptIni, &dUFin, &ptFin) ; + bLine ; + bLine = PL.GetNextULine( &dUIni, &ptIni, &dUFin, &ptFin)) { + // calcolo la distanza del punto dal segmento + DistPointLine dstPtLn( ptP, ptIni, ptFin) ; + if ( ! dstPtLn.GetSqDist( dSqDist)) + continue ; + // altro punto con la stessa minima distanza già trovata + if ( bFound && fabs( dSqDist - dSqMinDist) < 2 * dMinDist * LIN_TOL_APPROX) { + // salvo i dati nella struttura + approxMin.dDist = dMinDist ; + dstPtLn.GetPointMinDist( approxMin.ptQ) ; + dstPtLn.GetParamAtPointMinDist( dPar) ; + approxMin.dPar = ( 1 - dPar) * dUIni + dPar * dUFin ; + approxMin.dParMin = dUIni ; + approxMin.dParMax = dUFin ; + // setto che il punto è alla fine + bOnEnd = ( AreSamePointNear( approxMin.ptQ, ptFin)) ; + // aggiungo alla lista + vApproxMin.push_back( approxMin) ; + } + // primo punto o punto con minima distanza più bassa + else if ( ! bFound || dSqDist < dSqMinDist) { + // aggiorno i minimi + bFound = true ; + dSqMinDist = dSqDist ; + dMinDist = sqrt( dSqMinDist) ; + // salvo i dati nella struttura + approxMin.dDist = dMinDist ; + dstPtLn.GetPointMinDist( approxMin.ptQ) ; + dstPtLn.GetParamAtPointMinDist( dPar) ; + approxMin.dPar = ( 1 - dPar) * dUIni + dPar * dUFin ; + approxMin.dParMin = dUIni ; + approxMin.dParMax = dUFin ; + // setto che il punto è alla fine + bOnEnd = ( AreSamePointNear( approxMin.ptQ, ptFin)) ; + // il nuovo vettore deve contenere solo quest'ultimo minimo + vApproxMin.clear() ; + vApproxMin.push_back( approxMin) ; + } + // il minimo era alla fine, devo allargare l'intervallo di raffinamento del parametro + else if ( bOnEnd) { + bOnEnd = false ; + vApproxMin.back().dParMax = dUFin ; + } + } + if ( ! bFound) + return ; + // verifico presenza singolarità agli estremi degli intervalli trovati + double dSingP ; + if ( CrvBez.GetSingularParam( dSingP) == 0) + dSingP = - 1 ; + for ( Iter = vApproxMin.begin() ; Iter != vApproxMin.end() ; ++Iter) { + // imposto flag per singolarità agli estremi + (*Iter).bParMinSing = fabs( (*Iter).dParMin - dSingP) < EPS_SMALL ; + (*Iter).bParMaxSing = fabs( (*Iter).dParMax - dSingP) < EPS_SMALL ; + } + + // raffino i punti trovati + double dPolishedPar ; + Point3d ptPolishedQ ; + for ( Iter = vApproxMin.begin() ; Iter != vApproxMin.end() ; ++Iter) { + // eseguo raffinamento + if ( PolishMinDistPointCurve( ptP, CrvBez, *Iter, dPolishedPar, ptPolishedQ)) { + (*Iter).dDist = Dist( ptP, ptPolishedQ) ; + (*Iter).dPar = dPolishedPar ; + (*Iter).ptQ = ptPolishedQ ; + } + else + (*Iter).dDist = INFINITO ; + } + + // determino i minimi raffinati da tenere + MDIVECTOR vQmin ; + bFound = false ; + for ( Iter = vApproxMin.begin() ; Iter != vApproxMin.end() ; ++Iter) { + // altro punto con la stessa minima distanza + if ( bFound && fabs( (*Iter).dDist - dMinDist) < EPS_SMALL) { + // se abbastanza lontano lo aggiungo + if ( SqDist( (*Iter).ptQ, vQmin.back().ptQ) > 1) + vQmin.push_back( MinDistInfo( MDI_NORMAL, (*Iter).dPar, (*Iter).ptQ)) ; + // altrimenti lo sostituisco se distanza minore + else if ( (*Iter).dDist < dMinDist) + vQmin.back() = MinDistInfo( MDI_NORMAL, (*Iter).dPar, (*Iter).ptQ) ; + } + // primo punto o punto con minima distanza più bassa + else if ( ! bFound || (*Iter).dDist < dMinDist) { + // aggiorno i minimi + bFound = true ; + dMinDist = (*Iter).dDist ; + // il nuovo vettore deve contenere solo quest'ultimo minimo + vQmin.clear() ; + vQmin.push_back( MinDistInfo( MDI_NORMAL, (*Iter).dPar, (*Iter).ptQ)) ; + } + } + if ( vQmin.empty()) + return ; + + // se 2 o più minimi, verifico se tratto continuo + if ( vQmin.size() >= 2) { + bool bCont = true ; + double dU ; + Point3d ptQ ; + // se tutti i punti intermedi hanno la stessa distanza, è una zona continua + for ( int i = 1 ; i < (int) vQmin.size() ; ++ i) { + dU = 0.5 * ( vQmin[i-1].dPar + vQmin[i].dPar) ; + CrvBez.GetPointD1D2( dU, ICurve::FROM_MINUS, ptQ) ; + if ( fabs( SqDist( ptP, ptQ) - dMinDist * dMinDist) > 2 * dMinDist * EPS_SMALL) { + bCont = false ; + break ; + } + } + // se zona continua è un arco di circonferenza, tengo solo primo e ultimo punto e imposto opportuni flag + if ( bCont) { + // se è praticamente tutta la curva, la faccio diventare tutta + if ( ( vQmin.back().dPar - vQmin.front().dPar) > 0.8) { + // primo elemento == inizio della curva + vQmin[0].dPar = 0 ; + vQmin[0].nFlag = MDI_START_CONT ; + CrvBez.GetStartPoint( vQmin[0].ptQ) ; + // ultimo elemento == fine curva + vQmin[1].dPar = 1 ; + vQmin[1].nFlag = MDI_END_CONT ; + CrvBez.GetEndPoint( vQmin[1].ptQ) ; + } + else { + // sistemo flag primo elemento + vQmin[0].nFlag = MDI_START_CONT ; + // sposto ultimo elemento al secondo posto e ne sistemo il flag + vQmin[1] = vQmin.back() ; + vQmin[1].nFlag = MDI_END_CONT ; + } + // cancello tutti gli altri elementi + vQmin.erase( vQmin.begin() + 2, vQmin.end()) ; + } + } // assegno i dati - m_bDet = true ; - m_dDist = Dist( ptP, ptQ) ; - m_dParam = dPrevPar ; - m_ptMinDist = ptQ ; + m_bDet = ( vQmin[0].nFlag == MDI_NORMAL) ; + m_dDist = dMinDist ; + m_dParam = vQmin[0].dPar ; + m_ptMinDist = vQmin[0].ptQ ; } //---------------------------------------------------------------------------- diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index da8bb591d14e9e73b455f88f7c19ea173ebad923..e01ef677d694ff4db788f1cc4058f704e476e1eb 100644 GIT binary patch delta 85 zcmexl{m6R5H#SD2&G-4vGfl1&(wW@A$uaq#pd6#oWJO`!&0a!ojFWq~A52yf7Mt8B e@@jGqYry6syi$x{Ew;>XWok?yC7b(1QaAxd`y6%v delta 85 zcmaE4{mFX6H#SDY&EMG+nHdcy9~9Kx?8c?W#IDC+zyQQR0rAOoBCj^*@MbXrm0jnv UMJRJ)oZKP&09ongI*}Ak0H5s{=l}o! diff --git a/EgtGeomKernel.vcxproj b/EgtGeomKernel.vcxproj index e7a722b..3896a0f 100644 --- a/EgtGeomKernel.vcxproj +++ b/EgtGeomKernel.vcxproj @@ -59,6 +59,7 @@ WIN32;_WINDOWS;I_AM_EGK;_DEBUG;_USRDLL;%(PreprocessorDefinitions) true CompileAsCpp + false Windows @@ -95,7 +96,7 @@ copy $(TargetPath) \EgtProg\Dll Windows - true + false true true @@ -112,8 +113,7 @@ copy $(TargetPath) \EgtProg\Dll $(IntDir);%(AdditionalIncludeDirectories) - copy $(TargetDir)$(TargetName).pdb \EgtDev\Lib\ -copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\ + copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\ copy $(TargetPath) \EgtProg\Dll @@ -142,6 +142,7 @@ copy $(TargetPath) \EgtProg\Dll + Create Create @@ -176,7 +177,9 @@ copy $(TargetPath) \EgtProg\Dll + + @@ -198,6 +201,7 @@ copy $(TargetPath) \EgtProg\Dll + diff --git a/EgtGeomKernel.vcxproj.filters b/EgtGeomKernel.vcxproj.filters index ccf4fb4..96e4f91 100644 --- a/EgtGeomKernel.vcxproj.filters +++ b/EgtGeomKernel.vcxproj.filters @@ -93,9 +93,6 @@ File di origine\Geo - - File di origine\Gdb - File di origine\Base @@ -111,6 +108,12 @@ File di origine\Distanze + + File di origine\Distanze + + + File di origine\Base + @@ -266,6 +269,15 @@ File di intestazione + + File di intestazione + + + File di intestazione + + + File di intestazione + diff --git a/GdbExecutor.cpp b/GdbExecutor.cpp index e7f2ae0..04fb09c 100644 --- a/GdbExecutor.cpp +++ b/GdbExecutor.cpp @@ -45,7 +45,9 @@ CreateGdbExecutor( void) GdbExecutor::GdbExecutor( void) { m_pGDB = nullptr ; - // alias predefiniti + // dimensioni mappa dei nomi + m_NameMap.rehash( 100) ; + // nomi predefiniti m_NameMap.insert( pair< string, int>( "$ROOT", GDB_ID_ROOT)) ; m_NameMap.insert( pair< string, int>( "$NN", ID_NO)) ; } @@ -73,13 +75,16 @@ GdbExecutor::Execute( const string& sCmd1, const string& sCmd2, const STRVECTOR& // output di debug - sOut = "Cmd = [" + sCmd1 + "][" + sCmd2 + "] Par = [" ; + sOut = " " + sCmd1 ; + if ( ! sCmd2.empty()) + sOut += "." + sCmd2 ; + sOut += "( " ; for ( theConstIter = vsParams.begin() ; theConstIter != vsParams.end() ; ++theConstIter) { if ( theConstIter != vsParams.begin()) - sOut += "][" ; + sOut += ", " ; sOut += *theConstIter ; } - sOut += "]" ; + sOut += ")" ; LOG_DBG_INFO( GetEGkLogger(), sOut.c_str()) // esecuzione comando @@ -1061,13 +1066,29 @@ GdbExecutor::ExecuteOutScl( const string& sCmd2, const STRVECTOR& vsParams) } // emetto gruppo else if ( sCmd2 == "PUTGR") { + int nFlag ; int nId ; - // un parametro - if ( vsParams.size() != 1) + STRVECTOR vsNames ; + STRVECTOR::iterator Iter ; + // almeno un parametro + if ( vsParams.size() < 1) return false ; - nId = GetIdParam( vsParams[0]) ; - // emetto gruppo e suoi sottoposti - return OutGroupScl( nId) ; + // se esiste il secondo + if ( vsParams.size() == 2) + FromString( vsParams[1], nFlag) ; + else + nFlag = 0 ; + // recupero lista nomi + if ( ! GetNamesParam( vsParams[0], vsNames)) + return false ; + // esecuzione + for ( Iter = vsNames.begin() ; Iter != vsNames.end() ; ++Iter) { + // recupero l'oggetto ed eseguo l'output + nId = GetIdParam( *Iter) ; + if ( ! OutGroupScl( nId, nFlag)) + return false ; + } + return true ; } // emetto oggetto geometrico else if ( sCmd2 == "PUT") { @@ -1101,7 +1122,7 @@ GdbExecutor::ExecuteOutScl( const string& sCmd2, const STRVECTOR& vsParams) //---------------------------------------------------------------------------- bool -GdbExecutor::OutGroupScl( int nId) +GdbExecutor::OutGroupScl( int nId, int nFlag) { bool bNext ; int nParentId ; @@ -1126,11 +1147,11 @@ GdbExecutor::OutGroupScl( int nId) while ( bNext) { nGdbType = Iter.GetGdbType() ; if ( nGdbType == GDB_GEO) { - if ( ! m_OutScl.PutCurve( Iter.GetGeoObj(), 0)) + if ( ! m_OutScl.PutCurve( Iter.GetGeoObj(), nFlag)) return false ; } else if ( nGdbType == GDB_GROUP) { - if ( ! OutGroupScl( Iter.GetId())) + if ( ! OutGroupScl( Iter.GetId(), nFlag)) return false ; } bNext = Iter.GoToNext() ; diff --git a/GdbExecutor.h b/GdbExecutor.h index 05a0407..4c0d134 100644 --- a/GdbExecutor.h +++ b/GdbExecutor.h @@ -54,7 +54,7 @@ class GdbExecutor : public IGdbExecutor bool ExecuteLoad( const STRVECTOR& vsParams) ; bool ExecuteSave( const STRVECTOR& vsParams) ; bool ExecuteOutScl( const std::string& sCmd2, const STRVECTOR& vsParams) ; - bool OutGroupScl( int nId) ; + bool OutGroupScl( int nId, int nFlag) ; private : IGeomDB* m_pGDB ; diff --git a/OutScl.cpp b/OutScl.cpp index 1139e49..13a664f 100644 --- a/OutScl.cpp +++ b/OutScl.cpp @@ -305,6 +305,34 @@ OutScl::ArcCurvOrTgOrNone( const CrvPointDiffGeom& oDiffG) return true ; } +//---------------------------------------------------------------------------- +bool +OutScl::NormalOrNone( const CrvPointDiffGeom& oDiffG) +{ + const double NORM_LEN = 10 ; + + + // normale + if ( oDiffG.nStatus == CrvPointDiffGeom::NCRV && fabs( oDiffG.dCurv) > EPS_ZERO) { + double dLen = __min( NORM_LEN, 1 / oDiffG.dCurv) ; + Line2P( oDiffG.ptP - NORM_LEN * oDiffG.vtN, oDiffG.ptP + dLen * oDiffG.vtN) ; + } + // segmento perpendicolare alla tangente + else if ( oDiffG.nStatus == CrvPointDiffGeom::TANG && ! oDiffG.vtT.IsSmall()) { + Vector3d vtN = oDiffG.vtT ^ Z_AX ; + if ( ! vtN.Normalize()) { + vtN = oDiffG.vtT ^ Y_AX ; + vtN.Normalize() ; + } + Line2P( oDiffG.ptP - NORM_LEN * vtN, oDiffG.ptP + NORM_LEN * vtN) ; + } + // altrimenti cerchietto + else + CircleCR( oDiffG.ptP, 1) ; + + return true ; +} + //---------------------------------------------------------------------------- bool OutScl::PutCurve( const IGeoObj* pCurve, int nFlag) @@ -314,7 +342,7 @@ OutScl::PutCurve( const IGeoObj* pCurve, int nFlag) switch ( pCurve->GetType()) { case CRV_LINE : - return PutCurveLine( *GetCurveLine( pCurve)) ; + return PutCurveLine( *GetCurveLine( pCurve), nFlag) ; case CRV_ARC : return PutCurveArc( *GetCurveArc( pCurve), nFlag) ; case CRV_BEZ : @@ -328,12 +356,42 @@ OutScl::PutCurve( const IGeoObj* pCurve, int nFlag) //---------------------------------------------------------------------------- bool -OutScl::PutCurveLine( const ICurveLine& CrvLine) +OutScl::PutCurveLine( const ICurveLine& CrvLine, int nFlag) { - // scrittura comandi SCL + bool bFound ; + double dU ; + PolyLine PL ; + CrvPointDiffGeom oDiffG ; + + + // scrittura della linea Remark( "CurveLine") ; Line2P( CrvLine.GetStart(), CrvLine.GetEnd()) ; + // se richieste tangenti e curvature + if ( ( nFlag & 1) != 0) { + // ciclo per disegnare le derivate + Remark( "LineTangents+Der2") ; + for ( bFound = PL.GetFirstU( dU) ; bFound ; bFound = PL.GetNextU( dU)) { + // ricavo il punto, la tangente, la normale e la curvatura + CrvLine.GetPointDiffGeom( dU, ICurve::FROM_MINUS, oDiffG) ; + // curvatura o tangente o niente + ArcCurvOrTgOrNone( oDiffG) ; + } + } + + // se richieste normali + if ( ( nFlag & 2) != 0) { + // ciclo per disegnare le normali + Remark( "LineNormals") ; + for ( bFound = PL.GetFirstU( dU) ; bFound ; bFound = PL.GetNextU( dU)) { + // ricavo il punto, la tangente, la normale e la curvatura + CrvLine.GetPointDiffGeom( dU, ICurve::FROM_MINUS, oDiffG) ; + // curvatura o tangente o niente + NormalOrNone( oDiffG) ; + } + } + return true ; } @@ -355,7 +413,8 @@ OutScl::PutCurveArc( const ICurveArc& CrvArc, int nFlag) for ( bFound = PL.GetFirstLine( ptIni, ptFin) ; bFound ; bFound = PL.GetNextLine( ptIni, ptFin)) { Line2P( ptIni, ptFin) ; } - // se richieste derivate e curvature + + // se richieste tangenti e curvature if ( ( nFlag & 1) != 0) { // ciclo per disegnare le derivate Remark( "ArcTangents+Der2") ; @@ -367,6 +426,18 @@ OutScl::PutCurveArc( const ICurveArc& CrvArc, int nFlag) } } + // se richieste normali + if ( ( nFlag & 2) != 0) { + // ciclo per disegnare le normali + Remark( "ArcNormals") ; + for ( bFound = PL.GetFirstU( dU) ; bFound ; bFound = PL.GetNextU( dU)) { + // ricavo il punto, la tangente, la normale e la curvatura + CrvArc.GetPointDiffGeom( dU, ICurve::FROM_MINUS, oDiffG) ; + // curvatura o tangente o niente + NormalOrNone( oDiffG) ; + } + } + return true ; } @@ -383,7 +454,7 @@ OutScl::PutCurveBez( const ICurveBezier& CrvBez, int nFlag) // se richiesto anche il poligono di controllo - if ( ( nFlag & 2) != 0) + if ( ( nFlag & 4) != 0) PutPolygBez( CrvBez) ; // ciclo per disegnare i segmenti @@ -392,7 +463,8 @@ OutScl::PutCurveBez( const ICurveBezier& CrvBez, int nFlag) for ( bFound = PL.GetFirstLine( ptIni, ptFin) ; bFound ; bFound = PL.GetNextLine( ptIni, ptFin)) { Line2P( ptIni, ptFin) ; } - // se richieste derivate e curvature + + // se richieste tangenti e curvature if ( ( nFlag & 1) != 0) { // ciclo per disegnare le derivate e le curvature Remark( "BezierTangents+Der2") ; @@ -404,6 +476,18 @@ OutScl::PutCurveBez( const ICurveBezier& CrvBez, int nFlag) } } + // se richieste normali + if ( ( nFlag & 2) != 0) { + // ciclo per disegnare le derivate e le curvature + Remark( "BezierNormals") ; + for ( bFound = PL.GetFirstU( dU) ; bFound ; bFound = PL.GetNextU( dU)) { + // ricavo il punto, la tangente, la normale e la curvatura + CrvBez.GetPointDiffGeom( dU, ICurve::FROM_MINUS, oDiffG) ; + // curvatura o tangente o niente + NormalOrNone( oDiffG) ; + } + } + return true ; } diff --git a/OutScl.h b/OutScl.h index 5f2b85e..e380823 100644 --- a/OutScl.h +++ b/OutScl.h @@ -35,7 +35,7 @@ class OutScl bool SetPartLay( std::string sPart, std::string sLay) ; bool SetPartLayRef( std::string sPart, std::string sLay, const Frame3d& frFrame) ; bool PutCurve( const IGeoObj* pCurve, int nFlag) ; - bool PutCurveLine( const ICurveLine& CrvLine) ; + bool PutCurveLine( const ICurveLine& CrvLine, int nFlag) ; bool PutCurveArc( const ICurveArc& CrvArc, int nFlag) ; bool PutCurveBez( const ICurveBezier& CrvBez, int nFlag) ; bool PutPolygBez( const ICurveBezier& CrvBez) ; @@ -51,6 +51,7 @@ class OutScl bool ArcCPA( const Point3d& ptCen, const Point3d& ptMed, double dAngCenDeg) ; bool CircleCR( const Point3d& ptCen, double dRad) ; bool ArcCurvOrTgOrNone( const CrvPointDiffGeom& oDiffG) ; + bool NormalOrNone( const CrvPointDiffGeom& oDiffG) ; private : std::ofstream m_ofFile ; diff --git a/PolynomialPoint3d.cpp b/PolynomialPoint3d.cpp new file mode 100644 index 0000000..4ba4c1d --- /dev/null +++ b/PolynomialPoint3d.cpp @@ -0,0 +1,225 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2014 +//---------------------------------------------------------------------------- +// File : PolynomialPoint3d.cpp Data : 12.01.14 Versione : 1.5a2 +// Contenuto : Implementazione classe polinomio con coefficienti Point3d. +// +// +// +// Modifiche : 12.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +//--------------------------- Include ---------------------------------------- +#include "stdafx.h" +#include "PolynomialPoint3d.h" + + +//---------------------------------------------------------------------------- +bool +PolynomialPoint3d::SetDegree( int nDegree) +{ + // gradi negativi non hanno senso + if ( nDegree < 0) + return false ; + // pulisco, alloco e inizializzo a 0 + try { + m_Coeff.clear() ; + m_Coeff.reserve( nDegree + 1) ; + m_nDegree = nDegree ; + for ( int i = 0 ; i <= m_nDegree ; ++ i) + m_Coeff.push_back( Point3d( 0, 0, 0)) ; + return true ; + } + catch (...) { + return false ; + } +} + +//---------------------------------------------------------------------------- +bool +PolynomialPoint3d::EnsureDegree( int nDegree) +{ + // se il grado è già adeguato non devo fare alcunché + if ( nDegree <= m_nDegree) + return true ; + // alloco e inizializzo a 0 i nuovi coefficienti + try { + m_Coeff.reserve( nDegree + 1) ; + for ( int i = m_nDegree + 1 ; i <= nDegree ; ++ i) + m_Coeff.push_back( Point3d()) ; + m_nDegree = nDegree ; + return true ; + } + catch (...) { + return false ; + } +} + +//---------------------------------------------------------------------------- +bool +PolynomialPoint3d::SetCoeff( int nPower, Point3d& ptP) +{ + if ( nPower < 0 || nPower > m_nDegree) + return false ; + + m_Coeff[nPower] = ptP ; + return true ; +} + +//---------------------------------------------------------------------------- +bool +PolynomialPoint3d::Set( int nDegree, const PNTVECTOR& vP) +{ + if ( ! SetDegree( nDegree)) + return false ; + + for ( int i = 0 ; i <= nDegree ; i ++) + m_Coeff[i] = vP[i] ; + return true ; +} + +//---------------------------------------------------------------------------- +bool +PolynomialPoint3d::SetToConstant( const Point3d& ptP) +{ + if ( ! SetDegree( 0)) + return false ; + m_Coeff[0] = ptP ; + return true ; +} + +//---------------------------------------------------------------------------- +const PolynomialPoint3d& +PolynomialPoint3d::operator +=( const PolynomialPoint3d& pol3P) +{ + // mi assicuro che il polinomio risultante abbia grado sufficiente + EnsureDegree( pol3P.GetDegree()) ; + + // eseguo la somma + for ( int i = 0 ; i <= pol3P.GetDegree() ; ++ i) + m_Coeff[i] += pol3P.m_Coeff[i] ; + + return *this ; +} + +//---------------------------------------------------------------------------- +const PolynomialPoint3d& +PolynomialPoint3d::operator -=( const PolynomialPoint3d& pol3P) +{ + // mi assicuro che il polinomio risultante abbia grado sufficiente + EnsureDegree( pol3P.GetDegree()) ; + + // eseguo la somma + for ( int i = 0 ; i <= pol3P.GetDegree() ; ++ i) + m_Coeff[i] += ( - pol3P.m_Coeff[i]) ; + + return *this ; +} + +//---------------------------------------------------------------------------- +const PolynomialPoint3d& +PolynomialPoint3d::operator *=( const Polynomial& polP) +{ + // copio il polinomio corrente + PolynomialPoint3d pol3C = *this ; + + // pulisco e imposto il grado del polinomio risultante + SetDegree( pol3C.GetDegree() + polP.GetDegree()) ; + + // eseguo il prodotto + for ( int i = 0 ; i <= pol3C.GetDegree() ; ++ i) { + for ( int j = 0 ; j <= polP.GetDegree() ; ++ j) + m_Coeff[i+j] += pol3C.m_Coeff[i] * polP.GetCoeff(j) ; + } + + return *this ; +} + +//---------------------------------------------------------------------------- +void +PolynomialPoint3d::Derive( void) +{ + // polinomio non inizializzato + if ( m_nDegree < 0) + return ; + // polinomio costante + if ( m_nDegree == 0) { + m_Coeff[0] = Point3d( 0, 0, 0) ; + return ; + } + // caso normale + for ( int i = 0 ; i < m_nDegree ; ++ i) + m_Coeff[i] = ( i + 1) * m_Coeff[i+1] ; + m_Coeff.pop_back() ; + -- m_nDegree ; +} + +//---------------------------------------------------------------------------- +void +PolynomialPoint3d::Derive( const PolynomialPoint3d& pol3P) +{ + operator=( pol3P) ; + Derive() ; +} + +//---------------------------------------------------------------------------- +void +PolynomialPoint3d::AdjustDegree( void) +{ + // se il coefficiente del grado più alto è zero, diminuisco il grado + while ( m_nDegree >= 0 && + fabs( m_Coeff[m_nDegree].x) < DBL_EPSILON && + fabs( m_Coeff[m_nDegree].y) < DBL_EPSILON && + fabs( m_Coeff[m_nDegree].z) < DBL_EPSILON) { + m_Coeff.pop_back() ; + -- m_nDegree ; + } +} + +//---------------------------------------------------------------------------- +Point3d +PolynomialPoint3d::Evaluate( double dVal) +{ + // polinomio non inizializzato + if ( m_nDegree < 0) + return Point3d( 0, 0, 0) ; + // caso normale + Point3d ptRes = m_Coeff[m_nDegree] ; + for ( int i = m_nDegree - 1 ; i >= 0 ; -- i) + ptRes = ptRes * dVal + m_Coeff[i] ; + + return ptRes ; +} + +//---------------------------------------------------------------------------- +int +PolynomialPoint3d::FindMainComponentRoots( DBLVECTOR& vdRoot) +{ + // cerco la componente più significativa tra x, y e z + Point3d ptSumm ; + for ( int i = 0 ; i <= m_nDegree ; ++ i) { + ptSumm.x += fabs( m_Coeff[i].x) ; + ptSumm.y += fabs( m_Coeff[i].y) ; + ptSumm.z += fabs( m_Coeff[i].z) ; + } + // la copio in un polinomio numerico + Polynomial polP ; + if ( ! polP.SetDegree( m_nDegree)) + return 0 ; + if ( ptSumm.x > ptSumm.y && ptSumm.x > ptSumm.z) { + for ( int i = 0 ; i <= m_nDegree ; ++ i) + polP.SetCoeff( i, m_Coeff[i].x) ; + } + else if ( ptSumm.y > ptSumm.z) { + for ( int i = 0 ; i <= m_nDegree ; ++ i) + polP.SetCoeff( i, m_Coeff[i].y) ; + } + else { + for ( int i = 0 ; i <= m_nDegree ; ++ i) + polP.SetCoeff( i, m_Coeff[i].z) ; + } + // calcolo le radici reali di questo polinomio numerico + return polP.FindRoots( vdRoot) ; +} \ No newline at end of file diff --git a/PolynomialPoint3d.h b/PolynomialPoint3d.h new file mode 100644 index 0000000..d6c0fce --- /dev/null +++ b/PolynomialPoint3d.h @@ -0,0 +1,193 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2014 +//---------------------------------------------------------------------------- +// File : PolynomialPoint3d.h Data : 11.01.14 Versione : 1.5a2 +// Contenuto : Dichiarazione classe polinomio con coefficienti Point3d. +// +// +// +// Modifiche : 11.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +#pragma once + +#include "/EgtDev/Include/ENkPolynomial.h" +#include "/EgtDev/Include/EGkGeoCollection.h" + + +//---------------------------------------------------------------------------- +class PolynomialPoint3d +{ + public : + PolynomialPoint3d( void) : m_nDegree( -1) {} + bool SetDegree( int nDegree) ; + bool SetCoeff( int nPower, Point3d& ptP) ; + bool Set( int nDegree, const PNTVECTOR& vP) ; + bool SetToConstant( const Point3d& ptP) ; + const PolynomialPoint3d& operator =( const PolynomialPoint3d& pol3S) + { if ( &pol3S != this) { + SetDegree( pol3S.m_nDegree) ; + for ( int i = 0 ; i <= m_nDegree ; ++ i) + m_Coeff[i] = pol3S.m_Coeff[i] ;} + return *this ; } + + public : + int GetDegree( void) const { return m_nDegree ; } + Point3d GetCoeff( int nPower) const + { if ( nPower < 0 || nPower > m_nDegree) + return Point3d( 0, 0, 0) ; + return m_Coeff[nPower] ; } + const PolynomialPoint3d& operator +=( const PolynomialPoint3d& pol3P) ; + const PolynomialPoint3d& operator -=( const PolynomialPoint3d& pol3P) ; + const PolynomialPoint3d& operator *=( const Polynomial& polP) ; + void Derive( void) ; + void Derive( const PolynomialPoint3d& pol3P) ; + void AdjustDegree( void) ; + Point3d Evaluate( double dVal) ; + int FindMainComponentRoots( DBLVECTOR& vdRoot) ; + + private : + bool EnsureDegree( int nDegree) ; + + private : + int m_nDegree ; + PNTVECTOR m_Coeff ; +} ; + + +//---------------------------------------------------------------------------- +// Somma +//---------------------------------------------------------------------------- +inline PolynomialPoint3d +operator+( const PolynomialPoint3d& pol3P1, const PolynomialPoint3d& pol3P2) +{ + PolynomialPoint3d pol3Summ = pol3P1 ; + pol3Summ += pol3P2 ; + return pol3Summ ; +} + +//---------------------------------------------------------------------------- +// Differenza +//---------------------------------------------------------------------------- +inline PolynomialPoint3d +operator-( const PolynomialPoint3d& pol3P1, const PolynomialPoint3d& pol3P2) +{ + PolynomialPoint3d pol3Diff = pol3P1 ; + pol3Diff -= pol3P2 ; + return pol3Diff ; +} + +//---------------------------------------------------------------------------- +// Moltiplicazione con un polinomio di numeri +//---------------------------------------------------------------------------- +inline PolynomialPoint3d +operator*( const PolynomialPoint3d& pol3P1, const Polynomial& polP2) +{ + PolynomialPoint3d pol3Mul = pol3P1 ; + pol3Mul *= polP2 ; + return pol3Mul ; +} + + +#if 0 +//---------------------------------------------------------------------------- +void +PolynomialSumm( PNTVECTOR& vSou1, PNTVECTOR& vSou2, PNTVECTOR& vSumm) +{ + int nDeg1 = vSou1.size() - 1 ; + int nDeg2 = vSou2.size() - 1 ; + int nMin = (( nDeg1 <= nDeg2) ? nDeg1 : nDeg2) ; + int nMax = (( nDeg1 >= nDeg2) ? nDeg1 : nDeg2) ; + vSumm.clear() ; + vSumm.reserve( nMax + 1) ; + for ( int i = 0 ; i < nMin ; ++ i) + vSumm.push_back( vSou1[i] + vSou2[i]) ; + if ( nDeg1 > nDeg2) { + for ( int i = nMin ; i < nDeg1 ; ++ i) + vSumm.push_back( vSou1[i]) ; + } + else if ( nDeg1 < nDeg2) { + for ( int i = nMin ; i < nDeg2 ; ++ i) + vSumm.push_back( vSou2[i]) ; + } +} + +//---------------------------------------------------------------------------- +void +PolynomialDiff( PNTVECTOR& vSou1, PNTVECTOR& vSou2, PNTVECTOR& vSumm) +{ + int nDeg1 = vSou1.size() - 1 ; + int nDeg2 = vSou2.size() - 1 ; + int nMin = (( nDeg1 <= nDeg2) ? nDeg1 : nDeg2) ; + int nMax = (( nDeg1 >= nDeg2) ? nDeg1 : nDeg2) ; + vSumm.clear() ; + vSumm.reserve( nMax + 1) ; + for ( int i = 0 ; i < nMin ; ++ i) + vSumm.push_back( vSou1[i] + ( - vSou2[i])) ; + if ( nDeg1 > nDeg2) { + for ( int i = nMin ; i < nDeg1 ; ++ i) + vSumm.push_back( vSou1[i]) ; + } + else if ( nDeg1 < nDeg2) { + for ( int i = nMin ; i < nDeg2 ; ++ i) + vSumm.push_back( Point3d() + ( - vSou2[i])) ; + } +} + +//---------------------------------------------------------------------------- +void +PolynomialMult( PNTVECTOR& vSou1, DBLVECTOR& vSou2, PNTVECTOR& vMult) +{ + int nDeg1 = vSou1.size() - 1 ; + int nDeg2 = vSou2.size() - 1 ; + int nDim = nDeg1 + nDeg2 + 1 ; + vMult.clear() ; + vMult.reserve( nDim) ; + for ( int i = 0 ; i < nDim ; ++ i) + vMult.push_back( Point3d()) ; + for ( int i = 0 ; i <= nDeg1 ; ++ i) { + for ( int j = 0 ; j <= nDeg2 ; ++ j) + vMult[i+j] += vSou1[i] * vSou2[j] ; + } +} + +//---------------------------------------------------------------------------- +void +PolynomialDerive( PNTVECTOR& vSou, PNTVECTOR& vDer) +{ + int nDeg = vSou.size() - 1 ; + vDer.clear() ; + vDer.reserve( nDeg) ; + for ( int i = 0 ; i < nDeg ; ++ i) + vDer.push_back( ( i + 1) * vSou[i+1]) ; +} + +//---------------------------------------------------------------------------- +void +PolynomialGetMainCompo( PNTVECTOR& vSou, DBLVECTOR& vCompo) +{ + int nDim = vSou.size() ; + vCompo.clear() ; + vCompo.reserve( nDim) ; + Point3d ptSumm ; + for ( int i = 0 ; i < nDim ; ++ i) { + ptSumm.x += fabs( vSou[i].x) ; + ptSumm.y += fabs( vSou[i].y) ; + ptSumm.z += fabs( vSou[i].z) ; + } + if ( ptSumm.x > ptSumm.y && ptSumm.x > ptSumm.z) { + for ( int i = 0 ; i < nDim ; ++ i) + vCompo.push_back( vSou[i].x) ; + } + else if ( ptSumm.y > ptSumm.z) { + for ( int i = 0 ; i < nDim ; ++ i) + vCompo.push_back( vSou[i].y) ; + } + else { + for ( int i = 0 ; i < nDim ; ++ i) + vCompo.push_back( vSou[i].z) ; + } +} +#endif \ No newline at end of file diff --git a/stdafx.h b/stdafx.h index 06eabfd..4d461a4 100644 --- a/stdafx.h +++ b/stdafx.h @@ -29,3 +29,4 @@ #include "/EgtDev/Include/EgtLibVer.h" #pragma comment(lib, EGTLIBDIR "EgtGeneral" EGTLIBVER ".lib") +#pragma comment(lib, EGTLIBDIR "EgtNumKernel" EGTLIBVER ".lib")