Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
ef6300a3bc8a399bed221ec3fa0345e1535482c9
EgtGeomKernel/CurveBezier.cpp
T
Dario Sassi 2cbdfc1db9 EgtGeomKernel : 
- cambiato Load, con uso di NgeReader.
2014-04-15 06:48:39 +00:00

1412 lines
41 KiB
C++
Raw Blame History

//----------------------------------------------------------------------------
// EgalTech 2013-2014
//----------------------------------------------------------------------------
// File : CurveBezier.cpp Data : 26.03.14 Versione : 1.5c9
// Contenuto : Implementazione della classe Curva di Bezier.
//
//
//
// Modifiche : 28.12.12 DS Creazione modulo.
// 26.03.14 DS Modif. Save/Load, ora ogni punto su una linea.
//
//----------------------------------------------------------------------------
//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CurveBezier.h"
#include "GeoConst.h"
#include "DistPointLine.h"
#include "GeoObjFactory.h"
#include "NgeWriter.h"
#include "NgeReader.h"
#include "PolynomialPoint3d.h"
#include "/EgtDev/Include/EGkCurveArc.h"
#include "/EgtDev/Include/ENkPolynomial.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include <new>
using namespace std ;
//----------------------------------------------------------------------------
GEOOBJ_REGISTER( CRV_BEZ, NGE_C_BEZ, CurveBezier) ;
//----------------------------------------------------------------------------
CurveBezier::CurveBezier( void)
: m_nStatus( TO_VERIFY), m_nDeg( 0), m_bRat( false), m_dParSing( -2),
m_nDimArr( 0), m_aPtCtrl( nullptr), m_aWeCtrl( nullptr)
{
}
//----------------------------------------------------------------------------
CurveBezier::~CurveBezier( void)
{
if ( m_nDimArr > 0) {
delete [] m_aPtCtrl ;
delete [] m_aWeCtrl ;
}
m_nDimArr = 0 ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Init( int nDeg, bool bIsRational)
{
// verifico validità grado
if ( nDeg < 1 || nDeg > MAXDEG)
return false ;
// se spazio dinamico sufficiente
if ( ( nDeg + 1) <= m_nDimArr)
;
// se altrimenti spazio statico sufficiente
else if ( ( nDeg + 1) <= ST_PTC) {
m_aPtCtrl = m_aStPtCtrl ;
if ( bIsRational)
m_aWeCtrl = m_aStWeCtrl ;
}
// altrimenti
else {
// se necessaria, pulizia
if ( m_nDimArr > 0) {
delete [] m_aPtCtrl ;
delete [] m_aWeCtrl ;
}
// alloco punti
m_aPtCtrl = new Point3d [ nDeg + 1] ;
if ( m_aPtCtrl == nullptr)
return false ;
// se razionale, alloco pesi
if ( bIsRational) {
m_aWeCtrl = new double [ nDeg + 1] ;
if ( m_aWeCtrl == nullptr) {
delete [] m_aPtCtrl ;
return false ;
}
}
// salvo dimensione array allocati
m_nDimArr = nDeg + 1 ;
}
// salvo il grado
m_nDeg = nDeg ;
// 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) {
for ( int i = 0 ; i <= m_nDeg ; ++ i)
m_aWeCtrl[i] = 1 ;
}
m_nStatus = TO_VERIFY ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return Validate() ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::SetControlPoint( int nInd, const Point3d& ptCtrl)
{
// verifico validità indice
if ( m_nStatus != OK || m_bRat || nInd < 0 || nInd > m_nDeg)
return false ;
// assegno il valore
m_aPtCtrl[nInd] = ptCtrl ;
// se razionale, metto il peso a 1
if ( m_bRat)
m_aWeCtrl[nInd] = 1 ;
// annullo analisi presenza singolarità
m_dParSing = - 2 ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::SetControlPoint( int nInd, const Point3d& ptCtrl, double dW)
{
// verifico validità, razionalità e indice
if ( m_nStatus != OK || ! m_bRat || nInd < 0 || nInd > m_nDeg)
return false ;
// verifico che il peso non sia nullo o negativo
if ( dW < EPS_SMALL)
return false ;
// assegno il valore e il peso
m_aPtCtrl[nInd] = ptCtrl ;
m_aWeCtrl[nInd] = dW ;
// annullo analisi presenza singolarità
m_dParSing = - 2 ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::FromArc( const ICurveArc& crArc)
{
const double MAX_ANG_CEN_DEG = 120 ;
double dAngCen ;
double dCosAhalf ;
double dW ;
Point3d ptStart ;
Point3d ptEnd ;
Point3d ptEndNoDN ;
Point3d ptMed ;
Point3d ptNew ;
Point3d ptNew1 ;
Point3d ptNew2 ;
Vector3d vtDeltaN ;
Vector3d vtDir ;
if ( ! crArc.IsValid())
return false ;
dAngCen = crArc.GetAngCenter() ;
if ( fabs( dAngCen) > MAX_ANG_CEN_DEG)
return false ;
dCosAhalf = cos( 0.5 * dAngCen * DEGTORAD) ;
// se arco piatto, basta quadrica razionale
// algoritmo di Sederberg (BYU) CAGD cap. 2 pag. 33
if ( fabs( crArc.GetDeltaN()) < EPS_ZERO) {
// peso del punto di controllo intermedio
dW = dCosAhalf ;
// calcolo dei punti di controllo
crArc.GetStartPoint( ptStart) ;
crArc.GetEndPoint( ptEnd) ;
ptMed = Media( ptStart, ptEnd, 0.5) ;
vtDir = ptMed - crArc.GetCenter() ;
ptNew = crArc.GetCenter() + vtDir / ( dCosAhalf * dCosAhalf) ;
// inserimento nella curva dei punti di controllo con i pesi
Init( 2, true) ;
SetControlPoint( 0, ptStart, 1) ;
SetControlPoint( 1, ptNew, dW) ;
SetControlPoint( 2, ptEnd, 1) ;
}
// altrimenti per elica serve cubica razionale
// algoritmo di Juhasz HelixApproxByCubicBezier CAD 1995
else {
// peso dei punti di controllo intermedi
dW = ( 1 + 2 * dCosAhalf) / 3 ;
// calcolo dei punti di controllo
crArc.GetStartPoint( ptStart) ;
crArc.GetEndPoint( ptEnd) ;
vtDeltaN = crArc.GetDeltaN() * crArc.GetNormVersor() ;
ptEndNoDN = ptEnd - vtDeltaN ;
ptMed = Media( ptStart, ptEndNoDN, 0.5) ;
vtDir = ptMed - crArc.GetCenter() ;
ptNew = crArc.GetCenter() + vtDir / ( dCosAhalf * dCosAhalf) ;
ptNew1 = ( ptStart + 2 * dCosAhalf * ptNew) / ( 1 + 2 * dCosAhalf) +
0.5 * ( 1 - dCosAhalf / ( 2 + dCosAhalf)) * vtDeltaN ;
ptNew2 = ( ptEndNoDN + 2 * dCosAhalf * ptNew) / ( 1 + 2 * dCosAhalf) +
0.5 * ( 1 + ( dCosAhalf / ( 2 + dCosAhalf))) * vtDeltaN ;
// inserimento nella curva dei punti di controllo con i pesi
Init( 3, true) ;
SetControlPoint( 0, ptStart, 1) ;
SetControlPoint( 1, ptNew1, dW) ;
SetControlPoint( 2, ptNew2, dW) ;
SetControlPoint( 3, ptEnd, 1) ;
}
return true ;
}
//----------------------------------------------------------------------------
const Point3d&
CurveBezier::GetControlPoint( int nInd, bool* pbOk) const
{
// verifico validità e indice
if ( m_nStatus != OK || nInd < 0 || nInd > m_nDeg) {
if ( pbOk != NULL)
*pbOk = false ;
return ORIG ;
}
// ritorno i dati
if ( pbOk != NULL)
*pbOk = true ;
return m_aPtCtrl[nInd] ;
}
//----------------------------------------------------------------------------
double
CurveBezier::GetControlWeight( int nInd, bool* pbOk) const
{
// verifico validità, razionalità e indice
if ( m_nStatus != OK || ! m_bRat || nInd < 0 || nInd > m_nDeg) {
if ( pbOk != NULL)
*pbOk = false ;
return 0 ;
}
// ritorno i dati
if ( pbOk != NULL)
*pbOk = true ;
return m_aWeCtrl[nInd] ;
}
//----------------------------------------------------------------------------
CurveBezier*
CurveBezier::Clone( void) const
{
// alloco oggetto
CurveBezier* pCrv = new(nothrow) CurveBezier ;
if ( pCrv != nullptr) {
if ( ! pCrv->Copy( *this)) {
delete pCrv ;
return nullptr ;
}
}
return pCrv ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Copy( const IGeoObj* pGObjSrc)
{
const CurveBezier* pCB = dynamic_cast<const CurveBezier*>( pGObjSrc) ;
if ( pCB == nullptr)
return false ;
return Copy( *pCB) ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Copy( const CurveBezier& cbSrc)
{
if ( &cbSrc == this)
return true ;
if ( ! Init( cbSrc.m_nDeg, cbSrc.m_bRat))
return false ;
for ( int i = 0 ; i <= m_nDeg ; ++ i) {
m_aPtCtrl[i] = cbSrc.m_aPtCtrl[i] ;
if ( cbSrc.m_bRat)
m_aWeCtrl[i] = cbSrc.m_aWeCtrl[i] ;
}
return true ;
}
//----------------------------------------------------------------------------
GeoObjType
CurveBezier::GetType( void) const
{
return static_cast<GeoObjType>( GEOOBJ_GETTYPE( CurveBezier)) ;
}
//----------------------------------------------------------------------------
const string&
CurveBezier::GetTitle( void) const
{
static const string sTitle = "CBezier" ;
return sTitle ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Dump( string& sOut, const char* szNewLine) const
{
// flag razionale
sOut += ( m_bRat ? "Rat " : "Int ") ;
// grado
sOut += "Deg=" + ToString( m_nDeg) + szNewLine ;
// ciclo sui punti di controllo ( con pesi se razionale)
for ( int i = 0 ; i <= m_nDeg ; ++ i) {
sOut += "PC(" + ToString( m_aPtCtrl[i]) ;
if ( m_bRat)
sOut += "," + ToString( m_aWeCtrl[i]) ;
sOut += string( ")") + szNewLine ;
}
return true ;
}
//----------------------------------------------------------------------------
int
CurveBezier::GetNgeId( void) const
{
return GEOOBJ_GETNGEID( CurveBezier) ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Save( NgeWriter& ngeOut) const
{
// flag razionale
ngeOut.WriteBool( m_bRat, ";") ;
// grado
ngeOut.WriteInt( m_nDeg, ";", true) ;
// ciclo sui punti di controllo ( con pesi se razionale)
for ( int i = 0 ; i <= m_nDeg ; ++ i) {
if ( ! m_bRat)
ngeOut.WritePoint( m_aPtCtrl[i], ";", true) ;
else
ngeOut.WritePointW( m_aPtCtrl[i], m_aWeCtrl[i], ";", true) ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Load( NgeReader& ngeIn)
{
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
// leggo la prossima linea ( 2 parametri)
// recupero il flag razionale
bool bIsRat ;
if ( ! ngeIn.ReadBool( bIsRat, ";"))
return false ;
// recupero il grado
int nDeg ;
if ( ! ngeIn.ReadInt( nDeg, ";", true))
return false ;
// inizializzo la curva di Bezier
if ( ! Init( nDeg, bIsRat))
return false ;
// se integrale
if ( ! bIsRat) {
// recupero e setto punti di controllo
string sLine ;
Point3d ptP ;
for ( int i = 0 ; i <= nDeg ; ++ i) {
// leggo la prossima linea ( un punto)
if ( ! ngeIn.ReadPoint( ptP, ";", true))
return false ;
// lo assegno
if ( ! SetControlPoint( i, ptP))
return false ;
}
}
// altrimenti razionale
else {
// recupero e setto punti di controllo
string sLine ;
Point3d ptP ;
double dW ;
for ( int i = 0 ; i <= nDeg ; ++ i) {
// leggo la prossima linea ( un punto con peso)
if ( ! ngeIn.ReadPointW( ptP, dW, ";", true))
return false ;
// lo assegno
if ( ! SetControlPoint( i, ptP, dW))
return false ;
}
}
// eseguo validazione
return Validate() ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Validate( void)
{
if ( m_nStatus == TO_VERIFY)
m_nStatus = ( ( m_nDeg > 0 && m_aPtCtrl != nullptr) ? OK : ERR) ;
return ( m_nStatus == OK) ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetLocalBBox( BBox3d& b3Loc) const
{
// verifico lo stato
if ( m_nStatus != OK)
return false ;
// assegno il box in locale
b3Loc.Reset() ;
for ( int i = 0 ; i <= m_nDeg ; ++ i)
b3Loc.Add( m_aPtCtrl[i]) ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetBBox( const Frame3d& frRef, BBox3d& b3Ref) const
{
// verifico lo stato
if ( m_nStatus != OK)
return false ;
// verifico validità del frame
if ( frRef.GetType() == Frame3d::ERR)
return false ;
// assegno il box nel riferimento
b3Ref.Reset() ;
for ( int i = 0 ; i <= m_nDeg ; ++ i) {
Point3d ptTemp = m_aPtCtrl[i] ;
ptTemp.ToGlob( frRef) ;
b3Ref.Add( ptTemp) ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetStartPoint( Point3d& ptStart) const
{
// verifico lo stato
if ( m_nStatus != OK)
return false ;
// assegno il punto
ptStart = m_aPtCtrl[0] ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetEndPoint( Point3d& ptEnd) const
{
// verifico lo stato
if ( m_nStatus != OK)
return false ;
// assegno il punto
ptEnd = m_aPtCtrl[m_nDeg] ;
return true ;
}
//----------------------------------------------------------------------------
// Calcolo i polinomi di Bernstein e da questi il punto ( vedi Piegl-Tiller)
//----------------------------------------------------------------------------
bool
CurveBezier::GetPointD1D2( double dU, Side nS, Point3d& ptPos, Vector3d* pvtDer1, Vector3d* pvtDer2) const
{
int i ;
double dBern[MAXDEG] ;
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// il parametro U deve essere compreso tra 0 e 1
if ( dU < 0)
dU = 0 ;
else if ( dU > 1)
dU = 1 ;
// se forma polinomiale (o integrale)
if ( ! m_bRat) {
// calcolo dei polinomi di Bernstein di grado opportuno
dBern[0] = 1 ;
for ( i = 1 ; i <= m_nDeg - 2 ; ++ i)
IncreaseBernsteinOneDegree( dU, i, dBern) ;
// se richiesto, calcolo della derivata seconda
if ( pvtDer2 != nullptr && pvtDer1 != nullptr) {
pvtDer2->Set( 0, 0, 0) ;
for ( i = 0 ; i <= m_nDeg - 2 ; ++ i) {
*pvtDer2 += dBern[i] * ( m_aPtCtrl[i+2] + m_aPtCtrl[i] - 2 * m_aPtCtrl[i+1]) ;
}
*pvtDer2 *= m_nDeg * ( m_nDeg - 1) ;
}
// aumento il grado
IncreaseBernsteinOneDegree( dU, m_nDeg - 1, dBern) ;
// se richiesto, calcolo della derivata prima
if ( pvtDer1 != nullptr) {
pvtDer1->Set( 0, 0, 0) ;
for ( i = 0 ; i <= m_nDeg - 1 ; ++ i) {
*pvtDer1 += dBern[i] * ( m_aPtCtrl[i+1] - m_aPtCtrl[i]) ;
}
*pvtDer1 *= m_nDeg ;
}
// aumento il grado
IncreaseBernsteinOneDegree( dU, m_nDeg, dBern) ;
// calcolo del punto
ptPos.Set( 0, 0, 0) ;
for ( i = 0 ; i <= m_nDeg ; ++ i)
ptPos += dBern[i] * m_aPtCtrl[i] ;
}
// altrimenti forma razionale
else {
double dW ;
double dW1 ;
double dW2 ;
double dInvW ;
Point3d aPtWCtrl[MAXDEG] ;
Vector3d vtPos ;
// porto i punti in forma omogenea moltiplicandoli per i pesi
for ( i = 0 ; i <= m_nDeg ; ++ i)
aPtWCtrl[i] = m_aWeCtrl[i] * m_aPtCtrl[i] ;
// calcolo dei polinomi di Bernstein di grado opportuno
dBern[0] = 1 ;
for ( i = 1 ; i <= m_nDeg - 2 ; ++ i)
IncreaseBernsteinOneDegree( dU, i, dBern) ;
// se richiesto, calcolo della derivata seconda
if ( pvtDer2 != nullptr && pvtDer1 != nullptr) {
pvtDer2->Set( 0, 0, 0) ; dW2 = 0 ;
for ( i = 0 ; i <= m_nDeg - 2 ; ++ i) {
*pvtDer2 += dBern[i] * ( aPtWCtrl[i+2] + aPtWCtrl[i] - 2 * aPtWCtrl[i+1]) ;
dW2 += dBern[i] * ( m_aWeCtrl[i+2] + m_aWeCtrl[i] - 2 * m_aWeCtrl[i+1]) ;
}
*pvtDer2 *= m_nDeg * ( m_nDeg - 1) ;
dW2 *= m_nDeg * ( m_nDeg - 1) ;
}
// aumento il grado
IncreaseBernsteinOneDegree( dU, m_nDeg - 1, dBern) ;
// se richiesto, calcolo della derivata prima
if ( pvtDer1 != nullptr) {
pvtDer1->Set( 0, 0, 0) ; dW1 = 0 ;
for ( i = 0 ; i <= m_nDeg - 1 ; ++ i) {
*pvtDer1 += dBern[i] * ( aPtWCtrl[i+1] - aPtWCtrl[i]) ;
dW1 += dBern[i] * ( m_aWeCtrl[i+1] - m_aWeCtrl[i]) ;
}
*pvtDer1 *= m_nDeg ;
dW1 *= m_nDeg ;
}
// aumento il grado
IncreaseBernsteinOneDegree( dU, m_nDeg, dBern) ;
// calcolo del punto
ptPos.Set( 0, 0, 0) ; dW = 0 ;
for ( i = 0 ; i <= m_nDeg ; ++ i) {
ptPos += dBern[i] * aPtWCtrl[i] ;
dW += dBern[i] * m_aWeCtrl[i] ;
}
// ritrasformo da forma omogenea a forma standard
dInvW = 1 / ( ( dW > EPS_ZERO) ? dW : EPS_ZERO) ;
ptPos *= dInvW ;
if ( pvtDer1 != nullptr) {
vtPos.Set( ptPos.x, ptPos.y, ptPos.z) ;
*pvtDer1 = ( *pvtDer1 - dW1 * vtPos) * dInvW ;
if ( pvtDer2 != nullptr)
*pvtDer2 = ( *pvtDer2 - 2 * dW1 * *pvtDer1 - dW2 * vtPos) * dInvW ;
}
}
return true ;
}
//----------------------------------------------------------------------------
void
CurveBezier::IncreaseBernsteinOneDegree( double dU, int nDeg, double dBern[]) const
{
int j ;
double dU1 ;
double dTot ;
double dTmp ;
dU1 = 1 - dU ;
dTot = 0 ;
for ( j = 0 ; j < nDeg ; ++ j) {
dTmp = dBern[j] ;
dBern[j] = dTot + dU1 * dTmp ;
dTot = dU * dTmp ;
}
dBern[nDeg] = ( nDeg > 0 ? dTot : 1) ;
}
//----------------------------------------------------------------------------
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
{
// la curva deve essere valida
if ( m_nStatus != OK)
return false ;
// ciclo sui punti di controllo
dLen = 0 ;
for ( int i = 1 ; i <= m_nDeg ; ++ i)
dLen += Dist( m_aPtCtrl[i], m_aPtCtrl[i-1]) ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetLength( double& dLen) const
{
return ( GetLengthAtParam( 1, dLen) && dLen > EPS_SMALL) ;
}
//----------------------------------------------------------------------------
double
CurveBezier::GetSegmentLength( int nLev, double dU0, double dU1, double dU2,
const Point3d& ptP0, const Point3d& ptP1, const Point3d& ptP2) const
{
const int MAX_LEV = 10 ;
const double MAX_ARC = 1.05 ;
const double LEN_RATIO = 1.2 ;
double dD1 ;
double dDa ;
double dDb ;
double dD2 ;
// algoritmo dell'approssimazione con arco per tre punti di S. Vincent e D. Forsey
// verifica superamento del massimo livello di recursione per debug
if ( nLev >= MAX_LEV)
LOG_DBG_ERR( GetEGkLogger(), "ERROR : Exceeded recursions")
// calcolo delle distanze
dD1 = Dist( ptP0, ptP2) ;
dDa = Dist( ptP0, ptP1) ;
dDb = Dist( ptP1, ptP2) ;
dD2 = dDa + dDb ;
// distanza piccola o massimo livello di recursione, approx arco ok
if ( dD2 < EPS_SMALL || nLev >= MAX_LEV)
return ( dD2 + ( dD2 - dD1) / 3) ;
// se ci sono anomalie si suddivide
if ( dD1 < EPS_SMALL || dD2 / dD1 > MAX_ARC ||
dDa < 0.1 * EPS_SMALL || dDb / dDa > LEN_RATIO ||
dDb < 0.1 * EPS_SMALL || dDa / dDb > LEN_RATIO) {
double dUMid ;
Point3d ptMid ;
// prima metà
dUMid = 0.5 * ( dU0 + dU1) ;
GetPointD1D2( dUMid, ICurve::FROM_MINUS, ptMid) ;
double dD3 = GetSegmentLength( nLev +1, dU0, dUMid, dU1,
ptP0, ptMid, ptP1) ;
// seconda metà
dUMid = 0.5 * ( dU1 + dU2) ;
GetPointD1D2( dUMid, ICurve::FROM_MINUS, ptMid) ;
double dD4 = GetSegmentLength( nLev + 1, dU1, dUMid, dU2,
ptP1, ptMid, ptP2) ;
// la lunghezza è la somma delle due
return ( dD3 + dD4) ;
}
// caso normale
return ( dD2 + ( dD2 - dD1) / 3) ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetLengthAtParam( double dU, double& dLen) const
{
const int NUM_SEG = 16 ;
double dUIni ;
double dUMid ;
double dUFin ;
Point3d ptIni ;
Point3d ptMid ;
Point3d ptFin ;
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// fuori dominio del parametro -> errore
if ( dU < 0 - EPS_ZERO)
return false ;
if ( dU > ( 1 + EPS_ZERO))
return false ;
// inizio
if ( dU < 0 + EPS_ZERO) {
dLen = 0 ;
return true ;
}
// parto con un numero fisso di punti e poi affino dove serve
dLen = 0 ;
dUIni = 0 ;
GetPointD1D2( 0, ICurve::FROM_MINUS, ptIni) ;
for ( int i = 1 ; i <= NUM_SEG ; ++ i) {
// ricavo il punto
dUFin = i / (double) NUM_SEG * dU ;
GetPointD1D2( dUFin, ICurve::FROM_MINUS, ptFin) ;
// se punto dispari
if ( ( i % 2) == 1) {
// assegno nuovo medio
dUMid = dUFin ;
ptMid = ptFin ;
}
// se punto pari, calcolo la lunghezza
else {
dLen += GetSegmentLength( 0, dUIni, dUMid, dUFin, ptIni, ptMid, ptFin) ;
// nuovo iniziale prende i valori del finale
dUIni = dUFin ;
ptIni = ptFin ;
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetSegmentParam( double dLen, double& dCurrLen, double& dSegLen,
double& dUIni, double& dUFin) const
{
const int NUM_SEG = 16 ;
double dUMid ;
Point3d ptIni ;
Point3d ptMid ;
Point3d ptFin ;
double dUStart = dUIni ;
double dUEnd = dUFin ;
GetPointD1D2( dUIni, ICurve::FROM_MINUS, ptIni) ;
for ( int i = 1 ; i <= NUM_SEG ; ++ i) {
// ricavo il punto
dUFin = dUStart + i / (double) NUM_SEG * ( dUEnd - dUStart) ;
GetPointD1D2( dUFin, ICurve::FROM_MINUS, ptFin) ;
// se punto dispari
if ( ( i % 2) == 1) {
// assegno nuovo medio
dUMid = dUFin ;
ptMid = ptFin ;
}
// se punto pari, calcolo la lunghezza
else {
dSegLen = GetSegmentLength( 0, dUIni, dUMid, dUFin, ptIni, ptMid, ptFin) ;
if ( ( dCurrLen + dSegLen) >= dLen) {
return true ;
}
else
dCurrLen += dSegLen ;
// nuovo iniziale prende i valori del finale
dUIni = dUFin ;
ptIni = ptFin ;
}
}
return false ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::GetParamAtLength( double dLen, double& dU) const
{
const double MAX_SEG_LEN = 0.2 ;
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// inizio
if ( dLen < EPS_SMALL) {
dU = 0 ;
return true ;
}
// approfondisco nell'intervallo di interesse
bool bOk ;
double dSegLen ;
double dCurrLen = 0 ;
double dUIni = 0 ;
double dUFin = 1 ;
while ( ( bOk = GetSegmentParam( dLen, dCurrLen, dSegLen, dUIni, dUFin)) && dSegLen > MAX_SEG_LEN)
;
// aggiustamento finale parametro
if ( bOk)
dU = dUIni + ( dLen - dCurrLen) / dSegLen * ( dUFin - dUIni) ;
return bOk ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::FlatOrSplit( int nLev, const CurveBezier& crvBez, double dParStart, double dParEnd,
double dLinTol, double dAngTolDeg, PolyLine& PL) const
{
const int MAX_LEV = 10 ;
int i ;
double dMaxSqDist ;
double dSqDist ;
double dAngDeg ;
double dPolLen ;
double dLen ;
Point3d ptLast ;
Vector3d vtDirI ;
Vector3d vtDirF ;
Vector3d vtDir ;
// se raggiunto il massimo livello di recursione ...
if ( nLev >= MAX_LEV) {
// segnalo situazione per debug
LOG_DBG_ERR( GetEGkLogger(), "ERROR : Exceeded recursions")
// considero la curva piatta (inserisco il punto se abbastanza lontano dal precedente) ed esco
PL.GetLastPoint( ptLast) ;
if ( SqDist( ptLast, crvBez.m_aPtCtrl[m_nDeg]) > ( dLinTol * dLinTol))
PL.AddUPoint( dParEnd, crvBez.m_aPtCtrl[m_nDeg]) ;
return true ;
}
// massima distanza al quadrato dei punti di controllo intermedi dalla linea tra primo e ultimo
dMaxSqDist = 0 ;
DirDist( crvBez.m_aPtCtrl[0], crvBez.m_aPtCtrl[m_nDeg], vtDir, dLen) ;
for ( i = 1 ; i < m_nDeg ; i ++) {
DistPointLine dstPL( crvBez.m_aPtCtrl[i], crvBez.m_aPtCtrl[0], vtDir, dLen) ;
if ( dstPL.GetSqDist( dSqDist) && dSqDist > dMaxSqDist)
dMaxSqDist = dSqDist ;
}
// se distanza entro tolleranza
if ( dMaxSqDist <= ( dLinTol * dLinTol)) {
// deviazione angolare tra primo e ultimo tratto del poligono di controllo (grado >= 1)
vtDirI = crvBez.m_aPtCtrl[1] - crvBez.m_aPtCtrl[0] ;
vtDirF = crvBez.m_aPtCtrl[m_nDeg] - crvBez.m_aPtCtrl[m_nDeg-1] ;
vtDirI.GetAngle( vtDirF, dAngDeg) ;
// se deviazione angolare entro tolleranza oppure lunghezza poligono controllo entro tolleranza
if ( fabs( dAngDeg) <= dAngTolDeg ||
( crvBez.GetControlPolygonLength( dPolLen) && dPolLen <= dLinTol)) {
// considero la curva piatta (inserisco il punto se abbastanza lontano dal precedente) ed esco
PL.GetLastPoint( ptLast) ;
if ( SqDist( ptLast, crvBez.m_aPtCtrl[m_nDeg]) > ( dLinTol * dLinTol))
PL.AddUPoint( dParEnd, crvBez.m_aPtCtrl[m_nDeg]) ;
return true ;
}
}
// curva da dividere
{
double dParDiv ;
double dParMid ;
CurveBezier crvBez1 ;
// se prima suddivisione e c'è singolarità, divido su questa
if ( nLev == 0 && GetSingularParam( dParDiv) > 0)
;
// altrimenti divido a metà
else
dParDiv = 0.5 ;
dParMid = dParDiv * ( dParStart + dParEnd) ;
// prima metà
crvBez1 = crvBez ;
crvBez1.TrimEndAtParam( dParDiv) ;
if ( ! FlatOrSplit( nLev + 1, crvBez1, dParStart, dParMid, dLinTol, dAngTolDeg, PL))
return false ;
// seconda metà
crvBez1 = crvBez ;
crvBez1.TrimStartAtParam( dParDiv) ;
if ( ! FlatOrSplit( nLev + 1, crvBez1, dParMid, dParEnd, dLinTol, dAngTolDeg, PL))
return false ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::ApproxWithLines( double dLinTol, double dAngTolDeg, PolyLine& PL) const
{
// pulisco la polilinea
PL.Clear() ;
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// se di primo grado, basta inserire gli estremi
if ( m_nDeg == 1) {
PL.AddUPoint( 0, m_aPtCtrl[0]) ;
PL.AddUPoint( 1, m_aPtCtrl[m_nDeg]) ;
return true ;
}
// limiti minimi su tolleranza e deviazione angolare
dLinTol = max( dLinTol, LIN_TOL_MIN) ;
dAngTolDeg = max( dAngTolDeg, ANG_TOL_MIN_DEG) ;
// inserisco il punto iniziale
PL.AddUPoint( 0, m_aPtCtrl[0]) ;
// verifico se va divisa
FlatOrSplit( 0, *this, 0, 1, dLinTol, dAngTolDeg, PL) ;
// se è stato inserito un solo punto, aggiungo il finale
if ( PL.GetPointNbr() == 1)
PL.AddUPoint( 1, m_aPtCtrl[m_nDeg]) ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Invert( void)
{
// la curva deve essere valida
if ( m_nStatus != OK)
return false ;
// inverto i punti di controllo
int nMid = ( m_nDeg + 1) / 2 ;
for ( int i = 0 ; i < nMid ; ++ i) {
Point3d ptTemp = m_aPtCtrl[i] ;
m_aPtCtrl[i] = m_aPtCtrl[m_nDeg-i] ;
m_aPtCtrl[m_nDeg-i] = ptTemp ;
}
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::TrimStartAtParam( double dUTrim)
{
// la curva deve essere valida
if ( m_nStatus != OK)
return false ;
// se devo togliere dall'inizio o prima non tolgo niente
if ( dUTrim <= 0)
return true ;
// se devo togliere a partire dalla fine o dopo non rimane niente
if ( dUTrim >= 1)
return false ;
// eseguo il trim (algoritmo di de Casteljau)
for ( int k = 1 ; k <= m_nDeg ; k ++) {
for ( int i = 0 ; i <= m_nDeg - k ; i ++) {
if ( ! m_bRat) {
m_aPtCtrl[i] = ( 1 - dUTrim) * m_aPtCtrl[i] + dUTrim * m_aPtCtrl[i+1] ;
}
else {
m_aPtCtrl[i] = ( 1 - dUTrim) * m_aWeCtrl[i] * m_aPtCtrl[i] + dUTrim * m_aWeCtrl[i+1] * m_aPtCtrl[i+1] ;
m_aWeCtrl[i] = ( 1 - dUTrim) * m_aWeCtrl[i] + dUTrim * m_aWeCtrl[i+1] ;
m_aPtCtrl[i] = m_aPtCtrl[i] * ( 1 / m_aWeCtrl[i]) ;
}
}
}
// con i controlli sopra fatti rimane validata, ma la grafica va ricalcolata
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::TrimEndAtParam( double dUTrim)
{
// la curva deve essere valida
if ( m_nStatus != OK)
return false ;
// se devo togliere a partire dall'inizio o prima non rimane niente
if ( dUTrim <= 0)
return false ;
// se devo togliere dalla fine o dopo non tolgo niente
if ( dUTrim >= 1)
return true ;
// eseguo il trim (algoritmo di de Casteljau)
for ( int k = 1 ; k <= m_nDeg ; k ++) {
for ( int i = m_nDeg ; i >= k ; i --) {
if ( ! m_bRat) {
m_aPtCtrl[i] = ( 1 - dUTrim) * m_aPtCtrl[i-1] + dUTrim * m_aPtCtrl[ i] ;
}
else {
m_aPtCtrl[i] = ( 1 - dUTrim) * m_aWeCtrl[i-1] * m_aPtCtrl[i-1] + dUTrim * m_aWeCtrl[ i] * m_aPtCtrl[ i] ;
m_aWeCtrl[i] = ( 1 - dUTrim) * m_aWeCtrl[i-1] + dUTrim * m_aWeCtrl[ i] ;
m_aPtCtrl[i] = m_aPtCtrl[i] * ( 1 / m_aWeCtrl[i]) ;
}
}
}
// con i controlli sopra fatti rimane validata, ma la grafica va ricalcolata
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::TrimStartAtLen( double dLenTrim)
{
// lunghezze negative vengono considerate nulle
dLenTrim = __max( dLenTrim, 0) ;
// converto le lunghezze in valori parametrici
double dUTrim ;
if ( ! GetParamAtLength( dLenTrim, dUTrim))
return false ;
// utilizzo il trim sui parametri
return TrimStartAtParam( dUTrim) ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::TrimEndAtLen( double dLenTrim)
{
// lunghezze negative vengono considerate nulle
dLenTrim = __max( dLenTrim, 0) ;
// converto le lunghezze in valori parametrici
double dUTrim ;
if ( ! GetParamAtLength( dLenTrim, dUTrim))
return false ;
// utilizzo il trim sui parametri
return TrimEndAtParam( dUTrim) ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Translate( const Vector3d& vtMove)
{
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// traslo i punti di controllo
for ( int i = 0 ; i <= m_nDeg ; ++ i)
m_aPtCtrl[i].Translate( vtMove) ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dCosAng, double dSinAng)
{
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// verifico validità dell'asse di rotazione
if ( vtAx.IsSmall())
return false ;
// ruoto i punti di controllo
for ( int i = 0 ; i <= m_nDeg ; ++ i)
m_aPtCtrl[i].Rotate( ptAx, vtAx, dCosAng, dSinAng) ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double dCoeffZ)
{
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// verifico non sia nulla
if ( fabs( dCoeffX) < EPS_ZERO && fabs( dCoeffY) < EPS_ZERO && fabs( dCoeffZ) < EPS_ZERO)
return false ;
// scalo i punti di controllo
for ( int i = 0 ; i <= m_nDeg ; ++ i)
m_aPtCtrl[i].Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::Mirror( const Point3d& ptOn, const Vector3d& vtNorm)
{
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// verifico validità del piano di specchiatura
if ( vtNorm.IsSmall())
return false ;
// specchio i punti di controllo
for ( int i = 0 ; i <= m_nDeg ; ++ i)
m_aPtCtrl[i].Mirror( ptOn, vtNorm) ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::ToGlob( const Frame3d& frRef)
{
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// verifico validità del frame
if ( frRef.GetType() == Frame3d::ERR)
return false ;
// trasformo i punti di controllo
for ( int i = 0 ; i <= m_nDeg ; ++ i)
m_aPtCtrl[i].ToGlob( frRef) ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveBezier::ToLoc( const Frame3d& frRef)
{
// la curva deve essere validata
if ( m_nStatus != OK)
return false ;
// verifico validità del frame
if ( frRef.GetType() == Frame3d::ERR)
return false ;
// trasformo i punti di controllo
for ( int i = 0 ; i <= m_nDeg ; ++ i)
m_aPtCtrl[i].ToLoc( frRef) ;
// imposto ricalcolo della grafica
m_OGrMgr.Reset() ;
return true ;
}
Reference in New Issue View Git Blame Copy Permalink