Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
b4b996dac05aa2170812171f10a55586bb34d0a9
EgtGeomKernel/CurveByApprox.cpp
T
Dario Sassi b4b996dac0 EgtGeomKernel 1.6g7 : 
- aggiunta approssimazione di punti con archi e rette (CurveByApprox)
- fatte correzioni ad intersezioni rette/archi e archi/archi quasi tangenti
- correzioni ad offset di curve composite che non liberava memoria in caso di errore.
2015-07-27 13:32:02 +00:00

482 lines
17 KiB
C++
Raw Blame History

//----------------------------------------------------------------------------
// EgalTech 2015-2015
//----------------------------------------------------------------------------
// File : CurveByApprox.cpp Data : 23.07.15 Versione : 1.6g7
// Contenuto : Implementazione della classe CurveByApprox, per creare
// una curva mediante approssimazione di punti.
//
//
// Modifiche : 23.07.15 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------
//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CurveComposite.h"
#include "CalcDerivate.h"
#include "BiArcs.h"
#include "DistPointLine.h"
#include "/EgtDev/Include/EGkCurveByApprox.h"
#include "/EgtDev/Include/EGkPolyLine.h"
#include "/EgtDev/Include/EGkPolyArc.h"
#include "/EgtDev/Include/EgtPointerOwner.h"
#include <algorithm>
using namespace std ;
//----------------------------------------------------------------------------
bool
CurveByApprox::Reset( void)
{
// pulisco i diversi vettori
m_vPnt.clear() ;
m_vPar.clear() ;
m_vPrevDer.clear() ;
m_vNextDer.clear() ;
m_vSplits.clear() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveByApprox::AddPoint( const Point3d& ptP)
{
// se il punto coincide con il precedente, lo salto
if ( ! m_vPnt.empty() && AreSamePointApprox( ptP, m_vPnt.back()))
return false ;
// aggiungo il punto
try { m_vPnt.push_back( ptP) ; }
catch ( ...) { return false ; }
return true ;
}
//----------------------------------------------------------------------------
ICurve*
CurveByApprox::GetCurve( int nType, double dLinTol, double dAngTolDeg, double dLinFea)
{
// calcolo le tangenti
if ( ! CalcAkimaTangents( true))
return nullptr ;
// divido in tratti (rettilinei o tra spigoli)
if ( ! CalcSplitPoints( dLinTol, dAngTolDeg, dLinFea))
return nullptr ;
// approssimo ogni singolo tratto
PolyArc PA ;
int nPrev = 0 ;
int nSplits = int( m_vSplits.size()) ;
for ( int i = 0 ; i < nSplits ; ++ i) {
// creo la polilinea che unisce i punti
PolyLine PL ;
for ( int j = nPrev ; j <= m_vSplits[i] ; ++ j)
PL.AddUPoint( j, m_vPnt[j]) ;
// eseguo l'approssimazione con archi
if ( ! BiArcOrSplit( 0, PL, dLinTol, dAngTolDeg, PA))
return nullptr ;
// salvo fine come prox inizio
nPrev = m_vSplits[i] ;
}
// creo la composita formata da questa approssimazione
PtrOwner<CurveComposite> pCC( CreateBasicCurveComposite()) ;
if ( ! pCC->FromPolyArc( PA))
return nullptr ;
// eventuale fusione di curve compatibili
pCC->MergeCurves( dLinTol, dAngTolDeg) ;
return Release( pCC) ;
}
//----------------------------------------------------------------------------
bool
CurveByApprox::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 ;
}
//----------------------------------------------------------------------------
bool
CurveByApprox::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 ;
}
//----------------------------------------------------------------------------
bool
CurveByApprox::CalcSplitPoints( double dLinTol, double dAngTolDeg, double dLinFea)
{
// precalcolo funzioni dei limiti
double dSqLinTol = dLinTol * dLinTol ;
double dAngTolCos = cos( max( dAngTolDeg, 0.) * DEGTORAD) ;
double dSqLinFea = dLinFea * dLinFea ;
// cerco punti angolosi e punti di separazione di tratti lineari e non
// verifico i punti tranne primo e ultimo
int nSize = int( m_vPnt.size()) ;
for ( int i = 1 ; i < nSize - 1 ; ++ i) {
// se angolo, allora punto di divisione
if ( ( m_vPrevDer[i] * m_vNextDer[i]) < dAngTolCos) {
m_vSplits.push_back(i) ;
continue ;
}
// verifico linearità del tratto precedente
Vector3d vtPrev = m_vPnt[i] - m_vPnt[i-1] ;
bool bPrevLin = vtPrev.SqLen() > dSqLinFea &&
( m_vNextDer[i-1] * m_vPrevDer[i]) > dAngTolCos &&
( vtPrev ^ m_vNextDer[i-1]).SqLen() < dSqLinTol &&
( vtPrev ^ m_vPrevDer[i]).SqLen() < dSqLinTol ;
// verifico linearità del tratto successivo
Vector3d vtNext = m_vPnt[i+1] - m_vPnt[i] ;
bool bNextLin = vtNext.SqLen() > dSqLinFea &&
( m_vNextDer[i] * m_vPrevDer[i+1]) > dAngTolCos &&
( vtNext ^ m_vNextDer[i]).SqLen() < dSqLinTol &&
( vtNext ^ m_vPrevDer[i+1]).SqLen() < dSqLinTol ;
// se uno lineare e l'altro no, allora punto di divisione
if ( ( bPrevLin && ! bNextLin) || ( ! bPrevLin && bNextLin)) {
m_vSplits.push_back(i) ;
continue ;
}
// se entrambi non lineari, vado oltre
if ( ! bPrevLin && ! bNextLin)
continue ;
// sono entrambi tratti lineari, verifico siano allineati
DistPointLine dPL( m_vPnt[i], m_vPnt[i-1], m_vPnt[i+1]) ;
double dSqDist ;
if ( ! dPL.GetSqDist( dSqDist) || dSqDist > dSqLinTol) {
m_vSplits.push_back(i) ;
continue ;
}
}
// in ogni caso ci deve essere l'ultimo punto
m_vSplits.push_back(nSize - 1) ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveByApprox::BiArcOrSplit( int nLev, PolyLine& PL, double dLinTol, double dAngTolDeg, PolyArc& PA) const
{
const int MAX_LEV = 10 ;
// se non chiusa
if ( ! PL.IsClosed()) {
PtrOwner<ICurve> pCrv ;
double dMaxDist ;
// se la polilinea ha più di 2 punti
if ( PL.GetPointNbr() > 2) {
// calcolo punti e direzioni agli estremi della polilinea usando la curva di Bezier
int nI ;
double dU ;
Point3d ptP0, ptP1 ;
double dDir0Deg, dDir1Deg ;
if ( ! PL.GetFirstU( dU))
return false ;
nI = int( dU) ;
ptP0 = m_vPnt[nI] ;
m_vNextDer[nI].ToSpherical( nullptr, nullptr, &dDir0Deg) ;
if ( ! PL.GetLastU( dU))
return false ;
nI = int( dU) ;
ptP1 = m_vPnt[nI] ;
m_vPrevDer[nI].ToSpherical( nullptr, nullptr, &dDir1Deg) ;
// costruisco un biarco sulla polilinea (secondo metodo di Z. Sir)
pCrv.Set( GetBiArc( ptP0, dDir0Deg, ptP1, dDir1Deg, PL, dMaxDist)) ;
if ( IsNull( pCrv))
return false ;
}
// se la polilinea è formata da 2 punti
else if ( PL.GetPointNbr() == 2) {
// se molto vicini, esco
double dLen ;
if ( PL.GetLength( dLen) && dLen < EPS_SMALL)
return true ;
// costruisco la retta che li unisce
PtrOwner<CurveComposite> pCC( CreateBasicCurveComposite()) ;
if ( IsNull( pCC))
return false ;
if ( ! pCC->FromPolyLine( PL))
return false ;
pCrv.Set( Release( pCC)) ;
dMaxDist = 0 ;
}
// se la polilinea ha un solo punto, esco
else
return true ;
// se raggiunto il massimo livello di recursione, forzo l'accettazione del biarco
if ( nLev >= MAX_LEV) {
dMaxDist = 0 ;
// segnalo situazione per debug
LOG_DBG_ERR( GetEGkLogger(), "ERROR : Exceeded recursions")
}
// se lunghezza abbastanza picccola, forzo l'accettazione della curva
double dLen ;
if ( PL.GetApproxLength( dLen) && dLen < 10 * EPS_SMALL)
dMaxDist = 0 ;
// se in tolleranza
if ( dMaxDist <= dLinTol) {
// creo un nuovo poliarco
PolyArc PA2 ;
if ( ! pCrv->ApproxWithArcs( dLinTol, dAngTolDeg, PA2))
return false ;
// aggiusto la parametrizzazione
double dUStart, dUEnd ;
PL.GetFirstU( dUStart) ;
PL.GetLastU( dUEnd) ;
PA2.ParamLinearTransform( dUStart, dUEnd) ;
// accodo all'esistente
if ( ! PA.Join( PA2))
return false ;
// esco
return true ;
}
}
// se raggiunto il massimo livello di recursione, errore
if ( nLev >= MAX_LEV) {
// segnalo situazione per debug
LOG_DBG_ERR( GetEGkLogger(), "ERROR : Exceeded recursions")
return false ;
}
// spezzo l'intervallo in due parti a metà
double dParStart, dParEnd ;
if ( ! PL.GetFirstU( dParStart) || ! PL.GetLastU( dParEnd))
return false ;
double dParMid = 0.5 * ( dParStart + dParEnd) ;
PolyLine PL2 ;
if ( ! PL.Split( dParMid, PL2))
return false ;
// prima metà
if ( ! BiArcOrSplit( nLev + 1, PL, dLinTol, dAngTolDeg, PA))
return false ;
// seconda metà
return BiArcOrSplit( nLev + 1, PL2, dLinTol, dAngTolDeg, PA) ;
}
Reference in New Issue View Git Blame Copy Permalink