Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'NewApproxWithBezier'

Browse Source
This commit is contained in:
Daniele Bariletti
2026-02-06 17:28:46 +01:00
parent f111827545 2023ba3f7d
commit 64c3363426
7 changed files with 572 additions and 466 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CalcDerivate.cpp
+233
View File
@@ -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 ;
}
CalcDerivate.h
+5
View File
@@ -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) ;
CurveAux.cpp
+293 -29
View File
@@ -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<ICurveBezier> 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<ICurve> 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<ICurveBezier> 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<ICurve> 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<ICurveComposite> pCrvFit( CreateCurveComposite()) ;
PtrOwner<ICurveBezier> 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<ICurveComposite> 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<ICurveComposite> pCCApproxTot( CreateCurveComposite()) ;
for( INTINT iiSE : vConstCurv) {
nFirst = iiSE.first ;
nLast = iiSE.second ;
//definisco la bezier che vado a raffinare iterativamente
PtrOwner<ICurve> 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) ;
CurveByApprox.cpp
+2 -201
View File
@@ -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) ;
}
//----------------------------------------------------------------------------
CurveByInterp.cpp
+35 -236
View File
@@ -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<ICurve> 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<ICurve> 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) ;
}
EgtGeomKernel.vcxproj
+1
View File
@@ -281,6 +281,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
<ClCompile Include="BBox3d.cpp" />
<ClCompile Include="BiArcs.cpp" />
<ClCompile Include="CalcPocketing.cpp" />
<ClCompile Include="CalcDerivate.cpp" />
<ClCompile Include="CAvSilhouetteSurfTm.cpp" />
<ClCompile Include="CAvSimpleSurfFrMove.cpp" />
<ClCompile Include="CAvToolSurfTm.cpp" />
EgtGeomKernel.vcxproj.filters
+3
View File
@@ -567,6 +567,9 @@
<ClCompile Include="Trimming.cpp">
<Filter>File di origine\GeoStriping</Filter>
</ClCompile>
<ClCompile Include="CalcDerivate.cpp">
<Filter>File di origine\Geo</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<ClInclude Include="stdafx.h">