diff --git a/CalcDerivate.cpp b/CalcDerivate.cpp new file mode 100644 index 0000000..60f3582 --- /dev/null +++ b/CalcDerivate.cpp @@ -0,0 +1,233 @@ +//---------------------------------------------------------------------------- +// EgalTech 2026 +//---------------------------------------------------------------------------- +// File : CalcDerivate.cpp Data : 03.02.26 Versione : 1.5h1 +// Contenuto : Funzioni per calcolo derivate secondo Bessel e Akima. +// +// +// +// Modifiche : 03.02.26 DB Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +#pragma once + +#include "stdafx.h" +#include "CalcDerivate.h" +#include "/EgtDev/Include/EGkPoint3d.h" +#include "/EgtDev/Include/EgtNumCollection.h" +#include "/EgtDev/Include/EGkGeoCollection.h" + +//---------------------------------------------------------------------------- +bool +ComputeAkimaTangents( bool bDetectCorner, const DBLVECTOR& vPar, const PNTVECTOR& vPnt, VCT3DVECTOR& vPrevDer, VCT3DVECTOR& vNextDer) +{ + // pulisco i vettori dei parametri e delle tangenti + vPrevDer.clear() ; + vNextDer.clear() ; + + // numero di punti + int nSize = int( vPnt.size()) ; + + // sono necessari almeno due punti + if ( nSize < 2) + return false ; + + // calcolo le derivate + vPrevDer.reserve( nSize) ; + vNextDer.reserve( nSize) ; + // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce + if ( nSize == 2) { + // non esiste derivata prima del primo punto + vPrevDer.emplace_back( 0, 0, 0) ; + vNextDer.push_back( ( vPnt[1] - vPnt[0]) / ( vPar[1] - vPar[0])) ; + vPrevDer.push_back( vNextDer[0]) ; + // non esiste derivata dopo il secondo e ultimo punto + vNextDer.emplace_back( 0, 0, 0) ; + return true ; + } + // verifico se curva chiusa (primo e ultimo punto coincidono) + bool bClosed = AreSamePointApprox( vPnt.front(), vPnt.back()) ; + // calcolo le derivate + for ( int i = 0 ; i < nSize ; ++ i) { + Vector3d vtPrevDer ; + Vector3d vtNextDer ; + // primo punto + if ( i == 0) { + // se curva chiusa, come precedente uso il penultimo punto + if ( bClosed) { + // se non ci sono almeno 5 punti + if ( nSize < 5) { + if ( ! CalcCircleMidDer( vPar[nSize-2] - vPar[nSize-1], vPnt[nSize-2], vPar[i], vPnt[i], + vPar[i+1], vPnt[i+1], vtNextDer)) + return false ; + vtPrevDer = vtNextDer ; + } + // altrimenti + else { + if ( ! CalcAkimaMidDer( vPar[nSize-3] - vPar[nSize-1], vPnt[nSize-3], vPar[nSize-2] - vPar[nSize-1], vPnt[nSize-2], + vPar[i], vPnt[i], vPar[i+1], vPnt[i+1], + vPar[i+2], vPnt[i+2], bDetectCorner, + vtPrevDer, vtNextDer)) + return false ; + } + } + // altrimenti, uso arco sui primi tre punti + else { + if ( ! CalcCircleStartDer( vPar[i], vPnt[i], vPar[i+1], vPnt[i+1], + vPar[i+2], vPnt[i+2], vtNextDer)) + return false ; + vtPrevDer = Vector3d( 0, 0, 0) ; + } + } + // ultimo punto + else if ( i == nSize - 1) { + // se curva chiusa, le tg devono coincidere con quelle del primo + if ( bClosed) { + vtPrevDer = vPrevDer[0] ; + vtNextDer = vNextDer[0] ; + } + // altrimenti, uso arco sugli ultimi tre punti + else { + if ( ! CalcCircleEndDer( vPar[i-2], vPnt[i-2], vPar[i-1], vPnt[i-1], + vPar[i], vPnt[i], vtPrevDer)) + return false ; + vtNextDer = Vector3d( 0, 0, 0) ; + } + } + // punti intermedi + else { + // se secondo punto + if ( i == 1) { + // se curva aperta o non ci sono almeno 5 punti + if ( ! bClosed || nSize < 5) { + if ( ! CalcCircleMidDer( vPar[i-1], vPnt[i-1], vPar[i], vPnt[i], + vPar[i+1], vPnt[i+1], vtPrevDer)) + return false ; + vtNextDer = vtPrevDer ; + } + // altrimenti + else { + if ( ! CalcAkimaMidDer( vPar[nSize-2] - vPar[nSize-1], vPnt[nSize-2], vPar[i-1], vPnt[i-1], + vPar[i], vPnt[i], vPar[i+1], vPnt[i+1], + vPar[i+2], vPnt[i+2], bDetectCorner, + vtPrevDer, vtNextDer)) + return false ; + } + } + // se penultimo punto + else if ( i == nSize - 2) { + // se curva aperta o non ci sono almeno 5 punti + if ( ! bClosed || nSize < 5) { + if ( ! CalcCircleMidDer( vPar[i-1], vPnt[i-1], vPar[i], vPnt[i], + vPar[i+1], vPnt[i+1], vtPrevDer)) + return false ; + vtNextDer = vtPrevDer ; + } + // altrimenti + else { + if ( ! CalcAkimaMidDer( vPar[i-2], vPnt[i-2], vPar[i-1], vPnt[i-1], + vPar[i], vPnt[i], vPar[i+1], vPnt[i+1], + vPar[1] + vPar[i+1], vPnt[1], bDetectCorner, + vtPrevDer, vtNextDer)) + return false ; + } + } + // altrimenti + else { + if ( ! CalcAkimaMidDer( vPar[i-2], vPnt[i-2], vPar[i-1], vPnt[i-1], + vPar[i], vPnt[i], vPar[i+1], vPnt[i+1], + vPar[i+2], vPnt[i+2], bDetectCorner, + vtPrevDer, vtNextDer)) + return false ; + } + } + // salvo la derivata + vPrevDer.push_back( vtPrevDer) ; + vNextDer.push_back( vtNextDer) ; + } + + return true ; +} + +//---------------------------------------------------------------------------- +bool +ComputeBesselTangents( const DBLVECTOR& vPar, const PNTVECTOR& vPnt, VCT3DVECTOR& vPrevDer, VCT3DVECTOR& vNextDer) +{ + // pulisco i vettori dei parametri e delle tangenti + vPrevDer.clear() ; + vNextDer.clear() ; + + // numero di punti + int nSize = int( vPnt.size()) ; + + // sono necessari almeno due punti + if ( nSize < 2) + return false ; + + // calcolo le derivate + vPrevDer.reserve( nSize) ; + vNextDer.reserve( nSize) ; + // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce + if ( nSize == 2) { + // non esiste derivata prima del primo punto + vPrevDer.emplace_back( 0, 0, 0) ; + vNextDer.push_back( ( vPnt[1] - vPnt[0]) / ( vPar[1] - vPar[0])) ; + vPrevDer.push_back( vNextDer[0]) ; + // non esiste derivata dopo il secondo e ultimo punto + vNextDer.emplace_back( 0, 0, 0) ; + return true ; + } + // verifico se curva chiusa (primo e ultimo punto coincidono) + bool bClosed = AreSamePointApprox( vPnt.front(), vPnt.back()) ; + // calcolo le derivate + for ( int i = 0 ; i < nSize ; ++ i) { + Vector3d vtPrevDer ; + Vector3d vtNextDer ; + // primo punto + if ( i == 0) { + // se curva chiusa, come precedente uso il penultimo punto + if ( bClosed) { + if ( ! CalcBesselMidDer( vPar[nSize-2] - vPar[nSize-1], vPnt[nSize-2], vPar[i], vPnt[i], + vPar[i+1], vPnt[i+1], vtNextDer)) + return false ; + vtPrevDer = vtNextDer ; + } + // altrimenti, uso i primi tre punti + else { + if ( ! CalcBesselStartDer( vPar[i], vPnt[i], vPar[i+1], vPnt[i+1], + vPar[i+2], vPnt[i+2], vtNextDer)) + return false ; + vtPrevDer = Vector3d( 0, 0, 0) ; + } + } + // ultimo punto + else if ( i == nSize - 1) { + // se curva chiusa, le tg devono coincidere con quelle del primo + if ( bClosed) { + vtPrevDer = vPrevDer[0] ; + vtNextDer = vNextDer[0] ; + } + // altrimenti, uso gli ultimi tre punti + else { + if ( ! CalcBesselEndDer( vPar[i-2], vPnt[i-2], vPar[i-1], vPnt[i-1], + vPar[i], vPnt[i], vtPrevDer)) + return false ; + vtNextDer = Vector3d( 0, 0, 0) ; + } + } + // punti intermedi + else { + if ( ! CalcBesselMidDer( vPar[i-1], vPnt[i-1], vPar[i], vPnt[i], + vPar[i+1], vPnt[i+1], vtPrevDer)) + return false ; + vtNextDer = vtPrevDer ; + } + // salvo la derivata + vPrevDer.push_back( vtPrevDer) ; + vNextDer.push_back( vtNextDer) ; + } + + return true ; +} \ No newline at end of file diff --git a/CalcDerivate.h b/CalcDerivate.h index 1ba6904..d6a1bbd 100644 --- a/CalcDerivate.h +++ b/CalcDerivate.h @@ -14,6 +14,8 @@ #pragma once #include "/EgtDev/Include/EGkPoint3d.h" +#include "/EgtDev/Include/EgtNumCollection.h" +#include "/EgtDev/Include/EGkGeoCollection.h" //---------------------------------------------------------------------------- @@ -179,3 +181,6 @@ CalcAkimaMidDer( double dU0, const Point3d& ptP0, double dU1, const Point3d& ptP } return ( ! vtPrevDer.IsZero() && ! vtNextDer.IsZero()) ; } + +bool ComputeAkimaTangents( bool bDetectCorner, const DBLVECTOR& vPar, const PNTVECTOR& vPnt, VCT3DVECTOR& vPrevDer, VCT3DVECTOR& vNextDer) ; +bool ComputeBesselTangents( const DBLVECTOR& vPar, const PNTVECTOR& vPnt, VCT3DVECTOR& vPrevDer, VCT3DVECTOR& vNextDer) ; \ No newline at end of file diff --git a/CurveAux.cpp b/CurveAux.cpp index 5f4c1c3..7831265 100644 --- a/CurveAux.cpp +++ b/CurveAux.cpp @@ -13,6 +13,8 @@ //--------------------------- Include ---------------------------------------- #include "stdafx.h" +#include "CalcDerivate.h" +#include "Bernstein.h" #include "CurveAux.h" #include "GeoConst.h" #include "CurveLine.h" @@ -23,6 +25,7 @@ #include "IntersLineLine.h" #include "/EgtDev/Include/EGkDistPointCurve.h" #include "/EgtDev/Include/EGkStringUtils3d.h" +#include "/EgtDev/Include/EgtNumUtils.h" #include "/EgtDev/Include/EGkUiUnits.h" #include "/EgtDev/Include/EgtPointerOwner.h" #include "/EgtDev/Include/EGkCurveByInterp.h" @@ -30,6 +33,13 @@ #define EIGEN_NO_IO #include "/EgtDev/Extern/Eigen/Dense" +#define SAVEAPPROX 0 +#define SAVECURVEPASSED 0 +#define SAVELINEARAPPROX 0 +#if SAVEAPPROX || SAVECURVEPASSED || SAVELINEARAPPROX + #include "/EgtDev/Include/EGkGeoObjSave.h" +#endif + using namespace std ; static bool FindSpan( double dU, int nDeg, const DBLVECTOR& vKnots, int& nSpan) ; @@ -525,7 +535,8 @@ LineToBezierCurve( const ICurveLine* pCrvLine, int nDeg, bool bMakeRatOrNot) PtrOwner pCrvBezier( CreateCurveBezier()) ; // rendo tutte le curve di grado 2 e razionali così posso convertire anche archi e avere tutte curve dello stesso grado e razionali pCrvBezier->Init( nDeg, true) ; - pCrvBezier->FromLine( *pCrvLine) ; + if( ! pCrvBezier->FromLine( *pCrvLine)) + return nullptr ; if ( bMakeRatOrNot) pCrvBezier->MakeRational() ; return Release( pCrvBezier) ; @@ -1187,7 +1198,8 @@ CalcBasisFunc( double dU, int nSpan, int nDeg, const DBLVECTOR& vKnots, DBLVECTO //---------------------------------------------------------------------------- ICurve* -InterpolatePointSetWithBezierNoIntermedLines( const PNTVECTOR& vPnt, int nStart, int nEnd, int nDeg, const DBLVECTOR& vLen, double dLenTot) +InterpolatePointSetWithBezierNoIntermedLines( const PNTVECTOR& vPnt, int nStart, int nEnd, int nDeg, const DBLVECTOR& vLen, double dLenTot, + const Vector3d& vtStartDir = V_NULL, const Vector3d& vtEndDir = V_NULL) { PtrOwner pCrvInt ; @@ -1216,6 +1228,8 @@ InterpolatePointSetWithBezierNoIntermedLines( const PNTVECTOR& vPnt, int nStart, return nullptr ; } + bool bUseStartEndDir = vtStartDir.IsValid() && vtEndDir.IsValid() ; + DBLVECTOR vPntParam ; vPntParam.resize( nPoints) ; vPntParam[0] = 0 ; @@ -1224,13 +1238,15 @@ InterpolatePointSetWithBezierNoIntermedLines( const PNTVECTOR& vPnt, int nStart, vPntParam[i] = vPntParam[i-1] + vLen[i-1] / dLenTot ; DBLVECTOR vKnots ; - vKnots.resize( nPoints + nDeg - 1) ; + int nKnots = bUseStartEndDir ? nPoints + nDeg - 1 + 2 : nPoints + nDeg - 1 ; + vKnots.resize( nKnots) ; for ( int i = 0 ; i < nDeg ; ++i) { vKnots[i] = 0 ; vKnots.end()[-i-1] = 1 ; } - for ( int i = nDeg ; i < nPoints - 1 ; ++i) { + int nKnotsToEdit = bUseStartEndDir ? nPoints + 1 : nPoints - 1 ; + for ( int i = nDeg ; i < nKnotsToEdit ; ++i) { double dKnot = 0 ; for ( int j = i + 1 ; j < i + nDeg + 1 ; ++j) dKnot += vPntParam[j - nDeg] ; @@ -1238,13 +1254,14 @@ InterpolatePointSetWithBezierNoIntermedLines( const PNTVECTOR& vPnt, int nStart, vKnots[i] = dKnot ; } - Eigen::MatrixXd mA( nPoints, nPoints) ; + int nEq = bUseStartEndDir ? nPoints + 2 : nPoints ; + Eigen::MatrixXd mA( nEq, nEq) ; mA.fill( 0) ; - for ( int i = 0 ; i < nPoints ; ++i) { + for ( int i = 0 ; i < nEq ; ++i) { if ( i == 0) mA.row(0).col(0) << 1 ; - else if ( i == nPoints - 1) - mA.row(i).col(nPoints - 1) << 1 ; + else if ( i == nEq - 1) + mA.row(i).col( nEq - 1) << 1 ; else { int nSpan = 0 ; FindSpan( vPntParam[i], nDeg, vKnots, nSpan) ; DBLVECTOR vBasis ; vBasis.resize( nDeg + 1) ; @@ -1346,7 +1363,13 @@ InterpolatePointSetWithBezier( const PNTVECTOR& vPnt, double dLinTol, double dMa if ( vLen.size() != 0) { if ( nEnd == 0) nEnd = nPoints - 1 ; - pCrvInt->AddCurve( InterpolatePointSetWithBezierNoIntermedLines( vPnt, nStart, nEnd, nDeg, vLen, dLenTot)) ; + Vector3d vtStartDir = V_INVALID ; + Vector3d vtEndDir = V_INVALID ; + //if ( nStart != 0 && nEnd != nPoints - 1) { + // pCrvInt->GetEndDir( vtStartDir) ; + // vtEndDir = + //} + pCrvInt->AddCurve( InterpolatePointSetWithBezierNoIntermedLines( vPnt, nStart, nEnd, nDeg, vLen, dLenTot, vtStartDir, vtEndDir)) ; } if ( bFoundLine) { @@ -1358,6 +1381,7 @@ InterpolatePointSetWithBezier( const PNTVECTOR& vPnt, double dLinTol, double dMa nStart = nEnd ; } + //dErr = 0 ; CalcApproxError( pCrvOri, pCrvInt, dErr) ; if ( dErr > dLinTol && dMaxLen > 200 * EPS_SMALL) dMaxLen /= 2 ; @@ -1372,42 +1396,282 @@ InterpolatePointSetWithBezier( const PNTVECTOR& vPnt, double dLinTol, double dMa } //---------------------------------------------------------------------------- -ICurve* -ApproxCurveWithBezier( const ICurve* pCrv , double dTol, int nType) +static bool +ParamByLen( const PNTVECTOR& vPnt, DBLVECTOR& vParam, int nFirst, int nLast) { - PolyLine plApprox ; - double dAngTolFine = 2 ; - pCrv->ApproxWithLines( dTol, dAngTolFine, ICurve::APL_STD, plApprox) ; + int nPoints = nLast - nFirst + 1 ; + if ( nPoints < 2) + return false ; + if( ssize(vParam) == 0) + vParam.resize( nPoints) ; + if( vParam[nFirst] == 0 && vParam[nLast] == 1) + return true ; + vParam[nFirst] = 0 ; + for ( int i = nFirst + 1 ; i <= nLast ; ++i) { + double dDist = Dist( vPnt[i], vPnt[i-1]) ; + vParam[i] = vParam[i- 1] + dDist ; + } + for ( int i = nFirst + 1 ; i < nLast ; ++i) + vParam[i] /= vParam[nLast] ; + vParam[nLast] = 1 ; - PNTVECTOR vPnt ; - Point3d pt ; plApprox.GetFirstPoint( pt) ; - do { - vPnt.push_back( pt) ; - } while ( plApprox.GetNextPoint( pt)) ; + return true ; +} - // campiono punti lungo la curva e poi li interpolo +//---------------------------------------------------------------------------- +ICurveBezier* +ApproxPointSetWithSingleBezier( const PNTVECTOR& vPnt, int nFirst, int nLast, const Vector3d& vtStartDir, const Vector3d& vtEndDir, + const DBLVECTOR& vParam) +{ + // cerco di approssimare un set di punti con una sola bezier cubica non razionale + int nPoints = nLast - nFirst + 1 ; + // se ho solo quattro punti allora costruisco direttamente la curva + PtrOwner pCrvBez( CreateCurveBezier()) ; + int nDeg = 3 ; + bool bRat = false ; + pCrvBez->Init( nDeg, bRat) ; + const Point3d& pt0 = vPnt[nFirst] ; + const Point3d& pt3 = vPnt[nLast] ; + pCrvBez->SetControlPoint( 0, pt0) ; + pCrvBez->SetControlPoint( 3, pt3) ; + Eigen::Vector2d mA ; + if( nPoints > 4) { + // risoluzione sistema + Eigen::Matrix2d mC ; mC.setZero() ; + Eigen::Vector2d mX ; mX.setZero() ; + for ( int i = nFirst ; i <= nLast ; ++i) { + double dU = vParam[i] ; + DBLVECTOR vBern(4) ; + GetAllBernstein( dU, 3, vBern) ; - PtrOwner pCC( InterpolatePointSetWithBezier( vPnt, dTol, 100)) ; - if ( ! IsNull( pCC) && pCC->IsValid()) - return Release( pCC) ; - else - return nullptr ; + Vector3d A1 = vtStartDir * vBern[1] ; + Vector3d A2 = vtEndDir * vBern[2] ; + + Vector3d tmp = vPnt[i] - ( pt0 * ( vBern[0] + vBern[1])) - (( pt3 * ( vBern[2] + vBern[3])) - ORIG) ; // ORIG serve solo per trasformare l'ultimo termine in un vettore + + mC(0,0) += A1 * A1 ; + mC(0,1) += A1 * A2 ; + mC(1,0) += A1 * A2 ; + mC(1,1) += A2 * A2 ; + + mX(0) += A1 * tmp ; + mX(1) += A2 * tmp ; + } + mA = mC.fullPivLu().solve(mX) ; + } + // l'algoritmo è fatto in modo che alpha1 e alpha2 siano positivi ( se tutto va bene) + // io invece ho tenuto le tangenti con la direzione originale, quindi il primo dovrebbe essere positivo e il secondo negativo + if ( mA(0) < 0 || mA(1) > 0 || nPoints < 4) { + if( mA(0) < 0 || mA(1) > 0) + LOG_DBG_ERR( GetEGkLogger(), "valori di alfa sballati, potrebbe essere la spaziatura dismogenea tra punti") + double dDistCorr = Dist( pt3, pt0) / 3 ; + mA(0) = dDistCorr ; + mA(1) = - dDistCorr ; + } + + double alpah1 = mA(0) ; + double alpha2 = mA(1) ; + Point3d pt1 = pt0 + vtStartDir * mA(0) ; + Point3d pt2 = pt3 + vtEndDir * mA(1) ; + pCrvBez->SetControlPoint( 1, pt1) ; + pCrvBez->SetControlPoint( 2, pt2) ; + return Release( pCrvBez) ; +} + +//---------------------------------------------------------------------------- +bool +CalcPointSetApproxError( const PNTVECTOR& vPntOrig, const DBLVECTOR& vParam, int nFirst, int nLast, const ICurve* pCrvNew, double& dErr, int& nPointMaxErr) +{ + dErr = 0 ; + // calcolo l'errore di approssimazione + for ( int i = nFirst ; i <= nLast ; ++i) { + Point3d ptBez ; pCrvNew->GetPointD1D2( vParam[i], ICurve::Side::FROM_MINUS, ptBez) ; + double dErrTemp = Dist( vPntOrig[i], ptBez) ; + if( dErrTemp > dErr) { + dErr = dErrTemp ; + nPointMaxErr = i ; + } + } + + return true ; } //---------------------------------------------------------------------------- ICurve* -ApproxPointSetWithBezier( const ICurve* pCrv, double dTol) +FitWithBezier( const ICurve* pCrvOrig, const PNTVECTOR& vPnt, DBLVECTOR& vParam, int nFirst, int nLast, const VCT3DVECTOR& vPrevDer, const VCT3DVECTOR& vNextDer, double dTol) { - // campiono punti lungo la curva e poi li interpolo + ParamByLen(vPnt, vParam, nFirst, nLast) ; + PtrOwner pCrvFit( CreateCurveComposite()) ; + PtrOwner pCrvBez( CreateCurveBezier()) ; + double dErr = INFINITO ; + double dErrSplit = dTol * 25 ; + int nIter = 0 ; + while ( dErr > dTol && nIter < 10) { + if( dErr < INFINITO) { + // riparametrizzo i punti + for ( int i = nFirst + 1 ; i < nLast ; ++i) { + //questo potrebbe diventare un while appena capisco di quanto si aggiusta il parametro ad ogni iterazione + double dCorr = 1 ; + do { + Vector3d vtDer1, vtDer2 ; + Point3d ptBez ; pCrvBez->GetPointD1D2( vParam[i], ICurve::Side::FROM_MINUS, ptBez, &vtDer1, &vtDer2) ; + Vector3d vtLink = ptBez - vPnt[i] ; + double dNum = vtLink * vtDer1 ; + double dDen = vtLink * vtDer2 + vtDer1 * vtDer1 ; + dCorr = dDen > EPS_ZERO ? dNum / dDen : 0 ; + vParam[i] = vParam[i] - dCorr ; + } while ( abs( dCorr) > EPS_ZERO) ; + Clamp( vParam[i], 0., 1.) ; + } + } - PtrOwner pCC( CreateBasicCurveComposite()) ; - return Release( pCC) ; + // fit della curva + pCrvBez.Set( ApproxPointSetWithSingleBezier( vPnt, nFirst, nLast, vNextDer[nFirst], vPrevDer[nLast], vParam)) ; + if( IsNull( pCrvBez) || ! pCrvBez->IsValid()) + return nullptr ; +#if SAVEAPPROX + SaveGeoObj( pCrvBez->Clone(), "D:\\Temp\\bezier\\approxWithBezier\\first_approx.nge") ; +#endif + + int nPointMaxErr = 0 ; + CalcPointSetApproxError( vPnt, vParam, nFirst, nLast, pCrvBez, dErr, nPointMaxErr) ; + // se sto unendo due punti consecutivi e l'errore è oltre quello richiesto allora restituisco un segmento che unisce i punti + if( nLast - nFirst == 1 && dErr > dTol) { + CurveLine CL ; CL.Set( vPnt[nFirst], vPnt[nLast]) ; + pCrvBez.Set( GetCurveBezier( CurveToBezierCurve( &CL))) ; + dErr = 0 ; + } + + ++nIter ; + bool bSplit = false ; + if( nIter == 10 && dErr > dTol) + bSplit = true ; + // se la curva di approssimazione è ancora molto lontana dalla curva originale allora divido il set di punti in due + if ( dErr > dErrSplit || bSplit) { + if ( nLast - nFirst > 1) { + if( ! pCrvFit->AddCurve( FitWithBezier( pCrvOrig, vPnt, vParam, nFirst, nPointMaxErr, vPrevDer, vNextDer,dTol)) || + ! pCrvFit->AddCurve( FitWithBezier( pCrvOrig, vPnt, vParam, nPointMaxErr, nLast, vPrevDer, vNextDer, dTol))) + return nullptr ; + break ; + } + else + return nullptr ; + } + } + + if ( pCrvFit->GetCurveCount() > 0) + return Release( pCrvFit) ; + else if ( dErr < dTol && ! IsNull( pCrvBez) && pCrvBez->IsValid()) + return Release( pCrvBez) ; + else + return nullptr ; +} + +static int nCrvPassed = 0 ; + +//---------------------------------------------------------------------------- +ICurve* +ApproxCurveWithBezier( const ICurve* pCrv , double dTol, int nType) +{ + +#if SAVECURVEPASSED + SaveGeoObj( pCrv->Clone(), "D:\\Temp\\bezier\\approxWithBezier\\CurveDaApprossimare\\"+ToString(nCrvPassed) + ".nge") ; + ++nCrvPassed ; +#endif + + //// uso l'algoritmo di Schneider in Grafic Gems I + // mi aspetto che non ci siano angoli ( discontinuità della derivata prima) nel risultato desiderato + PolyLine plApprox ; + double dAngTolFine = 1 ; + double dLinTolFine = 0.05 ; + pCrv->ApproxWithLines( dLinTolFine, dAngTolFine, ICurve::APL_STD, plApprox) ; + +#if SAVELINEARAPPROX + CurveComposite CC ; CC.FromPolyLine(plApprox) ; + SaveGeoObj( CC.Clone(), "D:\\Temp\\bezier\\approxWithBezier\\approssimazione_lineare.nge") ; +#endif + + CurveByInterp cbi ; + PNTVECTOR vPnt ; + Point3d pt ; plApprox.GetFirstPoint( pt) ; + do { + if ( ssize(vPnt) != 0) { + vPnt.push_back( Media(vPnt.back(), pt,1./3.)) ; + vPnt.push_back( Media(vPnt.back(), pt,2./3.)) ; + } + vPnt.push_back( pt) ; + cbi.AddPoint( pt) ; + } while ( plApprox.GetNextPoint( pt)) ; + + // calcolo la curvatura nei vari punti per identificare zone a curvatura costante + DBLVECTOR vCrvVal( ssize(vPnt)) ; + DBLVECTOR vRad ( ssize( vPnt)) ; + for( int i = 1 ; i < ssize( vPnt) - 1 ; ++i) { + Vector3d vtA = vPnt[i] - vPnt[i-1] ; + Vector3d vtB = vPnt[i+1] - vPnt[i-1] ; + double dR = ( vtA.Len() * vtB.Len() * ( vtA - vtB).Len()) / ( 2 * (vtA ^ vtB).Len()) ; + double dCurvature = dR > 0 ? 1. / dR : 0 ; + vRad[i] = dR ; + vCrvVal[i] = dCurvature ; + } + vCrvVal[0] = vCrvVal[1] ; + vCrvVal.back() = vCrvVal.end()[-2] ; + // identifico le zone a curvatura costante // primo e ultimo punto degli intervalli non devono avere curvatura uguale agli altri perché sono di frontiera + INTINTVECTOR vConstCurv ; + double dCurvaturePrec = vCrvVal[0] ; + int nStart = 0 ; + int nEnd = 1 ; + while ( nStart < ssize( vPnt) - 1) { + dCurvaturePrec = vCrvVal[nEnd] ; + while ( nEnd < ssize( vPnt) - 1 && abs( vCrvVal[nEnd] - dCurvaturePrec) < 0.1) { + dCurvaturePrec = vCrvVal[nEnd] ; + ++nEnd ; + } + vConstCurv.emplace_back( nStart, nEnd) ; + nStart = nEnd ; + ++nEnd ; + } + if( ssize(vConstCurv) == 0) + vConstCurv.emplace_back( 0, ssize( vPnt) - 1) ; + + int nPoints = ssize( vPnt) ; + DBLVECTOR vParam( nPoints) ; + int nFirst = 0 ; + int nLast = ssize( vPnt) - 1 ; + ParamByLen( vPnt, vParam, nFirst, nLast) ; + + VCT3DVECTOR vPrevDer ; + VCT3DVECTOR vNextDer ; + ComputeAkimaTangents( false, vParam, vPnt, vPrevDer, vNextDer) ; + //normalizzo tutte le derivate + for( int i = 0 ; i < ssize( vPrevDer) ; ++i) { + vPrevDer[i].Normalize() ; + vNextDer[i].Normalize() ; + } + + PtrOwner pCCApproxTot( CreateCurveComposite()) ; + for( INTINT iiSE : vConstCurv) { + nFirst = iiSE.first ; + nLast = iiSE.second ; + //definisco la bezier che vado a raffinare iterativamente + PtrOwner pCApprox( FitWithBezier( pCrv, vPnt, vParam, nFirst, nLast, vPrevDer, vNextDer, dTol)) ; + if( IsNull( pCApprox) || ! pCApprox->IsValid()) + return nullptr ; + if( ! pCCApproxTot->AddCurve( Release( pCApprox))) + return nullptr ; + } + + return Release( pCCApproxTot) ; } //---------------------------------------------------------------------------- bool CalcApproxError( const ICurve* pCrvOri, const ICurve* pCrvNew, double& dErr, int nPoints) { + if( pCrvOri == nullptr || ! pCrvOri->IsValid() || pCrvNew == nullptr || ! pCrvNew->IsValid()){ + dErr = INFINITO ; + return false ; + } // controllo l'errore effettivo campionando più finemente double dLenOri = 0 ; pCrvOri->GetLength( dLenOri) ; double dLenNew = 0 ; pCrvNew->GetLength( dLenNew) ; diff --git a/CurveByApprox.cpp b/CurveByApprox.cpp index 836ccfe..2e3307c 100644 --- a/CurveByApprox.cpp +++ b/CurveByApprox.cpp @@ -202,213 +202,14 @@ CurveByApprox::CalcParameterization( void) bool CurveByApprox::CalcAkimaTangents( bool bDetectCorner) { - // pulisco i vettori delle tangenti - m_vPrevDer.clear() ; - m_vNextDer.clear() ; - - // numero di punti - int nSize = int( m_vPnt.size()) ; - - // sono necessari almeno due punti - if ( nSize < 2) - return false ; - - // calcolo le derivate - m_vPrevDer.reserve( nSize) ; - m_vNextDer.reserve( nSize) ; - // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce - if ( nSize == 2) { - // non esiste derivata prima del primo punto - m_vPrevDer.emplace_back( 0, 0, 0) ; - m_vNextDer.push_back( ( m_vPnt[1] - m_vPnt[0]) / ( m_vPar[1] - m_vPar[0])) ; - m_vPrevDer.push_back( m_vNextDer[0]) ; - // non esiste derivata dopo il secondo e ultimo punto - m_vNextDer.emplace_back( 0, 0, 0) ; - return true ; - } - // verifico se curva chiusa (primo e ultimo punto coincidono) - bool bClosed = AreSamePointApprox( m_vPnt.front(), m_vPnt.back()) ; - // calcolo le derivate - for ( int i = 0 ; i < nSize ; ++ i) { - Vector3d vtPrevDer ; - Vector3d vtNextDer ; - // primo punto - if ( i == 0) { - // se curva chiusa, come precedente uso il penultimo punto - if ( bClosed) { - // se non ci sono almeno 5 punti - if ( nSize < 5) { - if ( ! CalcCircleMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtNextDer)) - return false ; - vtPrevDer = vtNextDer ; - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[nSize-3] - m_vPar[nSize-1], m_vPnt[nSize-3], m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // altrimenti, uso arco sui primi tre punti - else { - if ( ! CalcCircleStartDer( m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], vtNextDer)) - return false ; - vtPrevDer = Vector3d( 0, 0, 0) ; - } - } - // ultimo punto - else if ( i == nSize - 1) { - // se curva chiusa, le tg devono coincidere con quelle del primo - if ( bClosed) { - vtPrevDer = m_vPrevDer[0] ; - vtNextDer = m_vNextDer[0] ; - } - // altrimenti, uso arco sugli ultimi tre punti - else { - if ( ! CalcCircleEndDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], vtPrevDer)) - return false ; - vtNextDer = Vector3d( 0, 0, 0) ; - } - } - // punti intermedi - else { - // se secondo punto - if ( i == 1) { - // se curva aperta o non ci sono almeno 5 punti - if ( ! bClosed || nSize < 5) { - if ( ! CalcCircleMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtPrevDer)) - return false ; - vtNextDer = vtPrevDer ; - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // se penultimo punto - else if ( i == nSize - 2) { - // se curva aperta o non ci sono almeno 5 punti - if ( ! bClosed || nSize < 5) { - if ( ! CalcCircleMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtPrevDer)) - return false ; - vtNextDer = vtPrevDer ; - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[1] + m_vPar[i+1], m_vPnt[1], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // salvo la derivata - m_vPrevDer.push_back( vtPrevDer) ; - m_vNextDer.push_back( vtNextDer) ; - } - - return true ; + return ComputeAkimaTangents( bDetectCorner, m_vPar, m_vPnt, m_vPrevDer, m_vNextDer) ; } //---------------------------------------------------------------------------- bool CurveByApprox::CalcBesselTangents( void) { - // pulisco i vettori delle tangenti - m_vPrevDer.clear() ; - m_vNextDer.clear() ; - - // numero di punti - int nSize = int( m_vPnt.size()) ; - - // sono necessari almeno due punti - if ( nSize < 2) - return false ; - - // calcolo le derivate - m_vPrevDer.reserve( nSize) ; - m_vNextDer.reserve( nSize) ; - // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce - if ( nSize == 2) { - // non esiste derivata prima del primo punto - m_vPrevDer.emplace_back( 0, 0, 0) ; - m_vNextDer.push_back( ( m_vPnt[1] - m_vPnt[0]) / ( m_vPar[1] - m_vPar[0])) ; - m_vPrevDer.push_back( m_vNextDer[0]) ; - // non esiste derivata dopo il secondo e ultimo punto - m_vNextDer.emplace_back( 0, 0, 0) ; - return true ; - } - // verifico se curva chiusa (primo e ultimo punto coincidono) - bool bClosed = AreSamePointApprox( m_vPnt.front(), m_vPnt.back()) ; - // calcolo le derivate - for ( int i = 0 ; i < nSize ; ++ i) { - Vector3d vtPrevDer ; - Vector3d vtNextDer ; - // primo punto - if ( i == 0) { - // se curva chiusa, come precedente uso il penultimo punto - if ( bClosed) { - if ( ! CalcBesselMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtNextDer)) - return false ; - vtPrevDer = vtNextDer ; - } - // altrimenti, uso i primi tre punti - else { - if ( ! CalcBesselStartDer( m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], vtNextDer)) - return false ; - vtPrevDer = Vector3d( 0, 0, 0) ; - } - } - // ultimo punto - else if ( i == nSize - 1) { - // se curva chiusa, le tg devono coincidere con quelle del primo - if ( bClosed) { - vtPrevDer = m_vPrevDer[0] ; - vtNextDer = m_vNextDer[0] ; - } - // altrimenti, uso gli ultimi tre punti - else { - if ( ! CalcBesselEndDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], vtPrevDer)) - return false ; - vtNextDer = Vector3d( 0, 0, 0) ; - } - } - // punti intermedi - else { - if ( ! CalcBesselMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtPrevDer)) - return false ; - vtNextDer = vtPrevDer ; - } - // salvo la derivata - m_vPrevDer.push_back( vtPrevDer) ; - m_vNextDer.push_back( vtNextDer) ; - } - - return true ; + return ComputeBesselTangents( m_vPar, m_vPnt, m_vPrevDer, m_vNextDer) ; } //---------------------------------------------------------------------------- diff --git a/CurveByInterp.cpp b/CurveByInterp.cpp index 0c7003c..6ba70b6 100644 --- a/CurveByInterp.cpp +++ b/CurveByInterp.cpp @@ -50,12 +50,42 @@ CurveByInterp::AddPoint( const Point3d& ptP) ICurve* CurveByInterp::GetCurve( int nMethod, int nType) { - // calcolo le tangenti + // se richieste curve di Bezier cubiche (ottenute da interpolazione con Nurbs) + if ( nType == CUBIC_BEZIERS_LONG) { + // creo la curva composita + PtrOwner pCrv ; + //pCrv.Set( InterpolatePointSetWithBezier( m_vPnt, 50 * EPS_SMALL, 50)) ; + //debug + pCrv.Set( InterpolatePointSetWithBezier( m_vPnt, 0.1, 100)) ; + if ( IsNull(pCrv) || ! pCrv->IsValid()) + return nullptr ; + return Release( pCrv) ; + } + + // numero di punti + int nSize = int( m_vPnt.size()) ; + + // sono necessari almeno due punti + if ( nSize < 2) + return nullptr ; + + // calcolo le distanze tra i punti per derivarne i parametri + m_vPar.reserve( nSize) ; + double dPar = 0 ; + m_vPar.push_back( dPar) ; + + for ( int i = 1 ; i < nSize ; ++ i) { + double dDist = Dist( m_vPnt[i-1], m_vPnt[i]) ; + dPar += dDist ; + m_vPar.push_back( dPar) ; + } + + // calcolo le tangenti if ( nMethod == BESSEL) { if ( ! CalcBesselTangents()) return nullptr ; } - else if ( nType != CUBIC_BEZIERS_LONG) { + else { if ( ! CalcAkimaTangents( nMethod == AKIMA_CORNER)) return nullptr ; } @@ -103,16 +133,6 @@ CurveByInterp::GetCurve( int nMethod, int nType) return ::Release( pCrvCompo) ; } - // se richieste curve di Bezier cubiche (ottenute da interpolazione con Nurbs) - if ( nType == CUBIC_BEZIERS_LONG) { - // creo la curva composita - PtrOwner pCrv ; - pCrv.Set( InterpolatePointSetWithBezier( m_vPnt, 50 * EPS_SMALL, 50)) ; - if ( IsNull(pCrv) || ! pCrv->IsValid()) - return nullptr ; - return Release( pCrv) ; - } - return nullptr ; } @@ -120,233 +140,12 @@ CurveByInterp::GetCurve( int nMethod, int nType) bool CurveByInterp::CalcAkimaTangents( bool bDetectCorner) { - // pulisco i vettori dei parametri e delle tangenti - m_vPar.clear() ; - m_vPrevDer.clear() ; - m_vNextDer.clear() ; - - // numero di punti - int nSize = int( m_vPnt.size()) ; - - // sono necessari almeno due punti - if ( nSize < 2) - return false ; - - // calcolo le distanze tra i punti per derivarne i parametri - m_vPar.reserve( nSize) ; - double dPar = 0 ; - m_vPar.push_back( dPar) ; - for ( int i = 1 ; i < nSize ; ++ i) { - double dDist = Dist( m_vPnt[i-1], m_vPnt[i]) ; - dPar += dDist ; - m_vPar.push_back( dPar) ; - } - - // calcolo le derivate - m_vPrevDer.reserve( nSize) ; - m_vNextDer.reserve( nSize) ; - // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce - if ( nSize == 2) { - // non esiste derivata prima del primo punto - m_vPrevDer.emplace_back( 0, 0, 0) ; - m_vNextDer.push_back( ( m_vPnt[1] - m_vPnt[0]) / ( m_vPar[1] - m_vPar[0])) ; - m_vPrevDer.push_back( m_vNextDer[0]) ; - // non esiste derivata dopo il secondo e ultimo punto - m_vNextDer.emplace_back( 0, 0, 0) ; - return true ; - } - // verifico se curva chiusa (primo e ultimo punto coincidono) - bool bClosed = AreSamePointApprox( m_vPnt.front(), m_vPnt.back()) ; - // calcolo le derivate - for ( int i = 0 ; i < nSize ; ++ i) { - Vector3d vtPrevDer ; - Vector3d vtNextDer ; - // primo punto - if ( i == 0) { - // se curva chiusa, come precedente uso il penultimo punto - if ( bClosed) { - // se non ci sono almeno 5 punti - if ( nSize < 5) { - if ( ! CalcCircleMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtNextDer)) - return false ; - vtPrevDer = vtNextDer ; - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[nSize-3] - m_vPar[nSize-1], m_vPnt[nSize-3], m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // altrimenti, uso arco sui primi tre punti - else { - if ( ! CalcCircleStartDer( m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], vtNextDer)) - return false ; - vtPrevDer = Vector3d( 0, 0, 0) ; - } - } - // ultimo punto - else if ( i == nSize - 1) { - // se curva chiusa, le tg devono coincidere con quelle del primo - if ( bClosed) { - vtPrevDer = m_vPrevDer[0] ; - vtNextDer = m_vNextDer[0] ; - } - // altrimenti, uso arco sugli ultimi tre punti - else { - if ( ! CalcCircleEndDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], vtPrevDer)) - return false ; - vtNextDer = Vector3d( 0, 0, 0) ; - } - } - // punti intermedi - else { - // se secondo punto - if ( i == 1) { - // se curva aperta o non ci sono almeno 5 punti - if ( ! bClosed || nSize < 5) { - if ( ! CalcCircleMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtPrevDer)) - return false ; - vtNextDer = vtPrevDer ; - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // se penultimo punto - else if ( i == nSize - 2) { - // se curva aperta o non ci sono almeno 5 punti - if ( ! bClosed || nSize < 5) { - if ( ! CalcCircleMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtPrevDer)) - return false ; - vtNextDer = vtPrevDer ; - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[1] + m_vPar[i+1], m_vPnt[1], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // altrimenti - else { - if ( ! CalcAkimaMidDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], bDetectCorner, - vtPrevDer, vtNextDer)) - return false ; - } - } - // salvo la derivata - m_vPrevDer.push_back( vtPrevDer) ; - m_vNextDer.push_back( vtNextDer) ; - } - - return true ; + return ComputeAkimaTangents( bDetectCorner, m_vPar, m_vPnt, m_vPrevDer, m_vNextDer) ; } //---------------------------------------------------------------------------- bool CurveByInterp::CalcBesselTangents( void) { - // pulisco i vettori dei parametri e delle tangenti - m_vPar.clear() ; - m_vPrevDer.clear() ; - m_vNextDer.clear() ; - - // numero di punti - int nSize = int( m_vPnt.size()) ; - - // sono necessari almeno due punti - if ( nSize < 2) - return false ; - - // calcolo le distanze tra i punti per derivarne i parametri - m_vPar.reserve( nSize) ; - double dPar = 0 ; - m_vPar.push_back( dPar) ; - for ( int i = 1 ; i < nSize ; ++ i) { - double dDist = Dist( m_vPnt[i-1], m_vPnt[i]) ; - dPar += dDist ; - m_vPar.push_back( dPar) ; - } - - // calcolo le derivate - m_vPrevDer.reserve( nSize) ; - m_vNextDer.reserve( nSize) ; - // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce - if ( nSize == 2) { - // non esiste derivata prima del primo punto - m_vPrevDer.emplace_back( 0, 0, 0) ; - m_vNextDer.push_back( ( m_vPnt[1] - m_vPnt[0]) / ( m_vPar[1] - m_vPar[0])) ; - m_vPrevDer.push_back( m_vNextDer[0]) ; - // non esiste derivata dopo il secondo e ultimo punto - m_vNextDer.emplace_back( 0, 0, 0) ; - return true ; - } - // verifico se curva chiusa (primo e ultimo punto coincidono) - bool bClosed = AreSamePointApprox( m_vPnt.front(), m_vPnt.back()) ; - // calcolo le derivate - for ( int i = 0 ; i < nSize ; ++ i) { - Vector3d vtPrevDer ; - Vector3d vtNextDer ; - // primo punto - if ( i == 0) { - // se curva chiusa, come precedente uso il penultimo punto - if ( bClosed) { - if ( ! CalcBesselMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtNextDer)) - return false ; - vtPrevDer = vtNextDer ; - } - // altrimenti, uso i primi tre punti - else { - if ( ! CalcBesselStartDer( m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1], - m_vPar[i+2], m_vPnt[i+2], vtNextDer)) - return false ; - vtPrevDer = Vector3d( 0, 0, 0) ; - } - } - // ultimo punto - else if ( i == nSize - 1) { - // se curva chiusa, le tg devono coincidere con quelle del primo - if ( bClosed) { - vtPrevDer = m_vPrevDer[0] ; - vtNextDer = m_vNextDer[0] ; - } - // altrimenti, uso gli ultimi tre punti - else { - if ( ! CalcBesselEndDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1], - m_vPar[i], m_vPnt[i], vtPrevDer)) - return false ; - vtNextDer = Vector3d( 0, 0, 0) ; - } - } - // punti intermedi - else { - if ( ! CalcBesselMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i], - m_vPar[i+1], m_vPnt[i+1], vtPrevDer)) - return false ; - vtNextDer = vtPrevDer ; - } - // salvo la derivata - m_vPrevDer.push_back( vtPrevDer) ; - m_vNextDer.push_back( vtNextDer) ; - } - - return true ; -} + return ComputeBesselTangents( m_vPar, m_vPnt, m_vPrevDer, m_vNextDer) ; +} \ No newline at end of file diff --git a/EgtGeomKernel.vcxproj b/EgtGeomKernel.vcxproj index c571c4d..202867f 100644 --- a/EgtGeomKernel.vcxproj +++ b/EgtGeomKernel.vcxproj @@ -281,6 +281,7 @@ copy $(TargetPath) \EgtProg\Dll64 + diff --git a/EgtGeomKernel.vcxproj.filters b/EgtGeomKernel.vcxproj.filters index 06c76e0..0c644f1 100644 --- a/EgtGeomKernel.vcxproj.filters +++ b/EgtGeomKernel.vcxproj.filters @@ -567,6 +567,9 @@ File di origine\GeoStriping + + File di origine\Geo +